Compare commits
9 Commits
a2bb97fe22
...
main
Author | SHA1 | Date | |
---|---|---|---|
4651e50dd1 | |||
224a704e84 | |||
4fa1acbb90 | |||
e11cc744d9 | |||
4173284800 | |||
4864ed8349 | |||
5b630571b3 | |||
06a2b217d4 | |||
2a687571c6 |
3
.gitignore
vendored
3
.gitignore
vendored
@ -3,3 +3,6 @@ fit_results_1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_201
|
||||
plot_1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_2017.11.27140704.pdf
|
||||
plot_1c5ca4b12ae2ddffc3960c1fe39a3cce35967ce23dbac57c010f450e796d01fd_2017.11.27140704.png
|
||||
naidis_fit.pdf
|
||||
__pycache__/Data.cpython-312.pyc
|
||||
__pycache__/Model.cpython-312.pyc
|
||||
__pycache__/fitter.cpython-312.pyc
|
||||
|
59
Model.py
59
Model.py
@ -33,46 +33,52 @@ class Model:
|
||||
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.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.Cai = 0.1 # uM, Cytoplasmic Ca2+ concentration fixed at point
|
||||
|
||||
self.LTRPN_tot: float = (
|
||||
70.0 # uM, Total myoplasmic troponin low-affinity site concentration
|
||||
70.0 # uM, Total myoplasmic troponin low-affinity site concentr.
|
||||
)
|
||||
self.HTRPN_tot: float = (
|
||||
140.0 # uM, Total myoplasmic troponin high-affinity site concentration
|
||||
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.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
|
||||
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-affinity sites
|
||||
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 concentration
|
||||
self.Km_CMDN: float = 0.238 # uM, Ca2 half-saturation constant for calmodulin
|
||||
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, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR
|
||||
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 the 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-saturation constant for SR Ca2+ -ATPase pump
|
||||
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
|
||||
)
|
||||
@ -84,7 +90,7 @@ class Model:
|
||||
)
|
||||
|
||||
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 inact.
|
||||
)
|
||||
self.K_pc_half: float = (
|
||||
20.0 # uM, Half-saturation constant for Ca2+-induced inactivation
|
||||
@ -97,7 +103,7 @@ class Model:
|
||||
|
||||
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
|
||||
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
|
||||
@ -117,9 +123,9 @@ class Model:
|
||||
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_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 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
|
||||
@ -160,7 +166,7 @@ class Model:
|
||||
|
||||
Kpcf = 13.0 * (1 - np.exp(-((V + 14.5) ** 2) / 100))
|
||||
|
||||
alpha = ( # dimensionless parameter for L-type Ca2+ channel closed states
|
||||
alpha = ( # dimensionless param. for L-type Ca2+ channel closed states
|
||||
0.4
|
||||
* np.exp((V + 12.0) / 10)
|
||||
* (
|
||||
@ -267,8 +273,10 @@ class Model:
|
||||
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
|
||||
-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
|
||||
@ -289,7 +297,8 @@ class Model:
|
||||
|
||||
"""
|
||||
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)))
|
||||
- (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
|
||||
@ -308,8 +317,6 @@ class Model:
|
||||
* ((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,
|
||||
@ -444,6 +451,8 @@ if __name__ == "__main__":
|
||||
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)
|
||||
states, V = model.solve(
|
||||
initial_values=initial_values, tspan=[0, 1000], dt=1.0, times=times
|
||||
)
|
||||
|
||||
model.plot_results(times, states, V)
|
||||
|
Binary file not shown.
@ -6,6 +6,45 @@ from Model import Model
|
||||
from fitter import Fitter
|
||||
from Data import Data
|
||||
|
||||
file = "ltcc_current.h5"
|
||||
|
||||
|
||||
def print_attrs(name, obj):
|
||||
print(f"\nAttributes for {name}:")
|
||||
for key, val in obj.attrs.items():
|
||||
print(f" {key}: {val}")
|
||||
|
||||
|
||||
with h5py.File(file, "r") as h5_file:
|
||||
h5_file.visititems(print_attrs)
|
||||
|
||||
# Dict to hold DFs 'sex', 'tag' & 'spid'
|
||||
dfs_by_sex_tag_spid = {}
|
||||
|
||||
with h5py.File(file, "r") as h5_file:
|
||||
for eid in h5_file.keys():
|
||||
attributes = h5_file[eid].attrs
|
||||
sex = attributes.get("sex")
|
||||
tag = attributes.get("tag")
|
||||
spid = attributes.get("spid")
|
||||
|
||||
key = f"{sex}_{tag}_{spid}"
|
||||
|
||||
if key not in dfs_by_sex_tag_spid:
|
||||
dfs_by_sex_tag_spid[key] = pd.DataFrame()
|
||||
|
||||
row_data = {"experiment_id": eid, "sex": sex, "tag": tag, "spid": spid}
|
||||
temp_df = pd.DataFrame([row_data])
|
||||
|
||||
dfs_by_sex_tag_spid[key] = pd.concat(
|
||||
[dfs_by_sex_tag_spid[key], temp_df], ignore_index=True
|
||||
)
|
||||
|
||||
for key, df in dfs_by_sex_tag_spid.items():
|
||||
print(f"DataFrame for {key}:")
|
||||
print(df)
|
||||
print()
|
||||
|
||||
|
||||
def fit_data():
|
||||
filename = "ltcc_current.h5"
|
||||
|
16
fitter.py
16
fitter.py
@ -61,16 +61,13 @@ class Fitter:
|
||||
model.K_pc_half = K_pc_half
|
||||
model.tau_xfer = tau_xfer
|
||||
|
||||
model.solve(
|
||||
states, V = model.solve(
|
||||
times=self.time_points,
|
||||
tspan=self.tspan,
|
||||
dt=self.dt,
|
||||
initial_values=self.initial_values,
|
||||
)
|
||||
|
||||
states = self.initial_values
|
||||
t = self.time_points
|
||||
V = model.mem_potential(t)
|
||||
_calc_curr = model.calculated_current(states, V)
|
||||
calculated_current = (
|
||||
self.convolve_current(_calc_curr, tau=tau_RC)[self.current_time_indecies]
|
||||
@ -100,6 +97,7 @@ class Fitter:
|
||||
def optimize(self, init_parameters=None):
|
||||
t0 = time.time()
|
||||
self.iteration = 0
|
||||
|
||||
if init_parameters is None:
|
||||
m = self.model()
|
||||
K_pc_half = m.K_pc_half
|
||||
@ -118,7 +116,15 @@ class Fitter:
|
||||
(10, 100, 100, 100, 10, 10),
|
||||
)
|
||||
|
||||
res = least_squares(self.cost_func, init_parameters, bounds=bounds, xtol=1e-10)
|
||||
res = least_squares(
|
||||
self.cost_func,
|
||||
init_parameters,
|
||||
bounds=bounds,
|
||||
xtol=1e-10,
|
||||
ftol=1e-10,
|
||||
# verbose=2,
|
||||
)
|
||||
|
||||
print()
|
||||
print(" Parameters: [gGaL, ECal, K_pc_half, tau_xfer, tau_RC, offset]")
|
||||
print(" Initial:", init_parameters.tolist())
|
||||
|
136
multiple_experiment_fitter.py
Normal file
136
multiple_experiment_fitter.py
Normal file
@ -0,0 +1,136 @@
|
||||
import h5py
|
||||
import pandas as pd
|
||||
import numpy as np
|
||||
import re
|
||||
|
||||
import matplotlib as plt
|
||||
|
||||
from Model import Model
|
||||
from fitter import Fitter
|
||||
from Data import Data
|
||||
|
||||
file = "ltcc_current.h5"
|
||||
|
||||
dfs_by_sex_tag_spid = {}
|
||||
|
||||
|
||||
def print_attrs(name, obj):
|
||||
# print(f"\nAttributes for {name}:")
|
||||
# for key, val in obj.attrs.items():
|
||||
# print(f" {key}: {val}")
|
||||
pass
|
||||
|
||||
|
||||
with h5py.File(file, "r") as h5_file:
|
||||
for eid in h5_file.keys():
|
||||
attributes = h5_file[eid].attrs
|
||||
sex = attributes.get("sex")
|
||||
tag = attributes.get("tag")
|
||||
spid = attributes.get("spid")
|
||||
|
||||
key = f"{sex}_{tag}_{spid}"
|
||||
|
||||
if key not in dfs_by_sex_tag_spid:
|
||||
dfs_by_sex_tag_spid[key] = pd.DataFrame()
|
||||
|
||||
row_data = {"experiment_id": eid, "sex": sex, "tag": tag, "spid": spid}
|
||||
temp_df = pd.DataFrame([row_data])
|
||||
|
||||
dfs_by_sex_tag_spid[key] = pd.concat(
|
||||
[dfs_by_sex_tag_spid[key], temp_df], ignore_index=True
|
||||
)
|
||||
|
||||
# for key, df in dfs_by_sex_tag_spid.items():
|
||||
# print(f"DataFrame for {key}:")
|
||||
# print(df)
|
||||
# print()
|
||||
|
||||
|
||||
def fit_and_plot_dataframes(dfs_by_sex_tag_spid):
|
||||
for key, df in dfs_by_sex_tag_spid.items():
|
||||
print(f"Fitting and plotting data for {key}...")
|
||||
|
||||
combined_current = []
|
||||
combined_time = []
|
||||
|
||||
for eid in df["experiment_id"].tolist():
|
||||
data = Data(file, group_key=eid)
|
||||
combined_current.append(data.current)
|
||||
combined_time.append(data.current_t)
|
||||
|
||||
combined_current = np.concatenate(combined_current)
|
||||
combined_time = np.concatenate(combined_time)
|
||||
|
||||
# Sort by time for consistency - inspired by data
|
||||
sorted_indices = np.argsort(combined_time)
|
||||
combined_time = combined_time[sorted_indices]
|
||||
combined_current = combined_current[sorted_indices]
|
||||
|
||||
combined_data = Data(file, group_key=eid)
|
||||
combined_data.current = combined_current
|
||||
combined_data.current_t = combined_time
|
||||
|
||||
fit = Fitter(Model, combined_data)
|
||||
fit.optimize()
|
||||
res, fig = fit.optimize()
|
||||
|
||||
plt.figure()
|
||||
plt.title(f"Fit results for {key}")
|
||||
|
||||
for i, eid in enumerate(df["experiment_id"].tolist()):
|
||||
plt.plot(combined_time, combined_current,
|
||||
label=f"Experiment {eid}")
|
||||
|
||||
plt.plot(
|
||||
combined_data.current_t, combined_data.current, "k-",
|
||||
label="Combined Fit"
|
||||
)
|
||||
|
||||
plt.legend()
|
||||
plt.xlabel("Time")
|
||||
plt.ylabel("Current")
|
||||
|
||||
key_cleaned = re.sub(r"[^\w.-]", "", key)
|
||||
plt.savefig(f"combined_plot_{key_cleaned}.png")
|
||||
plt.savefig(f"combined_plot_{key_cleaned}.pdf")
|
||||
plt.close()
|
||||
|
||||
fit_hist = pd.DataFrame.from_dict(fit.fit_results, orient="index").T
|
||||
fit_hist.index.name = "Iterations"
|
||||
|
||||
res_filename = f"combined_fit_results_{key}.csv"
|
||||
res_filename = res_filename.replace(" ", "_").replace(":", "-")
|
||||
fit_hist.to_csv(res_filename, index=True)
|
||||
|
||||
print(f"Finished fitting for {key}.")
|
||||
|
||||
|
||||
fit_and_plot_dataframes(dfs_by_sex_tag_spid)
|
||||
|
||||
"""
|
||||
def fit_data():
|
||||
filename = "ltcc_current.h5"
|
||||
with h5py.File(filename, "r") as h5:
|
||||
eids = list(h5.keys())
|
||||
|
||||
for eid in eids:
|
||||
data = Data(filename, group_key=eid)
|
||||
fit = Fitter(Model, data)
|
||||
fit.optimize()
|
||||
res, fig = fit.optimize()
|
||||
|
||||
fit_hist = pd.DataFrame.from_dict(fit.fit_results, orient="index").T
|
||||
fit_hist.index.name = "Iterations"
|
||||
|
||||
res_filename = f"fit_results_{eid}.csv"
|
||||
res_filename = res_filename.replace(" ", "_").replace(":", "-")
|
||||
fit_hist.to_csv(res_filename, index=True)
|
||||
|
||||
eid_cleaned = re.sub(r"[^w.-]", "", eid) # Eemaldab kõik eritähed
|
||||
fig.savefig(f"plot_{eid_cleaned}.png")
|
||||
fig.savefig(f"plot_{eid_cleaned}.pdf")
|
||||
|
||||
|
||||
fit_data()
|
||||
|
||||
"""
|
Reference in New Issue
Block a user