diff --git a/src/simcardems/ORdmm_Land.py b/src/simcardems/ORdmm_Land.py index b59c856c..8a096604 100644 --- a/src/simcardems/ORdmm_Land.py +++ b/src/simcardems/ORdmm_Land.py @@ -192,6 +192,8 @@ def default_parameters(disease_state: str = "healthy") -> Dict[str, float]: ("scale_drug_IpCa", 1.0), ("scale_drug_Isacns", 1.0), ("scale_drug_Isack", 1.0), + ("scale_drug_kws", 1.0), + ("scale_drug_kuw", 1.0), # Population factors ("scale_popu_GNa", 1.0), ("scale_popu_GCaL", 1.0), @@ -996,6 +998,8 @@ def F(self, v, s, time=None): scale_drug_IpCa = self._parameters["scale_drug_IpCa"] scale_drug_Isacns = self._parameters["scale_drug_Isacns"] scale_drug_Isack = self._parameters["scale_drug_Isack"] + scale_drug_kws = self._parameters["scale_drug_kws"] + scale_drug_kuw = self._parameters["scale_drug_kuw"] # Population factors scale_popu_GNa = self._parameters["scale_popu_GNa"] @@ -1585,10 +1589,12 @@ def F(self, v, s, time=None): zetas2 = (-1 - Zetas) * ufl.conditional(ufl.lt(Zetas, -1), 1, 0) gammasu = gammas * Max(zetas1, zetas2) - F_expressions[39] = kws * scale_popu_kws * XW - XS * gammasu - XS * ksu + F_expressions[39] = ( + kws * scale_drug_kws * scale_popu_kws * XW - XS * gammasu - XS * ksu + ) F_expressions[40] = ( - kuw * scale_popu_kuw * XU - - kws * scale_popu_kws * XW + kuw * scale_drug_kuw * scale_popu_kuw * XU + - kws * scale_drug_kws * scale_popu_kws * XW - XW * gammawu - XW * kwu ) diff --git a/src/simcardems/ep_model.py b/src/simcardems/ep_model.py index 7c5cf99c..57e370fb 100644 --- a/src/simcardems/ep_model.py +++ b/src/simcardems/ep_model.py @@ -54,8 +54,16 @@ def setup_ep_model(cellmodel, mesh): # Mark stimulation region defined as [0, L]^3 S1_marker = 1 - S1_markers = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()) - S1_markers.set_all(S1_marker) # Mark the whole mesh + + setstim = "stimall" + if setstim == "stimpart": # Stimpart + L = 2.0 + S1_subdomain = dolfin.CompiledSubDomain("x[0] <= L + DOLFIN_EPS", L=L) + S1_markers = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()) + S1_subdomain.mark(S1_markers, S1_marker) + elif setstim == "stimall": # Stimall + S1_markers = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()) + S1_markers.set_all(S1_marker) # Mark the whole mesh # Define stimulation (NB: region of interest carried by the mesh # and assumptions in cbcbeat) diff --git a/src/simcardems/postprocess.py b/src/simcardems/postprocess.py index 255b9a2b..253fb434 100644 --- a/src/simcardems/postprocess.py +++ b/src/simcardems/postprocess.py @@ -282,7 +282,8 @@ def plot_peaks(fname, data, threshold): def plot_state_traces(results_file): - fig, ax = plt.subplots(2, 2, figsize=(10, 8), sharex=True) + fig, ax = plt.subplots(2, 2, sharex=True) + fig_mech, ax_mech = plt.subplots(2, 2, sharex=True) results_file = Path(results_file) if not results_file.is_file(): raise FileNotFoundError(f"File {results_file} does not exist") @@ -292,7 +293,10 @@ def plot_state_traces(results_file): loader = DataLoader(results_file) bnd = {"ep": Boundary(loader.ep_mesh), "mechanics": Boundary(loader.mech_mesh)} - all_names = {"mechanics": ["lmbda", "Ta"], "ep": ["V", "Ca"]} + all_names = { + "mechanics": ["lmbda", "Ta", "zetas", "zetaw", "XS", "XW"], + "ep": ["V", "Ca", "CaTrpn", "TmB"], + } values = { group: {name: np.zeros(len(loader.time_stamps)) for name in names} @@ -348,6 +352,23 @@ def plot_state_traces(results_file): ) fig.savefig(outdir.joinpath("state_traces.png"), dpi=300) + ax_mech[0, 0].plot(times, values["mechanics"]["zetas"]) + ax_mech[0, 0].plot(times, values["mechanics"]["zetaw"]) + ax_mech[1, 0].plot(times, values["mechanics"]["XS"]) + ax_mech[1, 0].plot(times, values["mechanics"]["XW"]) + ax_mech[0, 1].plot(times, values["ep"]["CaTrpn"]) + ax_mech[1, 1].plot(times, values["ep"]["TmB"]) + + ax_mech[0, 0].set_title("zeta(s/w)") + ax_mech[1, 0].set_title("X(S/W)") + ax_mech[0, 1].set_title("CaTrpn") + ax_mech[1, 1].set_title("TmB") + ax_mech[1, 0].set_xlabel("Time (ms)") + ax_mech[1, 1].set_xlabel("Time (ms)") + for axi in ax_mech.flatten(): + axi.grid() + fig_mech.savefig(outdir.joinpath("state_traces_mechanics.png"), dpi=300) + def make_xdmffiles(results_file): @@ -476,9 +497,11 @@ def get_biomarkers(dict, outdir, num_models): time[np.where(u == np.min(u))[0][0]] % 1000 ) - ax.plot(PoMm, biomarker_dict[f"m{PoMm}"]["APD90"], "*") + if num_models > 1: + ax.plot(PoMm, biomarker_dict[f"m{PoMm}"]["APD90"], "*") - fig.savefig(outdir.joinpath("APD90_permodel.png"), dpi=300) + if num_models > 1: + fig.savefig(outdir.joinpath("APD90_permodel.png"), dpi=300) with open(outdir.joinpath("biomarkers_PoMcontrol.json"), "w") as f: json.dump(biomarker_dict, f) @@ -511,7 +534,10 @@ def save_popu_json(population_folder, num_models): dict = {} for PoMm in range(1, num_models + 1): print(f"Analyzing model {PoMm}") - results_file = population_folder.joinpath(f"m{PoMm}/results.h5") + if num_models == 1: + results_file = population_folder.joinpath("results.h5") + else: + results_file = population_folder.joinpath(f"m{PoMm}/results.h5") if not results_file.is_file(): raise FileNotFoundError(f"File {results_file} does not exist") diff --git a/src/simcardems/setup_models.py b/src/simcardems/setup_models.py index 8e915fb0..3181d6bb 100644 --- a/src/simcardems/setup_models.py +++ b/src/simcardems/setup_models.py @@ -391,6 +391,9 @@ def _setup_assigners(self): ) self._assign_displacement() + self.CaTrpn, self.CaTrpn_assigner = utils.setup_assigner(self._vs, 42) + self.TmB, self.TmB_assigner = utils.setup_assigner(self._vs, 43) + def _assign_displacement(self): self._u_assigner.assign( self._u, @@ -400,6 +403,8 @@ def _assign_displacement(self): def _assign_ep(self): self._v_assigner.assign(self._v, utils.sub_function(self._vs, 0)) self._Ca_assigner.assign(self._Ca, utils.sub_function(self._vs, 45)) + self.CaTrpn_assigner.assign(self.CaTrpn, utils.sub_function(self._vs, 42)) + self.TmB_assigner.assign(self.TmB, utils.sub_function(self._vs, 43)) def store(self): # Assign u, v and Ca for postprocessing @@ -422,6 +427,12 @@ def _setup_datacollector(self): ("ep", "Ca", self._Ca), ("mechanics", "lmbda", self.coupling.lmbda_mech), ("mechanics", "Ta", self.mech_heart.material.active.Ta_current), + ("ep", "CaTrpn", self.CaTrpn), + ("ep", "TmB", self.TmB), + ("mechanics", "zetas", self.coupling.Zetas_mech), + ("mechanics", "zetaw", self.coupling.Zetaw_mech), + ("mechanics", "XS", self.coupling.XS_mech), + ("mechanics", "XW", self.coupling.XW_mech), ]: self.collector.register(group, name, f) @@ -436,7 +447,7 @@ def _solve_mechanics_now(self) -> bool: # dt for the mechanics model should not be larger than 1 ms dt = self._t - self.mech_heart.material.active.t - return (XS_norm + XW_norm >= 0.1) or dt > 0.990 + return (XS_norm + XW_norm >= 0.05) or dt > 0.990 def _pre_mechanics_solve(self) -> None: self._preXS_assigner.assign(self._pre_XS, utils.sub_function(self._vs, 40)) @@ -497,7 +508,9 @@ def solve( # Store state every 5 beats if i > 0 and i % int(5000 / self._dt) == 0: io.save_state( - self._state_path.parent.joinpath(f"state_{int(i*self._dt/1000)}beat.h5"), + self._state_path.parent.joinpath( + f"state_{int(i*self._dt/1000)}beat.h5", + ), solver=self.ep_solver, mech_heart=self.mech_heart, coupling=self.coupling,