ralf/Calcium_Model
1
0
Fork 0
Code Issues 2 Pull Requests Actions Packages Projects Releases Wiki Activity
fb78cc2ddb
Calcium_Model/Model.py

469 lines
16 KiB
Python
Raw Blame History

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
R: float = 8.314 # ideal gas constant (J/(mol*K))
F: float = 96.485 # coulomb_per_millimole, Faraday constant
T: float = 298 # room temperature (K)
RT: float = R * T # J/mol
class Model:
def __init__(self):
self.A_cap: float = 1.534e-4 # cm^2, Capacitive membrane area
self.C_mem: float = 1.0 # uF/cm2 Specific membrane capacitance
self.V_myo: float = 25.84e-6 # uL, Myoplasmic volume
self.current2flux: float = self.A_cap * self.C_mem / 2 / F #
self.period: float = 1000 # ms, pulse period
self.V_mem_rest: float = -80.0 # mV, resting membrane potential
self.Nai: float = 11000 # uM, Myoplasmic Na+ concentration
self.Nao: float = 150000 # uM, Extracellular Na+ concentration
self.eta: float = 0.35 # Controls voltage dependance of Na/Ca2+ excng
self.km_Na: float = (
87500 # uM, Na+ half-saturation constant for Na+/Ca2+ exchange
)
self.k_sat: float = (
0.1 # Na+/Ca2+ exchange saturation factor at very neg potentials
)
self.km_Ca: float = (
1380 # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange
)
self.k_Na_Ca: float = 292.8 # pA/pF, Scaling factor of Na2+/Ca2+ exchange
self.I_max_pCa: float = 1.0 # pA/pF, Maximum Ca2+ pump current
# self.Cai = 0.1 # uM, Cytoplasmic Ca2+ concentration fixed at this point
self.LTRPN_tot: float = (
70.0 # uM, Total myoplasmic troponin low-affinity site concentration
)
self.HTRPN_tot: float = (
140.0 # uM, Total myoplasmic troponin high-affinity site concentration
)
# self.LTRPNCa = 11.2684 # uM, Concentration Ca2+ bound low-affinity troponin-binding sites
# self.HTRPNCa = 125.290 # uM, Concentration Ca2+ bound high-affinity troponin-binding sites
self.k_htrpn_positive: float = (
0.00237 # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affinity sites
)
self.k_htrpn_negative = (
3.2e-5 # ms, Ca2+ off rate const. for troponin high-affinity sites
)
self.k_ltrpn_positive = (
0.0327 # uM^(-1)/ms, Ca2+ on rate const. for troponin low-affinity sites
)
self.k_ltrpn_negative: float = (
0.0196 # ms, Ca2+ off rate const. for troponin low-affinity sites
)
self.CMDN_tot: float = 50 # uM, Total myoplasmic calmodulin concentration
self.Km_CMDN: float = 0.238 # uM, Ca2 half-saturation constant for calmodulin
self.nu_1: float = (
4.5 # 1/ms, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR
)
self.nu_2: float = 1.74e-5 # ms^(-1), Ca2+ leak rate const. from the NSR
self.nu_3: float = 0.45 # uM/ms, SR Ca2+ -ATPase maximum pupmp rate
self.Km_up: float = 0.5 # uM, Half-saturation constant for SR Ca2+ -ATPase pump
self.km_p_ca: float = (
0.5 # uM, Ca2+ half-saturation constant for Ca2+ pump current
)
# self.Ca_NSR = 1299.50 # uM,NSR Ca2+ concentration
self.tau_xfer: float = (
8.0 # ms, Time constant for transfer from subspace to myoplasm
)
self.K_pc_max: float = (
0.23324 # 1/ms, Maximum time constant for Ca2+-induced inactivation
)
self.K_pc_half: float = (
20.0 # uM, Half-saturation constant for Ca2+-induced inactivation
)
self.Kpcb: float = (
0.0005 # 1/ms, Voltage-insensitive rate constant for inactivation
)
self.Gcab: float = 0.000367 # mS/uF
self.Cao: float = 1130.0 # uM, Ca2+ outside the cell
self.gCaL: float = (
0.1729 # mS/uF, Specific maximum conductivity for L-type Ca2+ channel
)
self.ECaL: float = 43.0 # mV, Reversal potential for L-type Ca2
self.V_ss: float = 1.485e-9 # uL, Dyadic aka subspace volume
# self.Ca_JSR = 1299.50 # uM, JSR Ca2+ concentration
self.F_tot: float = 25 # uM, total concentration of Fluo-4
self.k_on: float = 0.1 # 1/uM * 1/ms, Fluo-4 reaction rate constant
self.k_off: float = 0.11 # 1/ms, Fluo-4 reaction rate constant 2
self.V_NSR: float = 2.098e-6 # ul, Network SR volume
self.tau_tr: float = 20 # ms, Time const for transfer from NSR to JSR
self.CSQN_tot: float = (
15000.0 # uM, total junctional SR calsequestrin concentration
)
self.Km_CSQN: float = (
800.0 # uM, Ca2 half-saturation constant for calsequestrin
)
self.k_a_positive: float = 0.006075 # uM^(-4)/ms, RyR Pc1 - Po1 rate constant
self.k_a_negative: float = 0.07125 # 1/ms, RyR Po1 - Pc1 rate constant
self.k_b_positive: float = 0.00405 # uM^(-3)/ms, RyR Po1 - Po2 rate constant
self.k_b_negative: float = 0.965 # 1/ms, RyR Po2 - Po1 rate constant
self.k_c_positive: float = 0.009 # 1/ms, RyR Po1 - Pc2 rate constant
self.k_c_negative: float = 0.0008 # 1/ms, RyR Pc2 - Po1 rate constant
self.n_ryr: int = 4 # RyR Ca2+ cooperativity parameter Pc1 - Po1
self.m_ryr: int = 3 # RyR Ca2+ cooperativity parameter Po1 - Po2
self.I_CaL_max: float = (
7.0 # pA/pF, normalization constant for L-type Ca2+ current
)
self.V_JSR: float = 0.12e-6 # ul, Junctional SR volume
def ode_system(self, t, states: list[float]) -> list[float]:
(
Cass,
Os,
C2,
C3,
C4,
I1,
I2,
I3,
LTRPNCa,
HTRPNCa,
Ca_NSR,
Ca_JSR,
P_RyR,
P_O1,
P_O2,
P_C2,
Cai,
FCa,
) = states
V = self.mem_potential(t)
VF = V * F
Kpcf = 13.0 * (1 - np.exp(-((V + 14.5) ** 2) / 100))
alpha = ( # dimensionless parameter for L-type Ca2+ channel closed states
0.4
* np.exp((V + 12.0) / 10)
* (
1
+ 0.7 * np.exp((-((V + 40.0) ** 2)) / 10.0)
- 0.75 * np.exp(-((V + 20.0) ** 2) / 400.0)
)
) / (1 + 0.12 * np.exp((V + 12.0) / 10.0))
beta = 0.05 * np.exp(
-(V + 12.0) / 13.0
) # parameter for L-type Ca2+ channel closed states
gamma = self.K_pc_max * Cass / (self.K_pc_half + Cass)
C1 = 1 - (Os + C2 + C3 + C4 + I1 + I2 + I3)
dC2dt = 4 * alpha * C1 - beta * C2 + 2 * beta * C3 - 3 * alpha * C2
dC3dt = 3 * alpha * C2 - 2 * beta * C3 + 3 * beta * C4 - 2 * alpha * C3
dC4dt = (
2 * alpha * C3
- 3 * beta * C4
+ 4 * beta * Os
- alpha * C4
+ 0.01 * (4 * self.Kpcb * beta * I1 - alpha * gamma * C4)
+ 0.002 * (4 * beta * I2 - Kpcf * C4)
+ 4 * beta * self.Kpcb * I3
- gamma * Kpcf * C4
)
dI1dt = (
gamma * Os
- self.Kpcb * I1
+ 0.001 * (alpha * I3 - Kpcf * I1)
+ 0.01 * (alpha * gamma * C4 - 4 * beta * self.Kpcb * I1)
)
dI2dt = 0.001 * (Kpcf * Os - alpha * I2)
+self.Kpcb * I3 - gamma * I2 + 0.002 * (Kpcf * C4 - 4 * beta * I2)
dI3dt = (
0.001 * (Kpcf * I1 - alpha * I3)
+ gamma * I2
- self.Kpcb * I3
+ gamma * Kpcf * C4
- 4 * beta * self.Kpcb * I3
)
Bi = 1 / (1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cai) ** 2)
Bss = 1 / (1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2)
J_xfer = (Cass - Cai) / self.tau_xfer
dLTRPNCadt = self.k_ltrpn_positive * Cai * (self.LTRPN_tot - LTRPNCa)
-self.k_ltrpn_negative * (LTRPNCa)
dHTRPNCadt = self.k_htrpn_positive * Cai * (self.HTRPN_tot - HTRPNCa)
-self.k_htrpn_negative * (HTRPNCa)
J_trpn = (
self.k_htrpn_positive * Cai * (self.HTRPN_tot - HTRPNCa)
- self.k_htrpn_negative * HTRPNCa
+ self.k_ltrpn_positive * Cai * (self.LTRPN_tot - LTRPNCa)
- self.k_ltrpn_negative * LTRPNCa
)
J_up = self.nu_3 * (Cai**2 / (self.Km_up**2 - Cai**2))
J_tr = (Ca_NSR - Ca_JSR) / self.tau_tr
J_leak = self.nu_2 * (Ca_NSR - Cai)
dCa_NSRdt = (J_up - J_leak) * (self.V_myo / self.V_NSR)
-J_tr * (self.V_JSR / self.V_NSR)
B_JSR = 1 / (1 + (self.CSQN_tot * self.Km_CSQN) / (self.Km_CSQN + Ca_JSR) ** 2)
J_rel = 0 # self.nu_1*(P_O1+P_O2)*(Ca_JSR-Cass)*P_RyR
dCa_JSRdt = B_JSR * (J_tr - J_rel)
b = 1 / (
(self.km_Na**3 + self.Nao**3)
* (self.km_Ca + self.Cao)
* (1 + self.k_sat * np.exp((self.eta - 1) * VF / RT))
)
p = (
np.exp(self.eta * VF / RT) * self.Nai**3 * self.Cao
- np.exp((self.eta - 1) * VF / RT) * self.Nao**3 * Cai
)
I_NaCa = self.k_Na_Ca * b * p
I_pCa = self.I_max_pCa * ((Cai) ** 2 / ((self.km_p_ca) ** 2 + (Cai) ** 2))
I_CaL = self.gCaL * Os * (V - self.ECaL)
dP_RyRdt = -0.04 * P_RyR - 0.1 * (I_CaL / self.I_CaL_max) * np.exp(
-((V - 5.0) ** 2) / 648
)
P_C1 = 1 - (P_C2 + P_O1 + P_O2)
dP_O1dt = self.k_a_positive * (Cass) ** self.n_ryr * P_C1
-self.k_a_negative * P_O1 - self.k_b_positive * (Cass) ** self.m_ryr * P_O1
+self.k_b_negative * P_O2 - self.k_c_positive * P_O1 + self.k_c_negative * P_C2
dP_O2dt = (
self.k_b_positive * (Cass) ** self.m_ryr * P_O1 - self.k_b_negative * P_O2
)
dP_C2dt = self.k_c_positive * P_O1 - self.k_c_negative * P_C2
dCassdt = Bss * (
-J_xfer * self.V_myo / self.V_ss - I_CaL * self.current2flux / self.V_ss
)
dOsdt = alpha * C4 - 4 * beta * Os + self.Kpcb * I1
-gamma * Os + 0.001 * (alpha * I2 - Kpcf * Os)
ECaN = (R * T) / (2 * F) * np.log(self.Cao / Cai)
I_Cab = self.Gcab * (V - ECaN)
"""
dCaidt = Bi * (J_leak + J_xfer - J_up - J_trpn
- (I_Cab - 2 * I_NaCa + I_pCa) * ((self.A_cap * self.C_mem)/(2 * self.V_myo * F)))
"""
# J_rel_caf = self.nu_1 * (self.Ca_JSR - Cass) #lisatud ette ennetavalt
JFCa = self.k_on * (self.F_tot - FCa) * Cai - self.k_off * FCa
dFCadt = JFCa
dCaidt = Bi * (
J_leak
+ J_xfer
- J_up
- J_trpn
- JFCa
- (I_Cab - 2 * I_NaCa + I_pCa)
* ((self.A_cap * self.C_mem) / (2 * self.V_myo * F))
)
# Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Ca_NSR, Ca_JSR, P_RyR, P_O1, P_O2, P_C2, Cai, FCa = states
return [
dCassdt,
dOsdt,
dC2dt,
dC3dt,
dC4dt,
dI1dt,
dI2dt,
dI3dt,
dLTRPNCadt,
dHTRPNCadt,
dCa_NSRdt,
dCa_JSRdt,
dP_RyRdt,
dP_O1dt,
dP_O2dt,
dP_C2dt,
dCaidt,
dFCadt,
]
def mem_potential(self, t: np.ndarray) -> np.ndarray:
t = np.asarray(t)
v1 = 0
t0 = 100
t1 = 350
n = (t // self.period) * self.period
return np.where((t0 + n <= t) & (t < t1 + n), v1, self.V_mem_rest)
def get_initial_values(self) -> list[float]:
Cai_0: float = 0.11712 # Cai
FCa_0: float = (self.k_on * self.F_tot * Cai_0) / (
self.k_on * Cai_0 + self.k_off
)
return [
Cai_0, # Cass
0.930308e-18, # Os,
0.124216e-3, # C2,
0.578679e-8, # C3,
0.119816e-12, # C4,
0.497923e-18, # I1,
0.345847e-13, # I2,
0.185106e-13, # I3,
11.2684, # LTRPNCa,
125.290, # HTRPNCa
1299.50, # Ca_NSR
1299.50, # Ca_JSR
0.0, # P_RyR
0.149102e-4, # P_O1
0.951726e-10, # P_02
0.167740e-3, # P_C2
Cai_0, # Cai
FCa_0, # FCa
]
def calculated_current(self, states: list[float], V: np.ndarray) -> np.ndarray:
Os = states[1]
I_CaL = self.gCaL * Os * (V - self.ECaL)
return I_CaL
def solve(self, initial_values: list[float], tspan, dt, times):
times = np.arange(*tspan, dt)
r = ode(self.ode_system)
r.set_integrator(
"lsoda", method="bdf", atol=1e-07, rtol=1e-07, max_step=0.1, nsteps=500
)
r.set_initial_value(initial_values, times[0])
states = np.zeros((len(initial_values), times.size))
states[:, 0] = initial_values
for i, t in enumerate(times[1:]):
if r.successful():
r.integrate(t)
states[:, i + 1] = r.y
else:
break
V = np.array([self.mem_potential(t) for t in times])
# plt.plot(times, V)
fig = plt.figure()
ax1 = fig.add_subplot(141)
ax2 = fig.add_subplot(142)
ax3 = fig.add_subplot(143)
ax4 = fig.add_subplot(144)
ax1.plot(times, states[0, :], label="Cass")
Cai = states[-2, :]
ax1.plot(times, (states[0, :] - Cai) / self.tau_xfer, label="Jxfer")
ax1.plot(
times,
-self.gCaL
* states[1, :]
* (V - self.ECaL)
* self.current2flux
/ self.V_ss
/ 1000,
label="JCaL mM/s",
)
ax1.legend(frameon=False)
ax1.set_xlabel("time [ms]")
ax1.set_ylabel(r"$\mu mol/l$")
ax2.plot(times, self.gCaL * states[1, :] * (V - self.ECaL), label="ICaL ")
ax2.legend(frameon=False)
ax2.set_xlabel("time [ms]")
ax2.set_ylabel("pA/pF")
ax3.plot(times, Cai, label="Cai ")
ax3.legend(frameon=False)
ax3.set_xlabel("time [ms]")
ax3.set_ylabel(r"$\mu mol/3l$")
ax4.plot(times, states[-1, :], label="Fca")
ax4.legend(frameon=False)
ax4.set_xlabel("time [ms]")
ax4.set_ylabel(r"$\mu mol/3l$")
return states, V
if __name__ == "__main__":
# Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Cai, FCa= states
model = Model()
Cai_0 = 0.11712 # Cai
FCa_0 = (model.k_on * model.F_tot * Cai_0) / (model.k_on * Cai_0 + model.k_off)
# print(FCa_0)
# Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Ca_NSR, Ca_JSR, P_RyR, P_O1, P_O2, P_C2, Cai, FCa = states
initial_values = [
Cai_0, # Cass
0.930308e-18, # Os,
0.124216e-3, # C2,
0.578679e-8, # C3,
0.119816e-12, # C4,
0.497923e-18, # I1,
0.345847e-13, # I2,
0.185106e-13, # I3,
11.2684, # LTRPNCa,
125.290, # HTRPNCa
1299.50, # Ca_NSR
1299.50, # Ca_JSR
0.0, # P_RyR
0.149102e-4, # P_O1
0.951726e-10, # P_02
0.167740e-3, # P_C2
Cai_0, # Cai
FCa_0, # FCa
]
model.solve(initial_values=initial_values, tspan=[0, 1000], dt=1.0, times=None)
plt.show()
plt.close('all')
Reference in New Issue View Git Blame Copy Permalink