Merge pull request #278 from biomass-dev/dev

Release v0.14.0

Author: Hiroaki Imoto

biomass/construction/reaction_rules.py

@@ -1811,13 +1811,17 @@ def user_defined(self, line_num: int, line: str) -> None:
         Examples
         --------
         >>> '@rxn Reactant --> Product: define rate equation here'
+        >>> '@rxn A --> 0: p[Vdeg] * u[A] * u[B] / (p[Kdeg] + u[A])'
+        >>> '@rxn E + S --> E + P: p[Vmax] * u[E] * u[S] / (p[Km] + u[S])'
 
         Notes
         -----
         * Use p[xxx] and u[xxx] for describing parameters and species, respectively.
 
         * Use '0' or '∅' for degradation/creation to/from nothing.
 
+        * Multiple reactants and products can be described by using '+'.
+
         * Differential equation
             .. math::
 
@@ -1833,52 +1837,79 @@ def user_defined(self, line_num: int, line: str) -> None:
         rate_equation = description[1].strip()
         for arrow in self.fwd_arrows:
             if arrow in balance:
-                reactant, product = balance.split(arrow)
-                if reactant in self.nothing:
-                    self._set_species(product.strip())
-                elif product in self.nothing:
-                    self._set_species(product.strip())
+                two_species = balance.split(arrow)
+                # reactant, product = balance.split(arrow)
+                if len(two_species[0].split(" ")) > 1 and "+" not in two_species[0]:
+                    reactants = None
+                else:
+                    reactants = [
+                        reactant
+                        for reactant in two_species[0].replace(" ", "").split("+")
+                        if reactant not in self.nothing
+                    ]
+                if len(two_species[1].split(" ")) > 1 and "+" not in two_species[1]:
+                    products = None
                 else:
-                    self._set_species(reactant.strip(), product.strip())
+                    products = [
+                        product
+                        for product in two_species[1].replace(" ", "").split("+")
+                        if product not in self.nothing
+                    ]
                 break
         else:
             raise ArrowError(f"line{line_num:d}: Use one of {', '.join(self.fwd_arrows)}.")
+        if reactants is None or products is None:
+            raise DetectionError(f"Unregistered words in line{line_num:d}: {line}")
+        if len(reactants) == 0:
+            self._set_species(*products)
+        elif len(products) == 0:
+            self._set_species(*reactants)
+        else:
+            self._set_species(*reactants, *products)
         rate_equation = (
             rate_equation.replace("p[", "x[C.").replace("u[", "y[V.").replace("^", "**")
         )
         self.reactions.append(f"v[{line_num:d}] = " + rate_equation.strip())
+        modulators = set(reactants) & set(products)
+        for modulator in modulators:
+            reactants.remove(modulator)
+            products.remove(modulator)
         modulators = (
             *list(
                 set(
                     [
                         ent
                         for ent in re.findall(r"(?<=\[V.)(.+?)(?=\])", rate_equation)
-                        if ent not in [reactant, product]
+                        if ent not in [*reactants, *products]
                     ]
                 )
+                | modulators
             ),
         )
         self.kinetics.append(
             KineticInfo(
-                () if reactant in self.nothing else (reactant,),
-                () if product in self.nothing else (product,),
-                () if modulators is None else (modulators),
+                tuple(reactants),
+                tuple(products),
+                tuple(modulators),
                 rate_equation.replace("x[C.", "").replace("y[V.", "").replace("]", ""),
             )
         )
