In [15]:

    

import sympy as sp
sp.init_printing(use_latex='png')


    

In [8]:

    

# Subspace calcium states for individual cicr states
CASS_INDEX = [2 * int(i == 2) + int(j == 2) for i in range(4) for j in range(10)]
CASS_INDEX_INV = [[k for k, val in enumerate(CASS_INDEX) if val == i] for i in range(4)]


    

In [9]:

    

# Sympy walkthrough


    

In [20]:

    

# one symbol
p = sp.symbols('p_0:4')
p


    

    
Out[20]:

In [11]:

    

# two symbols
r_ryr, r_xfer = sp.symbols('r_RyR, r_xfer')  # Separated by Comma or space 
r_ryr ** r_xfer  # Supports common math operations


    

    
Out[11]:





$$r_{RyR}^{r_{xfer}}$$

In [17]:

    

# functions
t = sp.Symbol('t')
f = sp.Function('f')(t)  # f = f(t)
g = sp.symbols('g', cls=sp.Function)(f)  # g = g(f(t))
sp.diff(g, t) # chain-rule


    

    
Out[17]:





d        ⎛ d        ⎞│       
──(f(t))⋅⎜───(g(ξ₁))⎟│       
dt       ⎝dξ₁       ⎠│ξ₁=f(t)

In [13]:

    

sp.Eq(sp.symbols('x'), 1)


    

    
Out[13]:





$$x = 1$$

In [14]:

    

ca_ss = [sp.Symbol("[Ca^{2+}]_{ss[%d]}" % (i), nonnegative=True) for i in range(4)]
ca_ss[0]**2


    

    
Out[14]:





$$[Ca^{2+}]_{ss[0]}^{2}$$

In [48]:

    

# Symbols names
simple_names = ["f", "g", "a", "b", "gamma_0", "omega", "alpha", "beta",
                "alpha_a", "beta_b", "y_inf", "tau_y", "k_b", "k_f", "omega_b", 
                "F", "R", "T", "t", "P_LCC", "r_RyR", "r_xfer"]

# symbols of RyR transition constants
RYR_NUM_PAIR = [(a, a + 1) for a in range(5)] + [(1, 4)] 
simple_names.append("K0_[%d,%d]" % (i, j) for (i, j) in RYR_NUM_PAIR)
simple_names.append("K0_[%d,%d]" % (j, i) for (i, j) in RYR_NUM_PAIR)

# Create symbols from names
syms = {n: sp.symbols(n, nonnegative=True) for n in simple_names}

# Indivisual names and symbols
vm = syms["V_m"] = sp.Symbol("V_m", real=True)
ca_ss = sp.Matrix([sp.Symbol("[Ca^{2+}]_{ss[%d]}" % (i), nonnegative=True) for i in range(4)])
ca_ss_sq = ca_ss.multiply_elementwise(ca_ss)



# Dependent relationships between symbols
deps = [(syms["alpha"], 0.78113 * sp.exp(0.08 * (vm - 8))),
        (syms["beta"],  0.34319 * sp.exp(-0.086676 * (vm - 8))),
        (syms["alpha_a"], syms["a"] * syms["alpha"]),
        (syms["beta_b"], syms["beta"] / syms["b"]),
        (syms["y_inf"], 0.25 + 0.75 / (1 + sp.exp((vm + 25) / 5))),
        (syms["tau_y"], 17 + 500 / ((1 + sp.exp((vm + 28) / 10)) * (1 + sp.exp(-(vm + 42) / 14)))),
        (syms["k_b"], syms["y_inf"] / syms["tau_y"]),
        (syms["k_f"], (1 - syms["y_inf"]) / syms["tau_y"]),
        (syms["omega_b"], syms["omega"] / syms["b"])]

# Equations for printing
equations = {str(lhs): sp.Eq(lhs, rhs, evaluate=False) for (lhs, rhs) in deps}


    

In [49]:

    

