- Author:
- Alan Garny <agarny@hellix.com>
- Date:
- 2021-06-21 13:19:38+12:00
- Desc:
- Python scripts for figures 3 to 7: various cleaning up and file renaming.
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/648/rawfile/ec3c101d5695fb4bfa95a972a40384168ca1f43a/Figure5.py
# To reproduce Figure 3 in the associated Physiome paper,
# execute this script from the command line:
#
# cd [PathToThisFile]
# [PathToOpenCOR]/pythonshell Figure5.py
import matplotlib
matplotlib.use('agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import opencor as opencor
# different values for y_shift
y_shift = [-15, -10, -5, 0, 5, 10, 15]
t = ["time"]
V_m = {}
# load the reference model
simulation = opencor.open_simulation("HumanSAN_Fabbri_Fantini_Wilders_Severi_2017.sedml")
data = simulation.data()
data.set_ending_point(1.9)
data.set_point_interval(0.001)
for y in y_shift:
# reset everything in case we are running interactively and have existing results
simulation.reset(True)
simulation.clear_results()
data.constants()["i_f/i_f_y_gate/y_shift"] = y
simulation.run()
ds = simulation.results().data_store()
V_m[y] = ds.voi_and_variables()["Membrane/V"].values()
simulation.reset(True)
Time = {}
simulation.run()
ds = simulation.results().data_store()
Time[t[0]] = ds.voi_and_variables()["environment/time"].values()
V_m.update(Time)
cl = []
for i in range(0, 7):
cl.append(((np.where(V_m[y_shift[i]] == V_m[y_shift[i]][700:1500].min()))[0]) -
(np.where(V_m[y_shift[i]] == V_m[y_shift[i]][200:700].min()))[0])
y_shift = [-15, -10, -5, 0, 5, 10, 15]
plt.figure(figsize=(14, 14))
plt.subplot(2, 2, 1)
plt.plot(y_shift, cl, 'navy', linestyle='', marker='D', markersize='14', label='', linewidth=3)
plt.grid()
plt.ylim(400, 1200)
plt.yticks(np.arange(400, 1300, 200))
plt.xlabel('y$_{\infty}$ shift', fontsize=16)
plt.tick_params(axis='both', labelsize=14)
plt.ylabel('CL (ms)', fontsize=16)
plt.title('A', loc='left', y=1.05, x=-0.06, fontsize='20')
DDR = []
for i in range(0, 7):
A = V_m[y_shift[i]][200:700].min()
end = (((np.where(V_m[y_shift[i]] == V_m[y_shift[i]][200:700].min()))[0]) + 100)
B = V_m[y_shift[i]][end]
DDR.append((B - A) * 10)
plt.subplot(2, 2, 2)
plt.plot(y_shift, DDR, 'navy', linestyle='', marker='D', markersize='14', label='', linewidth=3)
plt.grid()
plt.ylim(20, 80)
plt.yticks(np.arange(20, 90, 20))
plt.xlabel('y$_{\infty}$ shift', fontsize='16')
plt.ylabel('DDR$_{100}$ (mV/s)', fontsize='16')
plt.tick_params(axis='both', labelsize=14)
plt.title('B', loc='left', y=1.05, x=-0.06, fontsize='20')
MDP = []
for i in range(0, 7):
A = V_m[y_shift[i]][200:600].min()
MDP.append(A)
plt.subplot(2, 2, 3)
plt.plot(y_shift, MDP, 'navy', linestyle='', marker='D', markersize='14', label='')
plt.grid()
plt.ylim(-65, -50)
plt.yticks(np.arange(-65, -45, 5))
plt.xlabel('y$_{\infty}$ shift', fontsize='16')
plt.ylabel('MDP (mV)', fontsize='16')
plt.tick_params(axis='both', labelsize=14)
plt.title('C', loc='left', y=1.05, x=-0.06, fontsize='20')
def get_range(minx, maxx):
return 0.9 * (maxx - minx)
def find_min_index(data):
dy = np.gradient(data)
df = np.diff(dy)
df_max = np.where(df == max(df))
return df_max[-1][-1]
def find_max_index(data):
return np.where(data == min(data))[-1][-1]
def plot_apd90(apd):
plt.plot(y_shift, list(apd.values()), 'navy', linestyle='', marker='D', markersize='14', label='')
plt.grid()
plt.ylim(150, 180)
plt.yticks(np.arange(150, 190, 10))
plt.xlabel('y$_{\infty}$ shift', fontsize='16')
plt.ylabel('APD$_{90}$ (ms)', fontsize='16')
plt.tick_params(axis='both', labelsize=14)
plt.title('D', loc='left', y=1.05, x=-0.06, fontsize='20')
output = dict()
y_shift_list = [-15, -10, -5, 0, 5, 10, 15]
# find the start value
for i in range(0, 7):
d = V_m[y_shift[i]][200:700]
min_index = find_min_index(d)
max_index = find_max_index(d)
apd90 = get_range(min_index, max_index)
output[y_shift_list[i]] = apd90
plt.subplot(2, 2, 4)
plot_apd90(output)
plt.tight_layout(pad=0.5, w_pad=3, h_pad=3)
plt.savefig('Figure5.png')