-        counter_reactant = 0
-        counter_product = 0
-        for i, eq in enumerate(self.differential_equations):
-            if f"dydt[V.{reactant}]" in eq and reactant not in self.nothing:
-                counter_reactant += 1
-                self.differential_equations[i] = eq + f" - v[{line_num:d}]"
-            elif f"dydt[V.{product}]" in eq and product not in self.nothing:
-                counter_product += 1
-                self.differential_equations[i] = eq + f" + v[{line_num:d}]"
-        if counter_reactant == 0 and reactant not in self.nothing:
-            self.differential_equations.append(f"dydt[V.{reactant}] = - v[{line_num:d}]")
-        if counter_product == 0 and product not in self.nothing:
-            self.differential_equations.append(f"dydt[V.{product}] = + v[{line_num:d}]")
+        for reactant in reactants:
+            counter_reactant = 0
+            for i, eq in enumerate(self.differential_equations):
+                if f"dydt[V.{reactant}]" in eq:
+                    counter_reactant += 1
+                    self.differential_equations[i] = eq + f" - v[{line_num:d}]"
+            if counter_reactant == 0:
+                self.differential_equations.append(f"dydt[V.{reactant}] = - v[{line_num:d}]")
+        for product in products:
+            counter_product = 0
+            for i, eq in enumerate(self.differential_equations):
+                if f"dydt[V.{product}]" in eq:
+                    counter_product += 1
+                    self.differential_equations[i] = eq + f" + v[{line_num:d}]"
+            if counter_product == 0:
+                self.differential_equations.append(f"dydt[V.{product}] = + v[{line_num:d}]")
 
     def _extract_event(self, line_num: int, line: str):
         # About biochemical event

biomass/construction/template/observable.py

@@ -1,4 +1,4 @@
-from typing import List
+from typing import Dict, List
 
 import numpy as np
 
@@ -71,7 +71,7 @@ def simulate(self, x: list, y0: list, _perturbation: dict = {}):
     def set_data(self) -> None:
         pass
 
-    def get_timepoint(self, obs_name: str) -> List[int]:
+    def get_timepoint(self, obs_name: str) -> Dict[str, List[int]]:
         if obs_name in self.obs_names:
-            return []
+            return {condition: [] for condition in self.conditions}
         assert False

biomass/construction/template/problem.py

@@ -36,7 +36,7 @@ def _diff_sim_and_exp(sim_matrix, exp_dict, exp_timepoint, conditions, sim_norm_
 
         for idx, condition in enumerate(conditions):
             if condition in exp_dict.keys():
-                sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint))])
+                sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint[condition]))])
                 exp_val.extend(exp_dict[condition])
 
         return np.array(sim_val) / sim_norm_max, np.array(exp_val)

biomass/dynamics/temporal_dynamics.py

@@ -1,7 +1,7 @@
 import os
 from dataclasses import dataclass
 from math import isnan
-from typing import List, Optional
+from typing import Dict, List, Optional
 
 import numpy as np
 from matplotlib import pyplot as plt