ca_ss_sq


    

    
Out[49]:





$$\left[\begin{matrix}[Ca^{2+}]_{ss[0]}^{2}\\[Ca^{2+}]_{ss[1]}^{2}\\[Ca^{2+}]_{ss[2]}^{2}\\[Ca^{2+}]_{ss[3]}^{2}\end{matrix}\right]$$

In [70]:

    

# Symbols names
simple_names = ["f", "g", "a", "b", "gamma_0", "omega", "alpha", "beta",
                "alpha_a", "beta_b", "y_inf", "tau_y", "k_b", "k_f", "omega_b", 
                "F", "R", "T", "t", "P_LCC", "r_RyR", "r_xfer", "V_SS"]

syms = {n: sp.symbols(n, nonnegative=True) for n in simple_names}
vm = syms["V_m"] = sp.Symbol("V_m", real=True)
ca_ss = syms["ca_ss"] = sp.Matrix([sp.Symbol("[Ca^{2+}]_{ss[%d]}" % (i), nonnegative=True) for i in range(4)])

# Dependent relationships of symbols
deps = [(syms["alpha"], 0.78113 * sp.exp(0.08 * (vm - 8))),
        (syms["beta"],  0.34319 * sp.exp(-0.086676 * (vm - 8))),
        (syms["alpha_a"], syms["a"] * syms["alpha"]),
        (syms["beta_b"], syms["beta"] / syms["b"]),
        (syms["y_inf"], 0.25 + 0.75 / (1 + sp.exp((vm + 25) / 5))),
        (syms["tau_y"], 17 + 500 / ((1 + sp.exp((vm + 28) / 10)) * (1 + sp.exp(-(vm + 42) / 14)))),
        (syms["k_b"], syms["y_inf"] / syms["tau_y"]),
        (syms["k_f"], (1 - syms["y_inf"]) / syms["tau_y"]),
        (syms["omega_b"], syms["omega"] / syms["b"])]

equations = {str(lhs): sp.Eq(lhs, rhs, evaluate=False) for (lhs, rhs) in deps}


    

In [78]:

    

print(sp.latex(equations["k_f"]))


    

    




k_{f} = \frac{1}{\tau_{y}} \left(- y_{inf} + 1\right)

In [54]:

    

def add_rates(jac, i, j, kf, kb):
    """
    Adds 1st order reaction constants to
    the jacobian matrix of Markov-type ODE systems
    """
    jac[i, i] -= kf
    jac[j, j] -= kb
    jac[i, j] += kb
    jac[j, i] += kf
    return jac


# Jacobian matrix of the ODE system
jacobian = sp.zeros(40)

# Fill in the LCC part
for r_idx in range(0, 40, 10):
    g_cass = syms["gamma_0"] * syms["ca_ss"][CASS_INDEX[r_idx]]
    g_a_cass = syms["a"] * g_cass
    for i in range(r_idx, r_idx + 5):
        jacobian = add_rates(jacobian, i, i + 5, syms["k_f"], syms["k_b"])
    for i in (0 + r_idx, 5 + r_idx):
        jacobian = add_rates(jacobian, i, i + 1, syms["alpha"], syms["beta"])
        jacobian = add_rates(jacobian, i, i + 3, g_cass, syms["omega"])
    for i in (1 + r_idx, 6 + r_idx):
        jacobian = add_rates(jacobian, i, i + 1, syms["f"], syms["g"])
        jacobian = add_rates(jacobian, i, i + 3, g_a_cass, syms["omega_b"])
    for i in (3 + r_idx, 8 + r_idx):
        jacobian = add_rates(jacobian, i, i + 1, syms["alpha_a"], syms["beta_b"])


# RyR transition constants
for (i, j) in [(a, a + 1) for a in range(5)] + [(1, 4)]:
    syms[f"K0_{i}{j}"] = sp.Symbol(f"K_0[{i},{j}]", nonnegative=True)
    syms[f"K0_{j}{i}"] = sp.Symbol(f"K_0[{j},{i}]", nonnegative=True)


