From 0e504c01766b5e12ff32447d18a602bd549658f3 Mon Sep 17 00:00:00 2001 From: Ralf Magnus Palmiste Date: Thu, 11 Jul 2024 09:59:57 +0300 Subject: [PATCH] Upload files to "/" Contains all of the variables and collection of differential equations required to model calcium dynamics inside a mouse ventricular myocyte. --- Model.py | 407 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 407 insertions(+) create mode 100644 Model.py diff --git a/Model.py b/Model.py new file mode 100644 index 0000000..627e59c --- /dev/null +++ b/Model.py @@ -0,0 +1,407 @@ +# -*- coding: utf-8 -*- + +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import ode + +# %% Constants + +R = 8.314 # ideal gas constant (J/(mol*K)) +F = 96.485 # coulomb_per_millimole, Faraday constant +T = 298 # room temperature (K) +RT = R * T + + +# %% model +class Model: + def __init__(self): + + self.A_cap = 1.534e-4 # cm^2, Capacitive membrane area + self.C_mem = 1.0 # uF/cm2 Specific membrane capacitance + self.V_myo = 25.84e-6 # uL, Myoplasmic volume + self.current2flux = self.A_cap * self.C_mem / 2 / F # + self.period = 1000 # ms, pulse period + self.V_mem_rest = -80.0 # mV, resting membrane potential + self.Nai = 11000 # uM, Myoplasmic Na+ concentration + self.Nao = 150000 # uM, Extracellular Na+ concentration + self.eta = 0.35 # Controls voltage dependance of Na/Ca2+ exchange + self.km_Na = 87500 # uM, Na+ half-saturation constant for Na+/Ca2+ exchange + self.k_sat = ( + 0.1 # Na+/Ca2+ exchange saturation factor at very negative potentials + ) + self.km_Ca = (1380,) # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange + self.k_Na_Ca = 292.8 # pA/pF, Scaling factor of Na2+/Ca2+ exchange + self.I_max_pCa = 1.0 # pA/pF, Maximum Ca2+ pump current + # self.Cai = 0.1 # uM, Cytoplasmic Ca2+ concentration fixed at this point + self.LTRPN_tot = ( + 70.0 # uM, Total myoplasmic troponin low-affinity site concentration + ) + self.HTRPN_tot = ( + 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 = ( + 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 = ( + 0.0196 # ms, Ca2+ off rate const. for troponin low-affinity sites + ) + self.CMDN_tot = 50 # uM, Total myoplasmic calmodulin concentration + self.Km_CMDN = 0.238 # uM, Ca2 half-saturation constant for calmodulin + self.nu_1 = 4.5 # 1/ms, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR + self.nu_2 = 1.74e-5 # ms^(-1), Ca2+ leak rate const. from the NSR + self.nu_3 = 0.45 # uM/ms, SR Ca2+ -ATPase maximum pupmp rate + self.Km_up = 0.5 # uM, Half-saturation constant for SR Ca2+ -ATPase pump + self.km_p_ca = 0.5 # uM, Ca2+ half-saturation constant for Ca2+ pump current + # self.Ca_NSR = 1299.50 # uM,NSR Ca2+ concentration + self.tau_xfer = 8.0 # ms, Time constant for transfer from subspace to myoplasm + self.K_pc_max = ( + 0.23324 # 1/ms, Maximum time constant for Ca2+-induced inactivation + ) + self.K_pc_half = ( + 20.0 # uM, Half-saturation constant for Ca2+-induced inactivation + ) + self.Kpcb = 0.0005 # 1/ms, Voltage-insensitive rate constant for inactivation + self.km_Ca = 1380 # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange + self.Gcab = 0.000367 # mS/uF + self.Cao = 1130.0 # uM, Ca2+ outside the cell + self.gCaL = ( + 0.1729 # mS/uF, Specific maximum conductivity for L-type Ca2+ channel + ) + self.ECaL = 43.0 # mV, Reversal potential for L-type Ca2⫹ channel, kas arvutame voi jaab const? + self.V_ss = 1.485e-9 # uL, Dyadic aka subspace volume + # self.Ca_JSR = 1299.50 # uM, JSR Ca2+ concentration + self.F_tot = 25 # uM, total concentration of Fluo-4 + self.k_on = 0.1 # 1/uM * 1/ms, Fluo-4 reaction rate constant + self.k_off = 0.11 # 1/ms, Fluo-4 reaction rate constant 2 + self.V_NSR = 2.098e-6 # ul, Network SR volume + self.tau_tr = 20 # ms, Time const for transfer from NSR to JSR + self.CSQN_tot = 15000.0 # uM, total junctional SR calsequestrin concentration + self.Km_CSQN = 800.0 # uM, Ca2 half-saturation constant for calsequestrin + self.k_a_positive = 0.006075 # uM^(-4)/ms, RyR Pc1 - Po1 rate constant + self.k_a_negative = 0.07125 # 1/ms, RyR Po1 - Pc1 rate constant + self.k_b_positive = 0.00405 # uM^(-3)/ms, RyR Po1 - Po2 rate constant + self.k_b_negative = 0.965 # 1/ms, RyR Po2 - Po1 rate constant + self.k_c_positive = 0.009 # 1/ms, RyR Po1 - Pc2 rate constant + self.k_c_negative = 0.0008 # 1/ms, RyR Pc2 - Po1 rate constant + self.n_ryr = 4 # RyR Ca2+ cooperativity parameter Pc1 - Po1 + self.m_ryr = 3 # RyR Ca2+ cooperativity parameter Po1 - Po2 + self.I_CaL_max = 7.0 # pA/pF, normalization constant for L-type Ca2+ current + self.V_JSR = 0.12e-6 # ul, Junctional SR volume + + def ode_system(self, t, states): + + ( + 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 + (self.CMDN_tot * self.Km_CMDN)/(self.Km_CMDN+(Cai))**2)**(-1) + # Bss = 1 / ((self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2) + + Bi = 1 / ( + 1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cai) ** 2 + ) # A6 !!!!!!! korras V8 + + Bss = 1 / ( + 1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2 + ) # A7 !!!!!!!!! korras V8 + + 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) # !!!!!! korras v8 + - 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) # puudu CaNSR diff vorrand + + 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_a_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) / 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) + + 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): + v1 = 0 + t0 = 100 + t1 = 350 + n = (t // self.period) * self.period + return v1 if t0 + n <= t < t1 + n else self.V_mem_rest + + def solve(self, initial_values, tspan, dt): + 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.array([[0.0] * times.size] * len(initial_values)) + 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$") + + plt.show() + + +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)