@@ -93,14 +93,14 @@ def plot_timecourse(
                 if (
                     viz_type == "experiment" or singleplotting[i].exp_data
                 ) and self.model.problem.experiments[i] is not None:
-                    exp_t = self.model.problem.get_timepoint(obs_name)
+                    exp_t_dict = self.model.problem.get_timepoint(obs_name)
                     if self.model.problem.error_bars[i] is not None:
                         self._plot_experimental_data_with_error_bars(
-                            viz_type, exp_t, obs_name, mode, singleplotting, multiplotting
+                            viz_type, exp_t_dict, obs_name, mode, singleplotting, multiplotting
                         )
                     else:
                         self._plot_experimental_data_without_error_bars(
-                            viz_type, exp_t, obs_name, mode, singleplotting, multiplotting
+                            viz_type, exp_t_dict, obs_name, mode, singleplotting, multiplotting
                         )
                 if mode == 0:
                     self._save_mode_0(obs_name, singleplotting, viz_type)
@@ -360,7 +360,7 @@ def _plot_simulations(
     def _plot_experimental_data_with_error_bars(
         self,
         viz_type: str,
-        exp_t: Optional[List[int]],
+        exp_t_dict: Optional[Dict[str, List[int]]],
         obs_name: str,
         mode: int,
         singleplotting: List[SingleObservable],
@@ -377,6 +377,7 @@ def _plot_experimental_data_with_error_bars(
                 or (mode == 1 and condition == multiplotting.condition)
             ):
                 try:
+                    exp_t = exp_t_dict[condition]
                     exp_data = plt.errorbar(
                         np.array(exp_t) / singleplotting[i].divided_by,
                         self.model.problem.experiments[i][condition],
@@ -421,7 +422,7 @@ def _plot_experimental_data_with_error_bars(
     def _plot_experimental_data_without_error_bars(
         self,
         viz_type: str,
-        exp_t: Optional[List[int]],
+        exp_t_dict: Optional[Dict[str, List[int]]],
         obs_name: str,
         mode: int,
         singleplotting: List[SingleObservable],
@@ -438,6 +439,7 @@ def _plot_experimental_data_without_error_bars(
                 or (mode == 1 and condition == multiplotting.condition)
             ):
                 try:
+                    exp_t = exp_t_dict[condition]
                     plt.plot(
                         np.array(exp_t) / singleplotting[i].divided_by,
                         self.model.problem.experiments[i][condition],

biomass/models/Nakakuki_Cell_2010/observable.py

@@ -1,4 +1,4 @@
-from typing import List
+from typing import Dict, List
 
 import numpy as np
 
@@ -201,19 +201,23 @@ def set_data(self):
             "HRG": [sd / np.sqrt(3) for sd in [0.027, 0.059, 0.094, 0.124, 0.113, 0.108]],
         }
 
-    @staticmethod
-    def get_timepoint(obs_name) -> List[int]:
+    def get_timepoint(self, obs_name) -> Dict[str, List[int]]:
         if obs_name in [
             "Phosphorylated_MEKc",
             "Phosphorylated_ERKc",
             "Phosphorylated_RSKw",
             "Phosphorylated_cFos",
         ]:
-            return [0, 300, 600, 900, 1800, 2700, 3600, 5400]  # (Unit: sec.)
+            return {
+                condition: [0, 300, 600, 900, 1800, 2700, 3600, 5400]  # (Unit: sec.)
+                for condition in self.conditions
+            }
         elif obs_name == "Phosphorylated_CREBw":
-            return [0, 600, 1800, 3600, 5400]
+            return {condition: [0, 600, 1800, 3600, 5400] for condition in self.conditions}
         elif obs_name == "cfos_mRNA":
-            return [0, 600, 1200, 1800, 2700, 3600, 5400]
+            return {
+                condition: [0, 600, 1200, 1800, 2700, 3600, 5400] for condition in self.conditions
+            }
         elif obs_name in ["cFos_Protein", "dusp_mRNA"]:
-            return [0, 900, 1800, 2700, 3600, 5400]
+            return {condition: [0, 900, 1800, 2700, 3600, 5400] for condition in self.conditions}
         assert False

biomass/models/Nakakuki_Cell_2010/problem.py

@@ -36,7 +36,7 @@ def _diff_sim_and_exp(sim_matrix, exp_dict, exp_timepoint, conditions, sim_norm_
 
         for idx, condition in enumerate(conditions):
             if condition in exp_dict.keys():
-                sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint))])
+                sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint[condition]))])
                 exp_val.extend(exp_dict[condition])
 
         return np.array(sim_val) / sim_norm_max, np.array(exp_val)

biomass/models/_copy.py

@@ -19,6 +19,7 @@ def copy_to_current(
         "tgfb_smad",
         "g1s_transition",
         "prolif_quies",
+        "jak_stat_pathway",
     ]
 ) -> None:
     """Copy an example model to the current working directory.
@@ -46,6 +47,7 @@ def copy_to_current(
             "tgfb_smad",
             "g1s_transition",
             "prolif_quies",
+            "jak_stat_pathway",
         ]
     ):
         raise ValueError(f"model_name must be one of [{', '.join(available_models)}].")

biomass/models/circadian_clock/observable.py

@@ -79,5 +79,5 @@ def set_data(self):
 
     def get_timepoint(self, obs_name):
         if obs_name in self.obs_names:
-            return []
+            return {}
         assert False

biomass/models/circadian_clock/problem.py

@@ -36,7 +36,7 @@ def _diff_sim_and_exp(sim_matrix, exp_dict, exp_timepoint, conditions, sim_norm_
 
         for idx, condition in enumerate(conditions):
             if condition in exp_dict.keys():
-                sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint))])
+                sim_val.extend(sim_matrix[idx, list(map(int, exp_timepoint[condition]))])
                 exp_val.extend(exp_dict[condition])
 
         return np.array(sim_val) / sim_norm_max, np.array(exp_val)

biomass/models/g1s_transition/observable.py

@@ -82,5 +82,5 @@ def set_data(self):
 
     def get_timepoint(self, obs_name):
         if obs_name in self.obs_names:
-            return []
+            return {}
         assert False