def k_01(ca, syms):
    return syms["K0_01"] * (ca**2)


def k_10(ca, syms):
    return syms["K0_10"]


def k_12(ca, syms):
    return syms["K0_12"] * (ca**2)


def k_21(ca, syms):
    return syms["K0_21"] * syms["K0_32"] / (syms["K0_23"] * (ca**2) + syms["K0_32"])


def k_13(ca, syms):
    return syms["K0_14"] * (ca**2)


def k_31(ca, syms):
    return syms["K0_41"] * syms["K0_54"] / (syms["K0_45"] * (ca**2) + syms["K0_54"])


def k_23(ca, syms):
    return syms["K0_34"] * syms["K0_23"] * (ca**2) / (syms["K0_32"] + syms["K0_23"] * (ca**2))


def k_32(ca, syms):
    return syms["K0_54"] * syms["K0_43"] * (ca**2) / (syms["K0_54"] + syms["K0_45"] * (ca**2))


for l_idx in range(0, 10):
    ca = [syms["ca_ss"][CASS_INDEX[l_idx + i*10]] for i in range(4)]
    jacobian = add_rates(jacobian, 0 + l_idx, 10 + l_idx, k_01(ca[0], syms), k_10(ca[1], syms))
    jacobian = add_rates(jacobian, 10 + l_idx, 20 + l_idx, k_12(ca[1], syms), k_21(ca[2], syms))
    jacobian = add_rates(jacobian, 20 + l_idx, 30 + l_idx, k_23(ca[2], syms), k_32(ca[3], syms))
    jacobian = add_rates(jacobian, 10 + l_idx, 30 + l_idx, k_13(ca[1], syms), k_31(ca[3], syms))


    

In [79]:

    

x_states = sp.Matrix(sp.symbols('x0:40'), nonnegative=True)
dxdt = jacobian * x_states
time_sym = sp.Symbol('t', nonnegative=True)
for i, x in enumerate(x_states):
    x_states[i] = x(time_sym)


    

In [81]:

    

print(sp.latex(dxdt))


    

    




