- Author:
- Shelley Fong <sfon036@UoA.auckland.ac.nz>
- Date:
- 2024-09-27 09:52:37+12:00
- Desc:
- init
- Permanent Source URI:
- https://models.cellml.org/workspace/702/rawfile/841bbcd6641a430406471e20e4c7466aaea97112/fig4_CaL.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.1, 50.1, 2#2
Vm_plot = np.arange(start, stop + step, step)
Vm_plot = Vm_plot.reshape((1, Vm_plot.shape[0]))
file = load_sedml("Current_Ica.sedml")
def run_sim1(Vm_step, protocol, d_0, d_1, d_2, d_5, d_6, f_1, f_2, f_5, f_6, tau_d_const, tau_f_const):
Y = np.array((-70, 0.3, 0, 0, 1, 1)).reshape(6, 1)
period = [3000, 4000, 1000]
step1 = -40
stepf = -40
if protocol == 1: #Ma
period = [1000, 3000, 100, 100]
step1 = -80 #mV
step2 = -50
stepf = -50 #mV
elif protocol == 2: #Es-Salah_lamoureaux
period = [1000, 3000, 125, 100]
step1 = -80 # mV
step2 = -50
stepf = -50
else: #Veerman
period = [500, 500, 300, 100]
step1 = -80 # mV
step2 = -80
stepf = -80
col1 = step1 * np.ones(Vm_step.T.shape)
col2 = step2 * np.ones(Vm_step.T.shape)
col3 = Vm_step.T
col4 = stepf * np.ones(Vm_step.T.shape)
Vm_mat = np.hstack((col1, col2, col3, col4))
Time = sum(period) * np.ones((sum(period) * 100, len(Vm_step.T)))
ICaL = 0.01 * np.ones((sum(period) * 100, len(Vm_step.T)))
IV_cal = np.zeros(Vm_step.T.shape)
# sedml file
data = file.data()
x_act = np.zeros(Vm_step.T.shape)
x_inact = np.zeros(Vm_step.T.shape)
x_fca_inact = np.zeros(Vm_step.T.shape)
tau_act = np.zeros(Vm_step.T.shape)
tau_inact = np.zeros(Vm_step.T.shape)
tau_fca = 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_IcaL/Ki"] = 0 # Ma
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 5
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 5
elif protocol == 2:
data.constants()["Current_IcaL/Ki"] = 0 # Es Salah
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 5
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 5
elif protocol == 4:
data.constants()["Current_IcaL/Ki"] = 0
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 0
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 2.5
else:
data.constants()["Current_IcaL/Ki"] = 0 # Veerman
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 0
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 1.8
data.constants()["Current_IcaL/d_0"] = d_0
data.constants()["Current_IcaL/d_1"] = d_1
data.constants()["Current_IcaL/d_2"] = d_2
data.constants()["Current_IcaL/d_5"] = d_5
data.constants()["Current_IcaL/d_6"] = d_6
data.constants()["Current_IcaL/f_1"] = f_1
data.constants()["Current_IcaL/f_2"] = f_2
data.constants()["Current_IcaL/f_5"] = f_5
data.constants()["Current_IcaL/f_6"] = f_6
data.constants()["Current_IcaL/tau_d_const"] = tau_d_const
data.constants()["Current_IcaL/tau_f_const"] = tau_f_const
# for v in Vm_step:
data.states()["Current_IcaL/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_IcaL/x_ca_inf_act"].values()
x_inf_inact = ds.voi_and_variables()["Current_IcaL/x_ca_inf_inact"].values()
x_fca_inf_inact = ds.voi_and_variables()["Current_IcaL/x_fca_inf_inact"].values()
tau_ca_act = ds.voi_and_variables()["Current_IcaL/tau_ca_act"].values()
tau_ca_inact = ds.voi_and_variables()["Current_IcaL/tau_ca_inact"].values()
tau_fca_inact = ds.voi_and_variables()["Current_IcaL/tau_fca"].values()
x_act[j] = x_inf_act[-1]
x_inact[j] = x_inf_inact[-1]
x_fca_inact[j] = x_fca_inf_inact[-1]
tau_act[j] = tau_ca_act[-1]
tau_inact[j] = tau_ca_inact[-1]
tau_fca[j] = tau_fca_inact[-1]
Y_init = Y.flatten()
for i in range(len(period)):
Y_init[0] = Vm_mat[j, i].T
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_IcaL/Ki"] = 0 #Ma
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 5
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 5
elif protocol == 2:
data.constants()["Current_IcaL/Ki"] = 0 #Es Salah
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 5
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 5
elif protocol == 4:
data.constants()["Current_IcaL/Ki"] = 0
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 0
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 2.5
else:
data.constants()["Current_IcaL/Ki"] = 0 #Veerman
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 0
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 1.8
data.constants()["Current_IcaL/d_0"] = d_0
data.constants()["Current_IcaL/d_1"] = d_1
data.constants()["Current_IcaL/d_2"] = d_2
data.constants()["Current_IcaL/d_5"] = d_5
data.constants()["Current_IcaL/d_6"] = d_6
data.constants()["Current_IcaL/f_1"] = f_1
data.constants()["Current_IcaL/f_2"] = f_2
data.constants()["Current_IcaL/f_5"] = f_5
data.constants()["Current_IcaL/f_6"] = f_6
data.constants()["Current_IcaL/tau_d_const"] = tau_d_const
data.constants()["Current_IcaL/tau_f_const"] = tau_f_const
# data.constants()["main/k_new"] = k_new
data.states()["Current_IcaL/v"] = Y_init[0]
data.states()["Current_IcaL/ca_SR"] = Y_init[1]
data.states()["Current_IcaL/ca_i"] = Y_init[2]
data.states()["Current_IcaL/X_ca_act"] = Y_init[3]
data.states()["Current_IcaL/X_ca_inact"] = Y_init[4]
data.states()["Current_IcaL/X_fca_inact"] = Y_init[5]
data.set_starting_point(timespan[0])
data.set_ending_point(timespan[1])
data.set_point_interval(1)
file.run()
ds = file.results().data_store()
dv = ds.voi_and_variables()["Current_IcaL/v"].values()
dca_SR = ds.voi_and_variables()["Current_IcaL/ca_SR"].values()
dca_i = ds.voi_and_variables()["Current_IcaL/ca_i"].values()
dX_ca_act = ds.voi_and_variables()["Current_IcaL/X_ca_act"].values()
dX_ca_inact = ds.voi_and_variables()["Current_IcaL/X_ca_inact"].values()
dX_fca_inact = ds.voi_and_variables()["Current_IcaL/X_fca_inact"].values()
Time_hold = (ds.voi_and_variables()["Current_IcaL/t"].values()).T
Time_hold = Time_hold.reshape(Time_hold.shape[0], 1)
results = np.stack((dv, dca_SR, dca_i, dX_ca_act, dX_ca_inact, dX_fca_inact))
values_hold = results.T
Y_init = values_hold[-1, :]
ICaL_hold = np.zeros(Time_hold.shape)
for k in range(0,80):
timespan1 = [0,80]
# file.reset(True)
# file.clear_results()
if protocol == 1:
data.constants()["Current_IcaL/Ki"] = 0 # Ma
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 5
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 5
elif protocol == 2:
data.constants()["Current_IcaL/Ki"] = 0 # Es Salah
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 5
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 5
elif protocol == 4:
data.constants()["Current_IcaL/Ki"] = 0
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 0
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 2.5
else:
data.constants()["Current_IcaL/Ki"] = 0 # Veerman
data.constants()["Current_IcaL/Ko"] = 0
data.constants()["Current_IcaL/Na_i"] = 0
data.constants()["Current_IcaL/Nao"] = 0
data.constants()["Current_IcaL/Cao"] = 1.8
data.constants()["Current_IcaL/d_0"] = d_0
data.constants()["Current_IcaL/d_1"] = d_1
data.constants()["Current_IcaL/d_2"] = d_2
data.constants()["Current_IcaL/d_5"] = d_5
data.constants()["Current_IcaL/d_6"] = d_6
data.constants()["Current_IcaL/f_1"] = f_1
data.constants()["Current_IcaL/f_2"] = f_2
data.constants()["Current_IcaL/f_5"] = f_5
data.constants()["Current_IcaL/f_6"] = f_6
data.constants()["Current_IcaL/tau_d_const"] = tau_d_const
data.constants()["Current_IcaL/tau_f_const"] = tau_f_const
# data.constants()["main/k_new"] = k_new
data.states()["Current_IcaL/v"] = values_hold[k][0]
data.states()["Current_IcaL/ca_SR"] = values_hold[k][1]
data.states()["Current_IcaL/ca_i"] = values_hold[k][2]
data.states()["Current_IcaL/X_ca_act"] = values_hold[k][3]
data.states()["Current_IcaL/X_ca_inact"] = values_hold[k][4]
data.states()["Current_IcaL/X_fca_inact"] = values_hold[k][5]
data.set_starting_point(timespan1[0])
data.set_ending_point(timespan1[1])
data.set_point_interval(80)
file.run()
ds = file.results().data_store()
i_CaL = ds.voi_and_variables()["Current_IcaL/i_CaL"].values()
ICaL_hold[k] = i_CaL[0]
if i == 0:
Time_j = Time_hold
i_CaL_j = ICaL_hold
else:
Time_j = np.vstack((Time_j, Time_hold))
i_CaL_j = np.vstack((i_CaL_j, ICaL_hold))
if i == 2:
start = 30
index = np.count_nonzero((np.array(np.max(abs(ICaL_hold[start:-1])))))
IV_cal[j] = ICaL_hold[index+start-1]
return IV_cal, x_act, x_inact, x_fca_inact, tau_act, tau_inact, tau_fca
if __name__ == '__main__':
plt.figure(figsize=(15, 13))
d_0_list = [0.19016464505208927, 0.38010243347999983, 0.484048813248, 0.17779487373692016, 0.30802769137925234]
d_1_list = [7.724718062981113, 14.567438770685602, 8.770648120828913, 15.560795677650413, 12.966294189721642]
d_2_list = [6.388825030309134, 7.501436625497612, 7.8748150063959255, 5.861186262240001, 7.07914596471118]
d_5_list = [0.08803087507239495, 0.032552967698800464, 0.062165144462064205, 0.04001013436000526,
0.044909415506956644]
d_6_list = [-6.265773108811118, -7.305350518436171, - 7.6680726347258314, - 5.756217954563908, - 6.909880369241971]
f_1_list = [0.001644383040759749, 3.454960428030549E-4, 5.925134493322452E-4, 5.997599860362688E-4,
5.125898260571896E-4]
f_2_list = [-82.46034376474563, -43.41688999393955, - 52.77826478681527, - 52.32198132085614, - 49.50571203387032]
f_5_list = [818.2069824654249, 1580.8554165474573, 729.0226806467049, 3483.7555733487948, 1931.2112235143188]
f_6_list = [4.35226624639966, 6.131075895743708, 6.1849082640235915, 4.874098339328655, 5.730027499698651]
tau_d_const_list = [1.736167780874709, 1.6122826510533983, 1.6664595173798702, 1.6959986720577191,
1.6582469468303291]
tau_f_const_list = [97.74237390031138, 99.52320589734498, 97.14763412621564, 104.71683748974891, 100.46255917110318]
plt.subplot(2, 2, 1)
# protocol = 1
# ICaL_IV1 = run_sim1(Vm_plot, protocol, d_0_list[0], d_1_list[0], d_2_list[0], d_5_list[0], d_6_list[0], f_1_list[0],
# f_2_list[0], f_5_list[0], f_6_list[0],
# tau_d_const_list[0], tau_f_const_list[0])
# plt.plot(Vm_plot[0], ICaL_IV1[1], color='blue', label='Ma', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV1[2], color='blue', linewidth=4)
# protocol = 3
# ICaL_IV2 = run_sim1(Vm_plot, protocol, d_0_list[1], d_1_list[1], d_2_list[1], d_5_list[1], d_6_list[1], f_1_list[1],
# f_2_list[1], f_5_list[1], f_6_list[1],
# tau_d_const_list[1], tau_f_const_list[1])
# plt.plot(Vm_plot[0], ICaL_IV2[1], color='orangered', label='Veerman1', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV2[2], color='orangered', linewidth=4)
# protocol = 3
# ICaL_IV3 = run_sim1(Vm_plot, protocol, d_0_list[2], d_1_list[2], d_2_list[2], d_5_list[2], d_6_list[2], f_1_list[2],
# f_2_list[2], f_5_list[2], f_6_list[2],
# tau_d_const_list[2], tau_f_const_list[2])
# plt.plot(Vm_plot[0], ICaL_IV3[1], color='purple', label='Veerman2', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV3[2], color='purple', linewidth=4)
# protocol = 2
# ICaL_IV4 = run_sim1(Vm_plot, protocol, d_0_list[3], d_1_list[3], d_2_list[3], d_5_list[3], d_6_list[3], f_1_list[3],
# f_2_list[3], f_5_list[3], f_6_list[3],
# tau_d_const_list[3], tau_f_const_list[3])
# plt.plot(Vm_plot[0], ICaL_IV4[1], color='green', label='EsSalah', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV4[2], color='green', linewidth=4)
protocol = 4
ICaL_IV5 = run_sim1(Vm_plot, protocol, d_0_list[4], d_1_list[4], d_2_list[4], d_5_list[4], d_6_list[4], f_1_list[4],
f_2_list[4], f_5_list[4], f_6_list[4],
tau_d_const_list[4], tau_f_const_list[4])
plt.plot(Vm_plot[0], ICaL_IV5[1], color='black', label='Baseline', linewidth=4)
plt.plot(Vm_plot[0], ICaL_IV5[2], color='black', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize=18)
plt.ylabel('Normalized I$_{CaL}$', fontsize=18)
plt.tick_params(axis='both', labelsize='18')
plt.xticks(np.arange(-100, 51, 50))
plt.yticks(np.arange(0, 1.1, 0.5))
plt.xlim(-100, 50)
plt.ylim(0, 1)
plt.title('A', fontsize=18)
# plt.legend(fontsize='15')
plt.subplot(2, 2, 2)
# plt.plot(Vm_plot[0], ICaL_IV1[0].T[0], color='blue', label='Ma', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV2[0].T[0], color='orangered', label='Veerman1', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV3[0].T[0], color='purple', label='Veerman2', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV4[0].T[0], color='green', label='EsSalah', linewidth=4)
plt.plot(Vm_plot[0], ICaL_IV5[0].T[0], color='black', label='Baseline', linewidth=4)
plt.xticks(np.arange(-100, 51, 50))
plt.yticks(np.arange(-65, 21, 20))
plt.xlim(-100, 50)
plt.ylim(-70, 20)
plt.tick_params(axis='both', labelsize='18')
plt.xlabel('Voltage (mV)', size='18')
plt.ylabel('I$_{ICaL}$ (pA/pF)', size='18')
plt.title('B', fontsize=18)
plt.subplot(2, 2, 3)
# plt.plot(Vm_plot[0], ICaL_IV1[4], color='blue', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV2[4], color='orangered', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV3[4], color='purple', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV4[4], color='green', linewidth=4)
plt.plot(Vm_plot[0], ICaL_IV5[4], color='black', linewidth=4)
plt.xlabel('Voltage (mV)', fontsize=18)
plt.ylabel('Tau$_{act},I_{CaL}$ (ms)', fontsize=18)
plt.tick_params(axis='both', labelsize='18')
plt.xticks(np.arange(-100, 51, 50))
plt.yticks(np.arange(1.5, 3.6, 0.5))
plt.xlim(-100, 50)
plt.ylim(1.5, 3.5)
plt.title('C', fontsize=18)
plt.subplot(2, 2, 4)
# plt.plot(Vm_plot[0], ICaL_IV1[5], color='blue', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV2[5], color='orangered', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV3[5], color='purple', linewidth=4)
# plt.plot(Vm_plot[0], ICaL_IV4[5], color='green', linewidth=4)
plt.plot(Vm_plot[0], ICaL_IV5[5], color='black', linewidth=4)
# plt.xlabel('Voltage (mV)', fontsize=18)
# plt.ylabel('Tau$_{inact},I_{CaL}$ (ms)', fontsize=18)
# plt.tick_params(axis='both', labelsize='18')
# plt.xticks(np.arange(-100, 51, 50))
# plt.yticks(np.arange(0, 1001, 500))
# plt.xlim(-100, 50)
# plt.ylim(0, 1000)
# plt.title('D', fontsize=18)
plt.tick_params(axis='both', labelsize='18')
plt.xlabel('Voltage (mV)', size='18')
plt.ylabel('Tau$_{inact},I_{CaL}$ (ms)', size='18')
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.savefig('figures/Figure4_CaL.png')
