Calcium_Model/Model.py

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 point

        self.LTRPN_tot: float = (
            70.0  # uM, Total myoplasmic troponin low-affinity site concentr.
        )
        self.HTRPN_tot: float = (
            140.0  # uM, Total myoplasmic troponin high-affinity site concentr.
        )

        """
        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-affin 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-affin 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 concentr.
        self.Km_CMDN: float = 0.238  # uM, Ca2 half-sat const for calmodulin

        self.nu_1: float = (
            4.5  # 1/ms, Max RyR chan Ca2+ perm. Ca2+ leak rate const from NSR
        )
        self.nu_2: float = 1.74e-5  # ms^(-1), Ca2+ leak rate const. from NSR
        self.nu_3: float = 0.45  # uM/ms, SR Ca2+ -ATPase maximum pupmp rate

        self.Km_up: float = 0.5  # uM, Half-sat const 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 inact.
        )
        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 max conduct. 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
        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
        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 param. 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))
        )

        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])
        return states, V

    def plot_results(self, times, states, V):
        # 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()

    initial_values = model.get_initial_values()
    times = np.arange(0, 1000, 1.0)

    states, V = model.solve(
        initial_values=initial_values, tspan=[0, 1000], dt=1.0, times=times
    )

    model.plot_results(times, states, V)