- Author:
- Shelley Fong <sfon036@UoA.auckland.ac.nz>
- Date:
- 2024-09-26 14:33:29+12:00
- Desc:
- Fixing file content
- Permanent Source URI:
- https://models.cellml.org/workspace/702/rawfile/ecce9c2e3922ee4da420295236eb86e6d827a940/fig7_Ks.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 = -100, 100, 5
Vm_plot = np.arange(start, stop + step, step)
Vm_plot = Vm_plot.reshape((1, Vm_plot.shape[0]))
file = load_sedml("Current_Iks.sedml")
def run_sim1(Vm_step, protocol, ks_0, ks_1, ks_2, ks_5, ks_6, tau_ks_const):
Y = np.array((-70, 0)).reshape(2, 1)
col1 = -40 * np.ones(Vm_step.T.shape)
col2 = Vm_step.T
col3 = -40 * np.ones(Vm_step.T.shape)
period = [2000, 5000, 100]
Vm_mat = np.hstack((col1, col2, col3))
xy = 2 # variables in y
values = np.zeros((sum(period) * 100, len(Vm_step.T) * xy))
Time = sum(period) * np.ones((sum(period) * 100, len(Vm_step.T)))
IV_place = np.zeros(Vm_step.T.shape)
x_act = np.zeros(Vm_step.T.shape)
tau_act = np.zeros(Vm_step.T.shape)
# sedml file
data = file.data()
for j in range(len(Vm_step.T)):
print(j)
Y_init = Y.flatten()
for i in range(len(period)):
Y_init[0] = Vm_mat[j, i].T
if i == 0:
Time_forVmi = 0
timespan = [Time_forVmi, Time_forVmi + period[i]]
else:
timespan = [Time_forVmi[-1], Time_forVmi[-1] + period[i]]
# file = load_sedml(filename)
file.reset(True)
file.clear_results()
if protocol == 1:
data.constants()["Current_Iks/Ki"] = 130 #Ma Wei
else:
data.constants()["Current_Iks/Ki"] = 145 #Ma
data.constants()["Current_Iks/Ko"] = 5.4
data.constants()["Current_Iks/ks_0"] = ks_0
data.constants()["Current_Iks/ks_1"] = ks_1
data.constants()["Current_Iks/ks_2"] = ks_2
data.constants()["Current_Iks/ks_5"] = ks_5
data.constants()["Current_Iks/ks_6"] = ks_6
data.constants()["Current_Iks/tau_ks_const"] = tau_ks_const
data.states()["Current_Iks/X_ks_act"] = Y[1]
data.states()["Current_Iks/v"] = Y_init[0]
data.set_starting_point(timespan[0])
data.set_ending_point(timespan[1])
data.set_point_interval(0.992063)
file.run()
ds = file.results().data_store()
dX_ks_act = ds.voi_and_variables()["Current_Iks/X_ks_act"].values()
i_Ks = ds.voi_and_variables()["Current_Iks/i_Ks"].values()
dv = ds.voi_and_variables()["Current_Iks/v"].values()
Time_hold = (ds.voi_and_variables()["Current_Iks/t"].values()).T
Time_hold = Time_hold.reshape(Time_hold.shape[0], 1)
results = np.stack((dv, dX_ks_act))
values_hold = results.T
Y_init = results.T[-1, :]
if i == 0:
Time_forVmi = Time_hold
values_forVmi = values_hold
else:
values_forVmi = np.vstack((values_forVmi, values_hold))
Time_forVmi = np.vstack((Time_forVmi, Time_hold))
if i == 1:
timespan1 = [0, 1]
if protocol == 1:
data.constants()["Current_Iks/Ki"] = 140
else:
data.constants()["Current_Iks/Ki"] = 150
data.constants()["Current_Iks/Ko"] = 5.4
data.constants()["Current_Iks/ks_0"] = ks_0
data.constants()["Current_Iks/ks_1"] = ks_1
data.constants()["Current_Iks/ks_2"] = ks_2
data.constants()["Current_Iks/ks_5"] = ks_5
data.constants()["Current_Iks/ks_6"] = ks_6
data.constants()["Current_Iks/tau_ks_const"] = tau_ks_const
data.states()["Current_Iks/X_ks_act"] = values_hold[-1][1]
# data.constants()["main/k_new"] = k_new
data.states()["Current_Iks/v"] = values_hold[-1][0]
data.set_starting_point(timespan1[0])
data.set_ending_point(timespan1[1])
data.set_point_interval(1)
file.run()
ds = file.results().data_store()
x_ks_inf_act = ds.voi_and_variables()["Current_Iks/x_ks_inf_act"].values()
tau_ks_act = ds.voi_and_variables()["Current_Iks/tau_ks_act"].values()
i_ks = ds.voi_and_variables()["Current_Iks/i_Ks"].values()
#
IV_place[j] = i_ks[-1]
x_act[j] = x_ks_inf_act[-1]
tau_act[j] = tau_ks_act[-1]
Time[0: len(Time_forVmi), j] = Time_forVmi[:, 0]
values[0: len(Time_forVmi), ((j * xy) - (xy - 2)): (j * xy) + 2] = values_forVmi
return IV_place, x_act, tau_act
if __name__ == '__main__':
ks_0_list = [0.0030, 0.0096, 0.0105, 0.0077]
ks_1_list = [9.384e-4, 0.0013, 0.0013, 0.0012]
ks_2_list = [180.5160, 1e5, 1e5, 6.6727e4]
ks_5_list = [0.3532, 0.2614, 0.2268, 0.2805]
ks_6_list = [-14.6633, -22.0906, -19.847, -18.8670]
tau_ks_const_list = [1.4223e-5, 2e-11, 2e-11, 4.7411e-6]
plt.figure(figsize=(19, 7))
plt.subplot(1, 3, 1)
# protocol = 0
# Iks_IV1 = run_sim1(Vm_plot, protocol, ks_0_list[0], ks_1_list[0], ks_2_list[0], ks_5_list[0], ks_6_list[0],
# tau_ks_const_list[0])
# x_act = np.array(Iks_IV1[1])
# plt.plot(Vm_plot.T, x_act, color='blue', label= 'Ma et al.', linewidth=4)
# protocol = 1
# Iks_IV2 = run_sim1(Vm_plot, protocol, ks_0_list[1], ks_1_list[1], ks_2_list[1], ks_5_list[1], ks_6_list[1],
# tau_ks_const_list[1])
# x_act = np.array(Iks_IV2[1])
# plt.plot(Vm_plot.T, x_act, color='orangered', label='Ma, Wei et al. patient', linewidth=4)
# protocol = 1
# Iks_IV3 = run_sim1(Vm_plot, protocol, ks_0_list[2], ks_1_list[2], ks_2_list[2], ks_5_list[2], ks_6_list[2],
# tau_ks_const_list[2])
# x_act = np.array(Iks_IV3[1])
# plt.plot(Vm_plot.T, x_act, color='purple', label= 'Ma, Wie et al. iCell', linewidth=4)
protocol = 0
Iks_IV4 = run_sim1(Vm_plot, protocol, ks_0_list[3], ks_1_list[3], ks_2_list[3], ks_5_list[3], ks_6_list[3],
tau_ks_const_list[3])
x_act = np.array(Iks_IV4[1])
plt.plot(Vm_plot.T, x_act, color='black', label= 'Baseline Model', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize= 18)
plt.ylabel('Normalized I$_{Ks}$', fontsize= 18)
plt.xticks(np.arange(-100, 101, 50))
plt.yticks(np.arange(0, 1.1, 0.5))
plt.xlim(-100,100)
plt.ylim(0,1)
plt.tick_params(axis='both', labelsize='18')
plt.title('A', fontsize=18)
# plt.legend(fontsize= '14', loc = 'best')
plt.subplot(1, 3, 2)
# tau_act = np.array(Iks_IV1[2])
# plt.plot(Vm_plot.T, tau_act, color='blue', linewidth=4)
# tau_act = np.array(Iks_IV2[2])
# plt.plot(Vm_plot.T, tau_act, color='orangered', linewidth=4)
# tau_act = np.array(Iks_IV3[2])
# plt.plot(Vm_plot.T, tau_act, color='purple', linewidth=4)
tau_act = np.array(Iks_IV4[2])
plt.plot(Vm_plot.T, tau_act, color='black', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize= 18)
plt.ylabel('Tau$_{act},I_{Ks}$ (ms)', fontsize= 18)
plt.xticks(np.arange(-100, 101, 50))
plt.yticks(np.arange(0, 1001, 500))
plt.xlim(-100,100)
plt.ylim(0,1000)
plt.tick_params(axis= 'both', labelsize= '18')
plt.title('B', fontsize=18)
plt.subplot(1, 3, 3)
# plt.plot(Vm_plot.T, Iks_IV1[0], color='blue', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV2[0], color='orangered', linewidth=4)
# plt.plot(Vm_plot.T, Iks_IV3[0], color='purple', linewidth=4)
plt.plot(Vm_plot.T, Iks_IV4[0], color='black', linewidth=4)
#
plt.xticks(np.arange(-100, 101, 50))
plt.yticks(np.arange(0, 2.1, 0.5))
plt.xlim(-100, 100)
plt.ylim(0, 2)
plt.tick_params(axis='both', labelsize='18')
plt.xlabel('Voltage (mV)', size='18')
plt.ylabel('I$_{ks}$ (pA/pF)', size='18')
plt.title('C', fontsize=18)
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.9,
wspace=0.3,
hspace=0.3)
plt.savefig('figures/Figure7_Ks.png')
stop = timeit.default_timer()
print('Time: ', stop - start)
