- Author:
- Shelley Fong <sfon036@UoA.auckland.ac.nz>
- Date:
- 2024-11-12 14:16:20+13:00
- Desc:
- init
- Permanent Source URI:
- http://models.cellml.org/workspace/c17/rawfile/d44e06713a1dab8aa3c903ad6a6582e69adab1e1/fig5_Kr.py
import os
import matplotlib
matplotlib.use('agg')
import numpy as np
import matplotlib.pyplot as plt
import opencor as opencor
import timeit
start = timeit.default_timer()
def load_sedml(filename):
return opencor.open_simulation(filename)
start, stop, step = -150, 50, 5
Vm_plot = np.arange(start, stop + step, step)
Vm_plot = Vm_plot.reshape((1, Vm_plot.shape[0]))
file = load_sedml("Current_Ikr.sedml")
def run_sim1(Vm_step, protocol, Xr1_0, Xr1_1, Xr1_2, Xr1_5, Xr1_6, Xr2_1, Xr2_2, Xr2_5, Xr2_6, Xr1_7, Xr2_7):
period = [3000, 4000, 1000] #[30, 40, 10]
step1 = -40
stepf = -40
if protocol == 2: #Es-Salah_lamoureaux
period = [3000, 4000, 1000]
step1 = -50
stepf = -50
elif protocol == 3: #Wu
period = [3000, 2500, 1000]
step1 = -50
stepf = -40
col1 = step1 * np.ones(Vm_step.T.shape)
col2 = Vm_step.T
col3 = stepf * np.ones(Vm_step.T.shape)
Vm_mat = np.hstack((col1, col2, col3))
Y = np.array((-40, 0, 1)).reshape(3, 1)
IKr_tail = np.zeros(Vm_step.shape)
time_tail = np.zeros(Vm_step.shape)
# sedml file
data = file.data()
Ikr_IV = []
x_act = np.zeros(Vm_step.T.shape)
x_inact = np.zeros(Vm_step.T.shape)
tau_act = np.zeros(Vm_step.T.shape)
tau_inact = np.zeros(Vm_step.T.shape)
for j in range(len(Vm_step.T)):
print(j)
file.reset(True)
file.clear_results()
if protocol == 1:
data.constants()["Current_Ikr/Ki"] = 150 # Ma
data.constants()["Current_Ikr/Ko"] = 5.4
elif protocol == 2:
data.constants()["Current_Ikr/Ki"] = 145 # Es-Salah_lamoureaux
data.constants()["Current_Ikr/Ko"] = 4
elif protocol == 3:
data.constants()["Current_Ikr/Ki"] = 150 # WU
data.constants()["Current_Ikr/Ko"] = 5.4
elif protocol == 4:
data.constants()["Current_Ikr/Ki"] = 145 # Bellin
data.constants()["Current_Ikr/Ko"] = 5.4
data.constants()["Current_Ikr/Xr1_0"] = Xr1_0
data.constants()["Current_Ikr/Xr1_1"] = Xr1_1
data.constants()["Current_Ikr/Xr1_2"] = Xr1_2
data.constants()["Current_Ikr/Xr1_5"] = Xr1_5
data.constants()["Current_Ikr/Xr1_6"] = Xr1_6
data.constants()["Current_Ikr/Xr2_1"] = Xr2_1
data.constants()["Current_Ikr/Xr2_2"] = Xr2_2
data.constants()["Current_Ikr/Xr2_5"] = Xr2_5
data.constants()["Current_Ikr/Xr2_6"] = Xr2_6
data.constants()["Current_Ikr/Xr1_7"] = Xr1_7
data.constants()["Current_Ikr/Xr2_7"] = Xr2_7
# for v in Vm_step:
data.states()["Current_Ikr/v"] = Vm_step.T[j]
data.set_starting_point(0)
data.set_ending_point(10)
data.set_point_interval(10)
file.run()
ds = file.results().data_store()
x_inf_act = ds.voi_and_variables()["Current_Ikr/x_inf_act"].values()
x_inf_inact = ds.voi_and_variables()["Current_Ikr/x_inf_inact"].values()
tau_kr_act = ds.voi_and_variables()["Current_Ikr/tau_kr_act"].values()
tau_kr_inact = ds.voi_and_variables()["Current_Ikr/tau_kr_inact"].values()
x_act[j] = x_inf_act[-1]
x_inact[j] = x_inf_inact[-1]
tau_act[j] = tau_kr_act[-1]
tau_inact[j] = tau_kr_inact[-1]
Y_init = Y.flatten()
tail_max = 0
for i in range(len(period)):
Y_init[0] = Vm_mat[j, i].T
# if i == 0:
# Time_forVmi = 0
# timespan = [sum(period[0:i]), Time_forVmi + period[i]]
# else:
timespan = [sum(period[0:(i)]), sum(period[0:i+1])]
# file = load_sedml(filename)
file.reset(True)
file.clear_results()
if protocol == 1:
data.constants()["Current_Ikr/Ki"] = 150 #Ma
data.constants()["Current_Ikr/Ko"] = 5.4
elif protocol == 2:
data.constants()["Current_Ikr/Ki"] = 145 #Es-Salah_lamoureaux
data.constants()["Current_Ikr/Ko"] = 4
elif protocol == 3:
data.constants()["Current_Ikr/Ki"] = 150 #WU
data.constants()["Current_Ikr/Ko"] = 5.4
elif protocol == 4:
data.constants()["Current_Ikr/Ki"] = 145 #Bellin
data.constants()["Current_Ikr/Ko"] = 5.4
data.constants()["Current_Ikr/Xr1_0"] = Xr1_0
data.constants()["Current_Ikr/Xr1_1"] = Xr1_1
data.constants()["Current_Ikr/Xr1_2"] = Xr1_2
data.constants()["Current_Ikr/Xr1_5"] = Xr1_5
data.constants()["Current_Ikr/Xr1_6"] = Xr1_6
data.constants()["Current_Ikr/Xr2_1"] = Xr2_1
data.constants()["Current_Ikr/Xr2_2"] = Xr2_2
data.constants()["Current_Ikr/Xr2_5"] = Xr2_5
data.constants()["Current_Ikr/Xr2_6"] = Xr2_6
data.constants()["Current_Ikr/Xr1_7"] = Xr1_7
data.constants()["Current_Ikr/Xr2_7"] = Xr2_7
data.states()["Current_Ikr/X_kr_act"] = Y_init[2]
data.states()["Current_Ikr/X_kr_inact"] = Y_init[1]
# data.constants()["main/k_new"] = k_new
data.states()["Current_Ikr/v"] = Y_init[0]
data.set_starting_point(timespan[0])
data.set_ending_point(timespan[1])
data.set_point_interval(1)
file.run()
ds = file.results().data_store()
dX_kr_act = ds.voi_and_variables()["Current_Ikr/X_kr_act"].values()
dX_kr_inact = ds.voi_and_variables()["Current_Ikr/X_kr_inact"].values()
# i_to = ds.voi_and_variables()["Current_Ikr/i_to"].values()
dv = ds.voi_and_variables()["Current_Ikr/v"].values()
Time_hold = (ds.voi_and_variables()["Current_Ikr/t"].values()).T
Time_hold = Time_hold.reshape(Time_hold.shape[0], 1)
results = np.stack((dv, dX_kr_inact, dX_kr_act))
values_hold = results.T
Y_init = values_hold[-1, :]
IKr_hold = np.zeros(Time_hold.shape)
for k in range(0,60):
timespan1 = [0, 60] #60
# file.reset(True)
# file.clear_results()
if protocol == 1:
data.constants()["Current_Ikr/Ki"] = 150 # Ma
data.constants()["Current_Ikr/Ko"] = 5.4
elif protocol == 2:
data.constants()["Current_Ikr/Ki"] = 145 # Es-Salah_lamoureaux
data.constants()["Current_Ikr/Ko"] = 4
elif protocol == 3:
data.constants()["Current_Ikr/Ki"] = 150 # WU
data.constants()["Current_Ikr/Ko"] = 5.4
elif protocol == 4:
data.constants()["Current_Ikr/Ki"] = 145 # Bellin
data.constants()["Current_Ikr/Ko"] = 5.4
data.constants()["Current_Ikr/Xr1_0"] = Xr1_0
data.constants()["Current_Ikr/Xr1_1"] = Xr1_1
data.constants()["Current_Ikr/Xr1_2"] = Xr1_2
data.constants()["Current_Ikr/Xr1_5"] = Xr1_5
data.constants()["Current_Ikr/Xr1_6"] = Xr1_6
data.constants()["Current_Ikr/Xr2_1"] = Xr2_1
data.constants()["Current_Ikr/Xr2_2"] = Xr2_2
data.constants()["Current_Ikr/Xr2_5"] = Xr2_5
data.constants()["Current_Ikr/Xr2_6"] = Xr2_6
data.constants()["Current_Ikr/Xr1_7"] = Xr1_7
data.constants()["Current_Ikr/Xr2_7"] = Xr2_7
data.states()["Current_Ikr/X_kr_act"] = values_hold[k][2]
data.states()["Current_Ikr/X_kr_inact"] = values_hold[k][1]
# data.constants()["main/k_new"] = k_new
data.states()["Current_Ikr/v"] = values_hold[k][0]
data.set_starting_point(timespan1[0])
data.set_ending_point(timespan1[1])
data.set_point_interval(60)
file.run()
ds = file.results().data_store()
dX_kr_act = ds.voi_and_variables()["Current_Ikr/X_kr_act"].values()
dX_kr_inact = ds.voi_and_variables()["Current_Ikr/X_kr_inact"].values()
i_Kr = ds.voi_and_variables()["Current_Ikr/i_Kr"].values()
IKr_hold[k] = i_Kr[0]
if i == 2 and Time_hold[k] < sum(period[0:2])+10:
if abs(IKr_hold[k]) > tail_max:
print(j,i,k)
time_tail[0][j] = Time_hold[k]
IKr_tail[0][j] = IKr_hold[k]
tail_max = abs(IKr_hold[k])
if i == 0:
Time_j = Time_hold
Ikr_j = IKr_hold
else:
Time_j = np.vstack((Time_j, Time_hold))
Ikr_j = np.vstack((Ikr_j, IKr_hold))
return IKr_tail, x_act, x_inact, tau_act, tau_inact
if __name__ == '__main__':
plt.figure(figsize=(16, 13))
Xr1_0_list = [0.2346, 0.1173, 0.1346, 0.3856, 0.2180]
Xr1_1_list = [0.0024, 0.0051, 0.0028, 0.0126, 0.0057]
Xr1_2_list = [12.1564, 9.9823, 16.4121, 15.9432, 13.6235]
Xr1_5_list = [0.1168, 0.0192, 0.044, 0.0105, 0.0476]
Xr1_6_list = [-7.1218, -5.5038, -7.8373, -7.8094, -7.0681]
Xr2_1_list = [0.0125, 0.0125, 0.0125, 0.0125, 0.0125]
Xr2_2_list = [-25.9945, -25.9945, -25.9945, -25.9945, -25.9945]
Xr2_5_list = [37.3426, 37.3426, 37.3426, 37.3426, 37.3426]
Xr2_6_list = [22.0920, 22.0920, 22.0920, 22.0920, 22.0920]
Xr1_7_list = [50, 50, 50, 50, 50]
Xr2_7_list = [0, 0, 0, 0, 0]
plt.subplot(2, 2, 1)
# protocol = 3
# Ikr_IV1 = run_sim1(Vm_plot, protocol, Xr1_0_list[0], Xr1_1_list[0], Xr1_2_list[0], Xr1_5_list[0], Xr1_6_list[0],
# Xr2_1_list[0],
# Xr2_2_list[0], Xr2_5_list[0], Xr2_6_list[0], Xr1_7_list[0], Xr2_7_list[0])
# plt.plot(Vm_plot.T, Ikr_IV1[1], color='orangered', label='Wu Lab', linewidth=4)
# plt.plot(Vm_plot.T, Ikr_IV1[2], color='orangered', linewidth=4)
# protocol = 1
# Iks_IV2 = run_sim1(Vm_plot, protocol, Xr1_0_list[1], Xr1_1_list[1], Xr1_2_list[1], Xr1_5_list[1], Xr1_6_list[1],
# Xr2_1_list[1],
# Xr2_2_list[1], Xr2_5_list[1], Xr2_6_list[1], Xr1_7_list[1], Xr2_7_list[1])
# plt.plot(Vm_plot.T, Iks_IV2[1], color='blue', label='Ma et al.', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV2[2], color='blue', linewidth=4)
# protocol = 2
# Iks_IV3 = run_sim1(Vm_plot, protocol, Xr1_0_list[2], Xr1_1_list[2], Xr1_2_list[2], Xr1_5_list[2], Xr1_6_list[2],
# Xr2_1_list[2],
# Xr2_2_list[2], Xr2_5_list[2], Xr2_6_list[2], Xr1_7_list[2], Xr2_7_list[2])
# plt.plot(Vm_plot.T, Iks_IV3[1], color='purple', label='Es Salah', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV3[2], color='purple', linewidth=4)
# protocol = 4
# Iks_IV4 = run_sim1(Vm_plot, protocol, Xr1_0_list[3], Xr1_1_list[3], Xr1_2_list[3], Xr1_5_list[3], Xr1_6_list[3],
# Xr2_1_list[3],
# Xr2_2_list[3], Xr2_5_list[3], Xr2_6_list[3], Xr1_7_list[3], Xr2_7_list[3])
# plt.plot(Vm_plot.T, Iks_IV4[1], color='green', label='Bellin et al.', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV4[2], color='green', linewidth=4)
protocol = 1
Iks_IV5 = run_sim1(Vm_plot, protocol, Xr1_0_list[4], Xr1_1_list[4], Xr1_2_list[4], Xr1_5_list[4], Xr1_6_list[4], Xr2_1_list[4], Xr2_2_list[4], Xr2_5_list[4], Xr2_6_list[4], Xr1_7_list[4], Xr2_7_list[4])
plt.plot(Vm_plot.T, Iks_IV5[1], color='black', label='Baseline Model', linewidth=4)
plt.plot(Vm_plot.T, Iks_IV5[2], color='black', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize=18)
plt.ylabel('Normalized I$_{Kr}$', fontsize=18)
plt.tick_params(axis='both', labelsize='18')
plt.xticks(np.arange(-150, 51, 50))
plt.yticks(np.arange(0, 1.1, 0.5))
plt.xlim(-150, 50)
plt.ylim(0, 1)
plt.title('A', fontsize=18)
# plt.legend(fontsize='15', loc='best')
plt.subplot(2, 2, 2)
# plt.plot(Vm_plot.T, Ikr_IV1[0].T, color='orangered', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV2[0].T, color='blue', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV3[0].T, color='purple', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV4[0].T, color='green', linewidth=4)
plt.plot(Vm_plot.T, Iks_IV5[0].T, color='black', linewidth=4)
plt.xlabel('Voltage (mV)', size='18')
plt.ylabel('I$_{kr}$ (pA/pF)', size='18')
plt.tick_params(axis='both', labelsize='18')
plt.xticks(np.arange(-150, 51, 50))
plt.yticks(np.arange(0, 3.1, 1))
plt.xlim(-150, 50)
plt.ylim(-0.2, 3)
plt.title('B', fontsize=18)
plt.subplot(2, 2, 3)
# plt.plot(Vm_plot.T, Ikr_IV1[3], color='orangered', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV2[3], color='blue', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV3[3], color='purple', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV4[3], color='green', linewidth=4)
plt.plot(Vm_plot.T, Iks_IV5[3], color='black', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize=18)
plt.ylabel('Tau$_{act}$ (ms)', fontsize=18)
plt.tick_params(axis='both', labelsize='18')
plt.xticks(np.arange(-150, 51, 50))
plt.yticks(np.arange(0, 1001, 500))
plt.xlim(-150, 50)
plt.ylim(0, 1000)
plt.title('C', fontsize='18')
plt.subplot(2, 2, 4)
# plt.plot(Vm_plot.T, Ikr_IV1[4], color='orangered', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV2[4], color='blue', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV3[4], color='purple', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV4[4], color='green', linewidth=4)
plt.plot(Vm_plot.T, Iks_IV5[4], color='black', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize=18)
plt.ylabel('Tau$_{inact}$ (ms)', fontsize=18)
plt.tick_params(axis='both', labelsize='18')
plt.xticks(np.arange(-150, 51, 50))
plt.yticks(np.arange(0, 3.1, 1))
plt.xlim(-150, 50)
plt.ylim(0, 3)
plt.title('D', fontsize=18)
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.9,
wspace=0.3,
hspace=0.3)
stop = timeit.default_timer()
print('Time: ', stop - start)
# plt.show()
plt.savefig('figures/Figure5_Kr.png')
