Add scaling factors for omecamtiv mecarbil

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
         )

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)

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")

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,