\left[\begin{matrix}K_{0[1,0]} x_{10} + \beta x_{1} + k_{b} x_{5} + \omega x_{3} + x_{0} \left(- ^{2} K_{0[0,1]} - \gamma_{0} - \alpha - k_{f}\right)\\K_{0[1,0]} x_{11} + \alpha x_{0} + g x_{2} + k_{b} x_{6} + \omega_{b} x_{4} + x_{1} \left(- ^{2} K_{0[0,1]} - a \gamma_{0} - \beta - f - k_{f}\right)\\K_{0[1,0]} x_{12} + f x_{1} + k_{b} x_{7} + x_{2} \left(- ^{2} K_{0[0,1]} - g - k_{f}\right)\\\gamma_{0} x_{0} + K_{0[1,0]} x_{13} + \beta_{b} x_{4} + k_{b} x_{8} + x_{3} \left(- ^{2} K_{0[0,1]} - \alpha_{a} - k_{f} - \omega\right)\\a \gamma_{0} x_{1} + K_{0[1,0]} x_{14} + \alpha_{a} x_{3} + k_{b} x_{9} + x_{4} \left(- ^{2} K_{0[0,1]} - \beta_{b} - k_{f} - \omega_{b}\right)\\K_{0[1,0]} x_{15} + \beta x_{6} + k_{f} x_{0} + \omega x_{8} + x_{5} \left(- ^{2} K_{0[0,1]} - \gamma_{0} - \alpha - k_{b}\right)\\K_{0[1,0]} x_{16} + \alpha x_{5} + g x_{7} + k_{f} x_{1} + \omega_{b} x_{9} + x_{6} \left(- ^{2} K_{0[0,1]} - a \gamma_{0} - \beta - f - k_{b}\right)\\K_{0[1,0]} x_{17} + f x_{6} + k_{f} x_{2} + x_{7} \left(- ^{2} K_{0[0,1]} - g - k_{b}\right)\\\gamma_{0} x_{5} + K_{0[1,0]} x_{18} + \beta_{b} x_{9} + k_{f} x_{3} + x_{8} \left(- ^{2} K_{0[0,1]} - \alpha_{a} - k_{b} - \omega\right)\\a \gamma_{0} x_{6} + K_{0[1,0]} x_{19} + \alpha_{a} x_{8} + k_{f} x_{4} + x_{9} \left(- ^{2} K_{0[0,1]} - \beta_{b} - k_{b} - \omega_{b}\right)\\^{2} K_{0[0,1]} x_{0} + \frac{K_{0[2,1]} K_{0[3,2]} x_{20}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{30}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \beta x_{11} + k_{b} x_{15} + \omega x_{13} + x_{10} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - \gamma_{0} - K_{0[1,0]} - \alpha - k_{f}\right)\\^{2} K_{0[0,1]} x_{1} + \frac{K_{0[2,1]} K_{0[3,2]} x_{21}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{31}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \alpha x_{10} + g x_{12} + k_{b} x_{16} + \omega_{b} x_{14} + x_{11} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - a \gamma_{0} - K_{0[1,0]} - \beta - f - k_{f}\right)\\^{2} K_{0[0,1]} x_{2} + \frac{K_{0[2,1]} K_{0[3,2]} x_{22}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{32}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + f x_{11} + k_{b} x_{17} + x_{12} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - K_{0[1,0]} - g - k_{f}\right)\\^{2} K_{0[0,1]} x_{3} + \gamma_{0} x_{10} + \frac{K_{0[2,1]} K_{0[3,2]} x_{23}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{33}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \beta_{b} x_{14} + k_{b} x_{18} + x_{13} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - K_{0[1,0]} - \alpha_{a} - k_{f} - \omega\right)\\^{2} K_{0[0,1]} x_{4} + a \gamma_{0} x_{11} + \frac{K_{0[2,1]} K_{0[3,2]} x_{24}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{34}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \alpha_{a} x_{13} + k_{b} x_{19} + x_{14} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - K_{0[1,0]} - \beta_{b} - k_{f} - \omega_{b}\right)\\^{2} K_{0[0,1]} x_{5} + \frac{K_{0[2,1]} K_{0[3,2]} x_{25}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{35}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \beta x_{16} + k_{f} x_{10} + \omega x_{18} + x_{15} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - \gamma_{0} - K_{0[1,0]} - \alpha - k_{b}\right)\\^{2} K_{0[0,1]} x_{6} + \frac{K_{0[2,1]} K_{0[3,2]} x_{26}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{36}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \alpha x_{15} + g x_{17} + k_{f} x_{11} + \omega_{b} x_{19} + x_{16} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - a \gamma_{0} - K_{0[1,0]} - \beta - f - k_{b}\right)\\^{2} K_{0[0,1]} x_{7} + \frac{K_{0[2,1]} K_{0[3,2]} x_{27}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{37}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + f x_{16} + k_{f} x_{12} + x_{17} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - K_{0[1,0]} - g - k_{b}\right)\\^{2} K_{0[0,1]} x_{8} + \gamma_{0} x_{15} + \frac{K_{0[2,1]} K_{0[3,2]} x_{28}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{38}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \beta_{b} x_{19} + k_{f} x_{13} + x_{18} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - K_{0[1,0]} - \alpha_{a} - k_{b} - \omega\right)\\^{2} K_{0[0,1]} x_{9} + a \gamma_{0} x_{16} + \frac{K_{0[2,1]} K_{0[3,2]} x_{29}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \frac{K_{0[4,1]} K_{0[5,4]} x_{39}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \alpha_{a} x_{18} + k_{f} x_{14} + x_{19} \left(- ^{2} K_{0[1,2]} - ^{2} K_{0[1,4]} - K_{0[1,0]} - \beta_{b} - k_{b} - \omega_{b}\right)\\^{2} K_{0[1,2]} x_{10} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{30}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \beta x_{21} + k_{b} x_{25} + \omega x_{23} + x_{20} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \gamma_{0} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \alpha - k_{f}\right)\\^{2} K_{0[1,2]} x_{11} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{31}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \alpha x_{20} + g x_{22} + k_{b} x_{26} + \omega_{b} x_{24} + x_{21} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - a \gamma_{0} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \beta - f - k_{f}\right)\\^{2} K_{0[1,2]} x_{12} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{32}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + f x_{21} + k_{b} x_{27} + x_{22} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - g - k_{f}\right)\\^{2} K_{0[1,2]} x_{13} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{33}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \gamma_{0} x_{20} + \beta_{b} x_{24} + k_{b} x_{28} + x_{23} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \alpha_{a} - k_{f} - \omega\right)\\^{2} K_{0[1,2]} x_{14} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{34}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + a \gamma_{0} x_{21} + \alpha_{a} x_{23} + k_{b} x_{29} + x_{24} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \beta_{b} - k_{f} - \omega_{b}\right)\\^{2} K_{0[1,2]} x_{15} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{35}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \beta x_{26} + k_{f} x_{20} + \omega x_{28} + x_{25} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \gamma_{0} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \alpha - k_{b}\right)\\^{2} K_{0[1,2]} x_{16} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{36}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \alpha x_{25} + g x_{27} + k_{f} x_{21} + \omega_{b} x_{29} + x_{26} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - a \gamma_{0} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \beta - f - k_{b}\right)\\^{2} K_{0[1,2]} x_{17} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{37}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + f x_{26} + k_{f} x_{22} + x_{27} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - g - k_{b}\right)\\^{2} K_{0[1,2]} x_{18} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{38}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + \gamma_{0} x_{25} + \beta_{b} x_{29} + k_{f} x_{23} + x_{28} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \alpha_{a} - k_{b} - \omega\right)\\^{2} K_{0[1,2]} x_{19} + \frac{^{2} K_{0[4,3]} K_{0[5,4]} x_{39}}{^{2} K_{0[4,5]} + K_{0[5,4]}} + a \gamma_{0} x_{26} + \alpha_{a} x_{28} + k_{f} x_{24} + x_{29} \left(- \frac{^{2} K_{0[2,3]} K_{0[3,4]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \frac{K_{0[2,1]} K_{0[3,2]}}{^{2} K_{0[2,3]} + K_{0[3,2]}} - \beta_{b} - k_{b} - \omega_{b}\right)\\^{2} K_{0[1,4]} x_{10} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{20}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \beta x_{31} + k_{b} x_{35} + \omega x_{33} + x_{30} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \gamma_{0} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \alpha - k_{f}\right)\\^{2} K_{0[1,4]} x_{11} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{21}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \alpha x_{30} + g x_{32} + k_{b} x_{36} + \omega_{b} x_{34} + x_{31} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - a \gamma_{0} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \beta - f - k_{f}\right)\\^{2} K_{0[1,4]} x_{12} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{22}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + f x_{31} + k_{b} x_{37} + x_{32} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - g - k_{f}\right)\\^{2} K_{0[1,4]} x_{13} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{23}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \gamma_{0} x_{30} + \beta_{b} x_{34} + k_{b} x_{38} + x_{33} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \alpha_{a} - k_{f} - \omega\right)\\^{2} K_{0[1,4]} x_{14} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{24}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + a \gamma_{0} x_{31} + \alpha_{a} x_{33} + k_{b} x_{39} + x_{34} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \beta_{b} - k_{f} - \omega_{b}\right)\\^{2} K_{0[1,4]} x_{15} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{25}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \beta x_{36} + k_{f} x_{30} + \omega x_{38} + x_{35} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \gamma_{0} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \alpha - k_{b}\right)\\^{2} K_{0[1,4]} x_{16} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{26}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \alpha x_{35} + g x_{37} + k_{f} x_{31} + \omega_{b} x_{39} + x_{36} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - a \gamma_{0} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \beta - f - k_{b}\right)\\^{2} K_{0[1,4]} x_{17} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{27}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + f x_{36} + k_{f} x_{32} + x_{37} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - g - k_{b}\right)\\^{2} K_{0[1,4]} x_{18} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{28}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + \gamma_{0} x_{35} + \beta_{b} x_{39} + k_{f} x_{33} + x_{38} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \alpha_{a} - k_{b} - \omega\right)\\^{2} K_{0[1,4]} x_{19} + \frac{^{2} K_{0[2,3]} K_{0[3,4]} x_{29}}{^{2} K_{0[2,3]} + K_{0[3,2]}} + a \gamma_{0} x_{36} + \alpha_{a} x_{38} + k_{f} x_{34} + x_{39} \left(- \frac{^{2} K_{0[4,3]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \frac{K_{0[4,1]} K_{0[5,4]}}{^{2} K_{0[4,5]} + K_{0[5,4]}} - \beta_{b} - k_{b} - \omega_{b}\right)\end{matrix}\right]

In [91]:

    

def make_eq(lhs_name, rhs_expr, evaluate=False, **kwargs):
    """
    Make an equation from name
    """
    return sp.Eq(sp.Symbol(lhs_name, **kwargs), rhs_expr, evaluate=evaluate)

vm = sp.Symbol("V_m", real=True)
ca_out = sp.Symbol("[Ca^{2+}]_{o}", nonnegative=True)
ca_in = sp.Symbol("[Ca^{2+}]_{in}", nonnegative=True)
ca_jsr = sp.Symbol("[Ca^{2+}]_{JSR}", nonnegative=True)
ca_nsr = sp.Symbol("[Ca^{2+}]_{NSR}", nonnegative=True)
R_RYR = sp.Symbol("r_{RyR}", nonnegative=True)
R_XFER = sp.Symbol("r_{xfer}", nonnegative=True)
p_cicr = sp.symbols("p_{0:4}", nonnegative=True)
P_LCC = sp.Symbol("P_{LCC}", nonnegative=True)
V_SS = sp.Symbol("V_{SS}", nonnegative=True)
F = sp.Symbol("F", nonnegative=True)
R = sp.Symbol("R", nonnegative=True)
T = sp.Symbol("T", nonnegative=True)

v2FRT = 2 * vm * F / R / T
e2vfrt = sp.E ** v2FRT
ca_out_weight = P_LCC / V_SS * v2FRT / (e2vfrt - 1)
ca_out_value = 0.341 * ca_out * ca_out_weight
ca_out_weight *= e2vfrt

cain_weight = R_XFER
ca_in_value = ca_in * cain_weight
ca_jsr_weight = R_RYR
ca_jsr_value = ca_jsr * ca_jsr_weight

ca_ss = [0] * 4
ca_ss[0] = make_eq("[Ca^{2+}]_{ss[0]}", ca_in, nonnegative=True)
ca_ss[1] = make_eq("[Ca^{2+}]_{ss[1]}", (ca_out_value + ca_in_value) / (ca_out_weight + cain_weight), nonnegative=True)
ca_ss[2] = make_eq("[Ca^{2+}]_{ss[2]}", (ca_jsr_value + ca_in_value) / (ca_jsr_weight + cain_weight), nonnegative=True)
ca_ss[3] = make_eq("[Ca^{2+}]_{ss[3]}",
                (ca_out_value + ca_in_value + ca_jsr_value) / (ca_out_weight + cain_weight + ca_jsr_weight), nonnegative=True)


    

In [96]:

    

ca_ss[3]


    

    
Out[96]:

In [ ]: