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

Applied good practice tips; Fixed typos in expressions in model.py

Browse Source
This commit is contained in:
ralf 2024-08-08 11:01:19 +03:00
parent 89c492dd6c
commit bb67e34f73
3 changed files with 80 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

139
Model.py
View File

@ -11,41 +11,45 @@ RT = R * T
class Model: class Model:
def __init__(self): def __init__(self):
self.A_cap = 1.534e-4 # cm^2, Capacitive membrane area self.A_cap: float = 1.534e-4 # cm^2, Capacitive membrane area
self.C_mem = 1.0 # uF/cm2 Specific membrane capacitance self.C_mem: float = 1.0 # uF/cm2 Specific membrane capacitance
self.V_myo = 25.84e-6 # uL, Myoplasmic volume self.V_myo: float = 25.84e-6 # uL, Myoplasmic volume
self.current2flux = self.A_cap * self.C_mem / 2 / F # self.current2flux: float = self.A_cap * self.C_mem / 2 / F #
self.period = 1000 # ms, pulse period self.period: float = 1000 # ms, pulse period
self.V_mem_rest = -80.0 # mV, resting membrane potential self.V_mem_rest: float = -80.0 # mV, resting membrane potential
self.Nai = 11000 # uM, Myoplasmic Na+ concentration self.Nai: float = 11000 # uM, Myoplasmic Na+ concentration
self.Nao = 150000 # uM, Extracellular Na+ concentration self.Nao: float = 150000 # uM, Extracellular Na+ concentration
self.eta = 0.35 # Controls voltage dependance of Na/Ca2+ exchange self.eta: float = 0.35 # Controls voltage dependance of Na/Ca2+ exchange
self.km_Na = 87500 # uM, Na+ half-saturation constant for Na+/Ca2+ exchange self.km_Na: float = (
self.k_sat = ( 87500 # uM, Na+ half-saturation constant for Na+/Ca2+ exchange
)
self.k_sat: float = (
0.1 # Na+/Ca2+ exchange saturation factor at very negative potentials 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.km_Ca: float = (
self.k_Na_Ca = 292.8 # pA/pF, Scaling factor of Na2+/Ca2+ exchange 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 = 1.0 # pA/pF, Maximum Ca2+ pump current 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.Cai = 0.1 # uM, Cytoplasmic Ca2+ concentration fixed at this point
self.LTRPN_tot = ( self.LTRPN_tot: float = (
70.0 # uM, Total myoplasmic troponin low-affinity site concentration 70.0 # uM, Total myoplasmic troponin low-affinity site concentration
) )
self.HTRPN_tot = ( self.HTRPN_tot: float = (
140.0 # uM, Total myoplasmic troponin high-affinity site concentration 140.0 # uM, Total myoplasmic troponin high-affinity site concentration
) )
# self.LTRPNCa = 11.2684 # uM, Concentration Ca2+ bound low-affinity troponin-binding sites # 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.HTRPNCa = 125.290 # uM, Concentration Ca2+ bound high-affinity troponin-binding sites
self.k_htrpn_positive = ( self.k_htrpn_positive: float = (
0.00237 # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affinity sites 0.00237 # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affinity sites
) )
self.k_htrpn_negative = ( self.k_htrpn_negative = (
@ -55,65 +59,78 @@ class Model:
self.k_ltrpn_positive = ( self.k_ltrpn_positive = (
0.0327 # uM^(-1)/ms, Ca2+ on rate const. for troponin low-affinity sites 0.0327 # uM^(-1)/ms, Ca2+ on rate const. for troponin low-affinity sites
) )
self.k_ltrpn_negative = ( self.k_ltrpn_negative: float = (
0.0196 # ms, Ca2+ off rate const. for troponin low-affinity sites 0.0196 # ms, Ca2+ off rate const. for troponin low-affinity sites
) )
self.CMDN_tot = 50 # uM, Total myoplasmic calmodulin concentration self.CMDN_tot: float = 50 # uM, Total myoplasmic calmodulin concentration
self.Km_CMDN = 0.238 # uM, Ca2 half-saturation constant for calmodulin self.Km_CMDN: float = 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_1: float = (
self.nu_2 = 1.74e-5 # ms^(-1), Ca2+ leak rate const. from the NSR 4.5 # 1/ms, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR
self.nu_3 = 0.45 # uM/ms, SR Ca2+ -ATPase maximum pupmp rate )
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 = 0.5 # uM, Half-saturation constant for SR Ca2+ -ATPase pump self.Km_up: float = 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.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.Ca_NSR = 1299.50 # uM,NSR Ca2+ concentration
self.tau_xfer = 8.0 # ms, Time constant for transfer from subspace to myoplasm self.tau_xfer: float = (
8.0 # ms, Time constant for transfer from subspace to myoplasm
)
self.K_pc_max = ( self.K_pc_max: float = (
0.23324 # 1/ms, Maximum time constant for Ca2+-induced inactivation 0.23324 # 1/ms, Maximum time constant for Ca2+-induced inactivation
) )
self.K_pc_half = ( self.K_pc_half: float = (
20.0 # uM, Half-saturation constant for Ca2+-induced inactivation 20.0 # uM, Half-saturation constant for Ca2+-induced inactivation
) )
self.Kpcb = 0.0005 # 1/ms, Voltage-insensitive rate constant for inactivation self.Kpcb: float = (
self.km_Ca = 1380 # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange 0.0005 # 1/ms, Voltage-insensitive rate constant for inactivation
)
self.Gcab = 0.000367 # mS/uF self.Gcab: float = 0.000367 # mS/uF
self.Cao = 1130.0 # uM, Ca2+ outside the cell self.Cao: float = 1130.0 # uM, Ca2+ outside the cell
self.gCaL = ( self.gCaL: float = (
0.1729 # mS/uF, Specific maximum conductivity for L-type Ca2+ channel 0.1729 # mS/uF, Specific maximum conductivity for L-type Ca2+ channel
) )
self.ECaL = 43.0 # mV, Reversal potential for L-type Ca2 self.ECaL: float = 43.0 # mV, Reversal potential for L-type Ca2
self.V_ss = 1.485e-9 # uL, Dyadic aka subspace volume self.V_ss: float = 1.485e-9 # uL, Dyadic aka subspace volume
# self.Ca_JSR = 1299.50 # uM, JSR Ca2+ concentration # self.Ca_JSR = 1299.50 # uM, JSR Ca2+ concentration
self.F_tot = 25 # uM, total concentration of Fluo-4 self.F_tot: float = 25 # uM, total concentration of Fluo-4
self.k_on = 0.1 # 1/uM * 1/ms, Fluo-4 reaction rate constant self.k_on: float = 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.k_off: float = 0.11 # 1/ms, Fluo-4 reaction rate constant 2
self.V_NSR = 2.098e-6 # ul, Network SR volume self.V_NSR: float = 2.098e-6 # ul, Network SR volume
self.tau_tr = 20 # ms, Time const for transfer from NSR to JSR self.tau_tr: float = 20 # ms, Time const for transfer from NSR to JSR
self.CSQN_tot = 15000.0 # uM, total junctional SR calsequestrin concentration self.CSQN_tot: float = (
self.Km_CSQN = 800.0 # uM, Ca2 half-saturation constant for calsequestrin 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 = 0.006075 # uM^(-4)/ms, RyR Pc1 - Po1 rate constant self.k_a_positive: float = 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_a_negative: float = 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_positive: float = 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_b_negative: float = 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_positive: float = 0.009 # 1/ms, RyR Po1 - Pc2 rate constant
self.k_c_negative = 0.0008 # 1/ms, RyR Pc2 - Po1 rate constant self.k_c_negative: float = 0.0008 # 1/ms, RyR Pc2 - Po1 rate constant
self.n_ryr = 4 # RyR Ca2+ cooperativity parameter Pc1 - Po1 self.n_ryr: int = 4 # RyR Ca2+ cooperativity parameter Pc1 - Po1
self.m_ryr = 3 # RyR Ca2+ cooperativity parameter Po1 - Po2 self.m_ryr: int = 3 # RyR Ca2+ cooperativity parameter Po1 - Po2
self.I_CaL_max = 7.0 # pA/pF, normalization constant for L-type Ca2+ current self.I_CaL_max: float = (
self.V_JSR = 0.12e-6 # ul, Junctional SR volume 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): def ode_system(self, t, states):
@ -194,9 +211,6 @@ class Model:
- 4 * beta * self.Kpcb * I3 - 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) 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) Bss = 1 / (1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2)
@ -253,7 +267,7 @@ class Model:
P_C1 = 1 - (P_C2 + P_O1 + P_O2) P_C1 = 1 - (P_C2 + P_O1 + P_O2)
dP_O1dt = self.k_a_positive * (Cass) ** self.n_ryr * P_C1 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_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 +self.k_b_negative * P_O2 - self.k_c_positive * P_O1 + self.k_c_negative * P_C2
dP_O2dt = ( dP_O2dt = (
@ -269,7 +283,7 @@ class Model:
dOsdt = alpha * C4 - 4 * beta * Os + self.Kpcb * I1 dOsdt = alpha * C4 - 4 * beta * Os + self.Kpcb * I1
-gamma * Os + 0.001 * (alpha * I2 - Kpcf * Os) -gamma * Os + 0.001 * (alpha * I2 - Kpcf * Os)
ECaN = (R * T) / F * np.log(self.Cao / Cai) ECaN = (R * T) / (2 * F) * np.log(self.Cao / Cai)
I_Cab = self.Gcab * (V - ECaN) I_Cab = self.Gcab * (V - ECaN)
@ -326,8 +340,10 @@ class Model:
return np.where((t0 + n <= t) & (t < t1 + n), v1, self.V_mem_rest) return np.where((t0 + n <= t) & (t < t1 + n), v1, self.V_mem_rest)
def get_initial_values(self): def get_initial_values(self):
Cai_0 = 0.11712 # Cai Cai_0: float = 0.11712 # Cai
FCa_0 = (self.k_on * self.F_tot * Cai_0) / (self.k_on * Cai_0 + self.k_off) FCa_0: float = (self.k_on * self.F_tot * Cai_0) / (
self.k_on * Cai_0 + self.k_off
)
return [ return [
Cai_0, # Cass Cai_0, # Cass
0.930308e-18, # Os, 0.930308e-18, # Os,
@ -362,9 +378,6 @@ class Model:
) )
r.set_initial_value(initial_values, times[0]) r.set_initial_value(initial_values, times[0])
#states = np.array([[0.0] * times.size] * len(initial_values))
#states[:, 0] = initial_values
states = np.zeros((len(initial_values), times.size)) states = np.zeros((len(initial_values), times.size))
states[:, 0] = initial_values states[:, 0] = initial_values
@ -427,7 +440,7 @@ if __name__ == "__main__":
Cai_0 = 0.11712 # Cai Cai_0 = 0.11712 # Cai
FCa_0 = (model.k_on * model.F_tot * Cai_0) / (model.k_on * Cai_0 + model.k_off) FCa_0 = (model.k_on * model.F_tot * Cai_0) / (model.k_on * Cai_0 + model.k_off)
print(FCa_0) # 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 # 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 = [ initial_values = [

BIN
__pycache__/Model.cpython-311.pyc
View File

Binary file not shown.

8
fitter.py
View File

@ -23,13 +23,13 @@ class Fitter:
self.fit_results = {} self.fit_results = {}
self.tspan = [0, 1000] self.tspan = [0, 1000]
self.dt = 1 # 1.0 self.dt: float = 1 # 1.0
self.time_points = np.arange(*self.tspan, self.dt) self.time_points = np.arange(*self.tspan, self.dt)
model = Model() model = Model()
self.initial_values = model.get_initial_values() self.initial_values = model.get_initial_values()
self.iteration = 0 # least squares iteration counter self.iteration: int = 0 # least squares iteration counter
t0, t1 = self.current_fit_range = current_fit_range t0, t1 = self.current_fit_range = current_fit_range
self.current_time_indecies = (t0 <= self.time_points) & (self.time_points <= t1) self.current_time_indecies = (t0 <= self.time_points) & (self.time_points <= t1)
@ -104,8 +104,8 @@ class Fitter:
m = self.model() m = self.model()
K_pc_half = m.K_pc_half K_pc_half = m.K_pc_half
tau_xfer = m.tau_xfer tau_xfer = m.tau_xfer
tau_RC = 1.5 tau_RC: float = 1.5
offset = 0 offset: float = 0
d = self.data d = self.data
init_parameters = np.array( init_parameters = np.array(