abstochkin.de_calcs

Perform deterministic calculations on a set of processes. Construct the ordinary differential equations (ODEs) describing the system and obtain a numerical solution.

View Source

  1"""
  2Perform deterministic calculations on a set of processes.
  3Construct the ordinary differential equations (ODEs)
  4describing the system and obtain a numerical solution.
  5"""
  6
  7#  Copyright (c) 2024-2025, Alex Plakantonakis.
  8#
  9#  This program is free software: you can redistribute it and/or modify
 10#  it under the terms of the GNU General Public License as published by
 11#  the Free Software Foundation, either version 3 of the License, or
 12#  (at your option) any later version.
 13#
 14#  This program is distributed in the hope that it will be useful,
 15#  but WITHOUT ANY WARRANTY; without even the implied warranty of
 16#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 17#  GNU General Public License for more details.
 18#
 19#  You should have received a copy of the GNU General Public License
 20#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 21
 22from numpy import mean
 23from scipy.integrate import solve_ivp
 24from sympy import Add, Mul, Pow, sympify, lambdify, symbols
 25# We want all single-letter and Greek-letter variables to be symbols.
 26# We can then use the clashing-symbols dictionaries that have been defined
 27# as private variables in `_clash` (which includes both single and
 28# multi-letter names that are defined in `sympy.abc`).
 29from sympy.abc import _clash
 30
 31from .process import update_all_species, MichaelisMentenProcess, RegulatedProcess, \
 32    RegulatedMichaelisMentenProcess
 33
 34
 35class DEcalcs:
 36    """
 37    Perform deterministic calculations given the processes specified
 38    in an AbStochKin simulation object.
 39
 40    Attributes
 41    ----------
 42    p0 : dict[str: int]
 43        Dictionary specifying the initial population sizes of all
 44        species in the given processes.
 45    t_min : float or int
 46        Numerical value of the start of simulated time in the units
 47        specified in the class attribute `time_unit`.
 48    t_max : float or int
 49        Numerical value of the end of simulated time in the units
 50        specified in the class attribute `time_unit`.
 51    processes : list of Process objects
 52        A list of all the processes that define the system of ODEs.
 53        Each process is a `Process` object.
 54    ode_method: str
 55        Method for the ODE solver to use.
 56    time_unit: str
 57        A string of the time unit to be used for describing the
 58        kinetics of the given processes.
 59    """
 60
 61    def __init__(self,
 62                 p0: dict,
 63                 t_min: int | float,
 64                 t_max: int | float,
 65                 processes: list,
 66                 ode_method: str,
 67                 time_unit: str):
 68        self.p0 = p0
 69        self.t_min = t_min
 70        self.t_max = t_max
 71        self.processes = processes
 72        self.ode_method = ode_method
 73        self.time_unit = time_unit
 74
 75        # *** For computing and storing the solution to the system of ODEs. ***
 76        self._all_species, self._procs_by_reactant, self._procs_by_product = update_all_species(
 77            tuple(self.processes))
 78        self.odes = dict()  # store species-specific ODEs
 79        self.odes_sol = None  # numerical solution of species trajectories
 80
 81        self.species_with_ode = list()
 82        self.jac = None  # Jacobian matrix
 83        self.fixed_pts = dict()  # dictionary of fixed points
 84
 85        self.setup_ODEs()
 86        # self.get_fixed_pts()
 87
 88    def setup_ODEs(self, agent_based=True):
 89        """
 90        Set up the system of ODEs to be used for computing the
 91        deterministic trajectory of all the species in the given processes.
 92        The equations consist of sympy objects.
 93
 94        Parameters
 95        ----------
 96        agent_based : bool, default: True
 97            If True, set up the agent-based (or microscopic) form of ODEs.
 98            For instance, for the process `2X -> Y`, the ODE for species
 99            `X` would include an `X(X - 1)` term instead
100            of `X^2` (the canonical form). If False, the canonical form
101            of the ODEs is constructed.
102
103        Notes
104        -----
105        The rate constant, `k`, for a given process is taken to be the mean
106        of `k`, unless `k` was defined to be normally-distributed,
107        in which case `k` is a 2-tuple and `k[0]` is the specified mean.
108
109        Implementation note: Building up the ODE expressions by separately iterating
110        over the processes for products and reactants. This is to properly
111        handle 0th order processes for product species. For example, for the
112        0th order process ' --> X' with rate constant k1, the ODE is dX/dt = k1.
113        """
114        # Set up ODEs for product species
115        for spec, procs in self._procs_by_product.items():
116            terms = []
117            for proc in procs:
118                k_val = proc.k[0] if isinstance(proc.k, tuple) else mean(proc.k)
119
120                if agent_based:
121                    temp = [Mul(*[sympify(name, locals=_clash) - i for i in
122                                  range(order)]) if name != '' else 1
123                            for name, order in proc.reactants.items()]
124                else:
125                    temp = [Pow(sympify(name, locals=_clash), order) if name != '' else 1
126                            for name, order in proc.reactants.items()]
127
128                new_term = Mul(proc.products[spec], k_val, *temp)
129                terms.append(new_term * self.get_term_multiplier(proc))
130
131            if spec in self.odes:
132                self.odes[spec] += Add(*terms)
133            else:
134                self.odes[spec] = Add(*terms)
135
136        # Set up ODEs for reactant species
137        for spec, procs in self._procs_by_reactant.items():
138            terms = []
139            for proc in procs:
140                k_val = proc.k[0] if isinstance(proc.k, tuple) else mean(proc.k)
141
142                if agent_based:
143                    temp = [Mul(*[sympify(name, locals=_clash) - i for i in range(order)]) for
144                            name, order in proc.reactants.items()]
145                else:
146                    temp = [Pow(sympify(name, locals=_clash), order) for name, order in
147                            proc.reactants.items()]
148
149                new_term = Mul(proc.reactants[spec], -1 * k_val, *temp)
150                terms.append(new_term * self.get_term_multiplier(proc))
151
152            if spec in self.odes:
153                self.odes[spec] += Add(*terms)
154            else:
155                self.odes[spec] = Add(*terms)
156
157        self.species_with_ode = [sp for sp in self.odes.keys()]
158
159        # Handle species whose population remains constant (e.g., a catalyst)
160        if len(self._all_species) != len(self.species_with_ode):
161            missing_species = self._all_species - set(self.species_with_ode)
162            for m in missing_species:
163                self.odes[m] = sympify(0)  # dm/dt = 0
164                self.species_with_ode.append(m)
165
166    @staticmethod
167    def get_term_multiplier(proc):
168        """
169        Generate the multiplicative term (or terms) needed for generating the
170        correct algebraic expressions for specialized processes
171        (such as Michaelis-Menten and regulated processes).
172        """
173        if isinstance(proc, MichaelisMentenProcess):
174            Km_val = proc.Km[0] if isinstance(proc.Km, tuple) else mean(proc.Km)
175            c = symbols(proc.catalyst)
176            return c / (Km_val + sympify(proc.reacts_[0], locals=_clash))
177        elif issubclass(type(proc), RegulatedProcess):
178            if isinstance(proc.regulating_species, list):
179                term = 1
180                for sp, a, Kval, hc in zip(proc.regulating_species, proc.alpha, proc.K50, proc.nH):
181                    K50_val = Kval[0] if isinstance(Kval, tuple) else mean(Kval)
182                    ratio = (symbols(sp) / K50_val) ** hc
183                    term *= (1 + a * ratio) / (1 + ratio)
184            else:
185                K50_val = proc.K50[0] if isinstance(proc.K50, tuple) else mean(proc.K50)
186                ratio = (symbols(proc.regulating_species) / K50_val) ** proc.nH
187                term = (1 + proc.alpha * ratio) / (1 + ratio)
188
189            if isinstance(proc, RegulatedMichaelisMentenProcess):
190                Km_val = proc.Km[0] if isinstance(proc.Km, tuple) else mean(proc.Km)
191                c = symbols(proc.catalyst)
192                term *= c / (Km_val + sympify(proc.reacts_[0], locals=_clash))
193
194            return term
195        else:
196            return 1
197
198    def solve_ODEs(self):
199        """
200        Solve system of ordinary differential equations (ODEs).
201
202        Notes
203        -----
204        Using the solver `scipy.integrate.solve_ivp`,
205        whose method can be one of the following:
206
207        - RK45 : Explicit Runge-Kutta method of order 5(4).
208        - RK23 : Explicit Runge-Kutta method of order 3(2).
209        - DOP853 : Explicit Runge-Kutta method of order 8.
210        - Radau : Implicit Runge-Kutta method of the Radau IIA family of order 5
211        - BDF : Implicit multistep variable-order (1 to 5) method based on a
212                 backward differentiation formula for the derivative approximation.
213        - LSODA : Adams/BDF method with automatic stiffness detection and switching.
214
215        Documentation
216        -------------
217        https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
218        """
219        # Converting ODE expressions: sympy -> scipy/numpy
220        _odes_lambdified = lambdify([sympify(sp, locals=_clash) for sp in self.odes.keys()],
221                                    list(self.odes.values()))
222
223        # Converting initial population values to a list to ensure that the order of
224        # species-specific ODEs and initial values are correctly ordered.
225        p0_list = [self.p0[sp] for sp in self.odes.keys()]
226        self.odes_sol = solve_ivp(lambda t, S: _odes_lambdified(*S),
227                                  t_span=[self.t_min, self.t_max],
228                                  y0=p0_list,
229                                  method=self.ode_method,
230                                  t_eval=None,  # Specify points where the solution is desired
231                                  dense_output=True)  # Compute a continuous solution
232
233    # def get_fixed_pts(self):
234    #     """
235    #     Not currently implemented.
236    #     """
237    #     pass

class DEcalcs: View Source

 36class DEcalcs:
 37    """
 38    Perform deterministic calculations given the processes specified
 39    in an AbStochKin simulation object.
 40
 41    Attributes
 42    ----------
 43    p0 : dict[str: int]
 44        Dictionary specifying the initial population sizes of all
 45        species in the given processes.
 46    t_min : float or int
 47        Numerical value of the start of simulated time in the units
 48        specified in the class attribute `time_unit`.
 49    t_max : float or int
 50        Numerical value of the end of simulated time in the units
 51        specified in the class attribute `time_unit`.
 52    processes : list of Process objects
 53        A list of all the processes that define the system of ODEs.
 54        Each process is a `Process` object.
 55    ode_method: str
 56        Method for the ODE solver to use.
 57    time_unit: str
 58        A string of the time unit to be used for describing the
 59        kinetics of the given processes.
 60    """
 61
 62    def __init__(self,
 63                 p0: dict,
 64                 t_min: int | float,
 65                 t_max: int | float,
 66                 processes: list,
 67                 ode_method: str,
 68                 time_unit: str):
 69        self.p0 = p0
 70        self.t_min = t_min
 71        self.t_max = t_max
 72        self.processes = processes
 73        self.ode_method = ode_method
 74        self.time_unit = time_unit
 75
 76        # *** For computing and storing the solution to the system of ODEs. ***
 77        self._all_species, self._procs_by_reactant, self._procs_by_product = update_all_species(
 78            tuple(self.processes))
 79        self.odes = dict()  # store species-specific ODEs
 80        self.odes_sol = None  # numerical solution of species trajectories
 81
 82        self.species_with_ode = list()
 83        self.jac = None  # Jacobian matrix
 84        self.fixed_pts = dict()  # dictionary of fixed points
 85
 86        self.setup_ODEs()
 87        # self.get_fixed_pts()
 88
 89    def setup_ODEs(self, agent_based=True):
 90        """
 91        Set up the system of ODEs to be used for computing the
 92        deterministic trajectory of all the species in the given processes.
 93        The equations consist of sympy objects.
 94
 95        Parameters
 96        ----------
 97        agent_based : bool, default: True
 98            If True, set up the agent-based (or microscopic) form of ODEs.
 99            For instance, for the process `2X -> Y`, the ODE for species
100            `X` would include an `X(X - 1)` term instead
101            of `X^2` (the canonical form). If False, the canonical form
102            of the ODEs is constructed.
103
104        Notes
105        -----
106        The rate constant, `k`, for a given process is taken to be the mean
107        of `k`, unless `k` was defined to be normally-distributed,
108        in which case `k` is a 2-tuple and `k[0]` is the specified mean.
109
110        Implementation note: Building up the ODE expressions by separately iterating
111        over the processes for products and reactants. This is to properly
112        handle 0th order processes for product species. For example, for the
113        0th order process ' --> X' with rate constant k1, the ODE is dX/dt = k1.
114        """
115        # Set up ODEs for product species
116        for spec, procs in self._procs_by_product.items():
117            terms = []
118            for proc in procs:
119                k_val = proc.k[0] if isinstance(proc.k, tuple) else mean(proc.k)
120
121                if agent_based:
122                    temp = [Mul(*[sympify(name, locals=_clash) - i for i in
123                                  range(order)]) if name != '' else 1
124                            for name, order in proc.reactants.items()]
125                else:
126                    temp = [Pow(sympify(name, locals=_clash), order) if name != '' else 1
127                            for name, order in proc.reactants.items()]
128
129                new_term = Mul(proc.products[spec], k_val, *temp)
130                terms.append(new_term * self.get_term_multiplier(proc))
131
132            if spec in self.odes:
133                self.odes[spec] += Add(*terms)
134            else:
135                self.odes[spec] = Add(*terms)
136
137        # Set up ODEs for reactant species
138        for spec, procs in self._procs_by_reactant.items():
139            terms = []
140            for proc in procs:
141                k_val = proc.k[0] if isinstance(proc.k, tuple) else mean(proc.k)
142
143                if agent_based:
144                    temp = [Mul(*[sympify(name, locals=_clash) - i for i in range(order)]) for
145                            name, order in proc.reactants.items()]
146                else:
147                    temp = [Pow(sympify(name, locals=_clash), order) for name, order in
148                            proc.reactants.items()]
149
150                new_term = Mul(proc.reactants[spec], -1 * k_val, *temp)
151                terms.append(new_term * self.get_term_multiplier(proc))
152
153            if spec in self.odes:
154                self.odes[spec] += Add(*terms)
155            else:
156                self.odes[spec] = Add(*terms)
157
158        self.species_with_ode = [sp for sp in self.odes.keys()]
159
160        # Handle species whose population remains constant (e.g., a catalyst)
161        if len(self._all_species) != len(self.species_with_ode):
162            missing_species = self._all_species - set(self.species_with_ode)
163            for m in missing_species:
164                self.odes[m] = sympify(0)  # dm/dt = 0
165                self.species_with_ode.append(m)
166
167    @staticmethod
168    def get_term_multiplier(proc):
169        """
170        Generate the multiplicative term (or terms) needed for generating the
171        correct algebraic expressions for specialized processes
172        (such as Michaelis-Menten and regulated processes).
173        """
174        if isinstance(proc, MichaelisMentenProcess):
175            Km_val = proc.Km[0] if isinstance(proc.Km, tuple) else mean(proc.Km)
176            c = symbols(proc.catalyst)
177            return c / (Km_val + sympify(proc.reacts_[0], locals=_clash))
178        elif issubclass(type(proc), RegulatedProcess):
179            if isinstance(proc.regulating_species, list):
180                term = 1
181                for sp, a, Kval, hc in zip(proc.regulating_species, proc.alpha, proc.K50, proc.nH):
182                    K50_val = Kval[0] if isinstance(Kval, tuple) else mean(Kval)
183                    ratio = (symbols(sp) / K50_val) ** hc
184                    term *= (1 + a * ratio) / (1 + ratio)
185            else:
186                K50_val = proc.K50[0] if isinstance(proc.K50, tuple) else mean(proc.K50)
187                ratio = (symbols(proc.regulating_species) / K50_val) ** proc.nH
188                term = (1 + proc.alpha * ratio) / (1 + ratio)
189
190            if isinstance(proc, RegulatedMichaelisMentenProcess):
191                Km_val = proc.Km[0] if isinstance(proc.Km, tuple) else mean(proc.Km)
192                c = symbols(proc.catalyst)
193                term *= c / (Km_val + sympify(proc.reacts_[0], locals=_clash))
194
195            return term
196        else:
197            return 1
198
199    def solve_ODEs(self):
200        """
201        Solve system of ordinary differential equations (ODEs).
202
203        Notes
204        -----
205        Using the solver `scipy.integrate.solve_ivp`,
206        whose method can be one of the following:
207
208        - RK45 : Explicit Runge-Kutta method of order 5(4).
209        - RK23 : Explicit Runge-Kutta method of order 3(2).
210        - DOP853 : Explicit Runge-Kutta method of order 8.
211        - Radau : Implicit Runge-Kutta method of the Radau IIA family of order 5
212        - BDF : Implicit multistep variable-order (1 to 5) method based on a
213                 backward differentiation formula for the derivative approximation.
214        - LSODA : Adams/BDF method with automatic stiffness detection and switching.
215
216        Documentation
217        -------------
218        https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
219        """
220        # Converting ODE expressions: sympy -> scipy/numpy
221        _odes_lambdified = lambdify([sympify(sp, locals=_clash) for sp in self.odes.keys()],
222                                    list(self.odes.values()))
223
224        # Converting initial population values to a list to ensure that the order of
225        # species-specific ODEs and initial values are correctly ordered.
226        p0_list = [self.p0[sp] for sp in self.odes.keys()]
227        self.odes_sol = solve_ivp(lambda t, S: _odes_lambdified(*S),
228                                  t_span=[self.t_min, self.t_max],
229                                  y0=p0_list,
230                                  method=self.ode_method,
231                                  t_eval=None,  # Specify points where the solution is desired
232                                  dense_output=True)  # Compute a continuous solution
233
234    # def get_fixed_pts(self):
235    #     """
236    #     Not currently implemented.
237    #     """
238    #     pass

Perform deterministic calculations given the processes specified in an AbStochKin simulation object.

Attributes

p0 : dict[str (int]): Dictionary specifying the initial population sizes of all species in the given processes.
t_min (float or int): Numerical value of the start of simulated time in the units specified in the class attribute time_unit.
t_max (float or int): Numerical value of the end of simulated time in the units specified in the class attribute time_unit.
processes (list of Process objects): A list of all the processes that define the system of ODEs. Each process is a Process object.
ode_method (str): Method for the ODE solver to use.
time_unit (str): A string of the time unit to be used for describing the kinetics of the given processes.

DEcalcs( p0: dict, t_min: int | float, t_max: int | float, processes: list, ode_method: str, time_unit: str) View Source

62    def __init__(self,
63                 p0: dict,
64                 t_min: int | float,
65                 t_max: int | float,
66                 processes: list,
67                 ode_method: str,
68                 time_unit: str):
69        self.p0 = p0
70        self.t_min = t_min
71        self.t_max = t_max
72        self.processes = processes
73        self.ode_method = ode_method
74        self.time_unit = time_unit
75
76        # *** For computing and storing the solution to the system of ODEs. ***
77        self._all_species, self._procs_by_reactant, self._procs_by_product = update_all_species(
78            tuple(self.processes))
79        self.odes = dict()  # store species-specific ODEs
80        self.odes_sol = None  # numerical solution of species trajectories
81
82        self.species_with_ode = list()
83        self.jac = None  # Jacobian matrix
84        self.fixed_pts = dict()  # dictionary of fixed points
85
86        self.setup_ODEs()
87        # self.get_fixed_pts()

t_min

t_max

processes

ode_method

time_unit

odes

odes_sol

species_with_ode

jac

fixed_pts

def setup_ODEs(self, agent_based=True): View Source

 89    def setup_ODEs(self, agent_based=True):
 90        """
 91        Set up the system of ODEs to be used for computing the
 92        deterministic trajectory of all the species in the given processes.
 93        The equations consist of sympy objects.
 94
 95        Parameters
 96        ----------
 97        agent_based : bool, default: True
 98            If True, set up the agent-based (or microscopic) form of ODEs.
 99            For instance, for the process `2X -> Y`, the ODE for species
100            `X` would include an `X(X - 1)` term instead
101            of `X^2` (the canonical form). If False, the canonical form
102            of the ODEs is constructed.
103
104        Notes
105        -----
106        The rate constant, `k`, for a given process is taken to be the mean
107        of `k`, unless `k` was defined to be normally-distributed,
108        in which case `k` is a 2-tuple and `k[0]` is the specified mean.
109
110        Implementation note: Building up the ODE expressions by separately iterating
111        over the processes for products and reactants. This is to properly
112        handle 0th order processes for product species. For example, for the
113        0th order process ' --> X' with rate constant k1, the ODE is dX/dt = k1.
114        """
115        # Set up ODEs for product species
116        for spec, procs in self._procs_by_product.items():
117            terms = []
118            for proc in procs:
119                k_val = proc.k[0] if isinstance(proc.k, tuple) else mean(proc.k)
120
121                if agent_based:
122                    temp = [Mul(*[sympify(name, locals=_clash) - i for i in
123                                  range(order)]) if name != '' else 1
124                            for name, order in proc.reactants.items()]
125                else:
126                    temp = [Pow(sympify(name, locals=_clash), order) if name != '' else 1
127                            for name, order in proc.reactants.items()]
128
129                new_term = Mul(proc.products[spec], k_val, *temp)
130                terms.append(new_term * self.get_term_multiplier(proc))
131
132            if spec in self.odes:
133                self.odes[spec] += Add(*terms)
134            else:
135                self.odes[spec] = Add(*terms)
136
137        # Set up ODEs for reactant species
138        for spec, procs in self._procs_by_reactant.items():
139            terms = []
140            for proc in procs:
141                k_val = proc.k[0] if isinstance(proc.k, tuple) else mean(proc.k)
142
143                if agent_based:
144                    temp = [Mul(*[sympify(name, locals=_clash) - i for i in range(order)]) for
145                            name, order in proc.reactants.items()]
146                else:
147                    temp = [Pow(sympify(name, locals=_clash), order) for name, order in
148                            proc.reactants.items()]
149
150                new_term = Mul(proc.reactants[spec], -1 * k_val, *temp)
151                terms.append(new_term * self.get_term_multiplier(proc))
152
153            if spec in self.odes:
154                self.odes[spec] += Add(*terms)
155            else:
156                self.odes[spec] = Add(*terms)
157
158        self.species_with_ode = [sp for sp in self.odes.keys()]
159
160        # Handle species whose population remains constant (e.g., a catalyst)
161        if len(self._all_species) != len(self.species_with_ode):
162            missing_species = self._all_species - set(self.species_with_ode)
163            for m in missing_species:
164                self.odes[m] = sympify(0)  # dm/dt = 0
165                self.species_with_ode.append(m)

Set up the system of ODEs to be used for computing the deterministic trajectory of all the species in the given processes. The equations consist of sympy objects.

Parameters

agent_based : bool, default (True): If True, set up the agent-based (or microscopic) form of ODEs. For instance, for the process 2X -> Y, the ODE for species X would include an X(X - 1) term instead of X^2 (the canonical form). If False, the canonical form of the ODEs is constructed.

Notes

The rate constant, k, for a given process is taken to be the mean of k, unless k was defined to be normally-distributed, in which case k is a 2-tuple and k[0] is the specified mean.

Implementation note: Building up the ODE expressions by separately iterating over the processes for products and reactants. This is to properly handle 0th order processes for product species. For example, for the 0th order process ' --> X' with rate constant k1, the ODE is dX/dt = k1.

@staticmethod

def get_term_multiplier(proc): View Source

167    @staticmethod
168    def get_term_multiplier(proc):
169        """
170        Generate the multiplicative term (or terms) needed for generating the
171        correct algebraic expressions for specialized processes
172        (such as Michaelis-Menten and regulated processes).
173        """
174        if isinstance(proc, MichaelisMentenProcess):
175            Km_val = proc.Km[0] if isinstance(proc.Km, tuple) else mean(proc.Km)
176            c = symbols(proc.catalyst)
177            return c / (Km_val + sympify(proc.reacts_[0], locals=_clash))
178        elif issubclass(type(proc), RegulatedProcess):
179            if isinstance(proc.regulating_species, list):
180                term = 1
181                for sp, a, Kval, hc in zip(proc.regulating_species, proc.alpha, proc.K50, proc.nH):
182                    K50_val = Kval[0] if isinstance(Kval, tuple) else mean(Kval)
183                    ratio = (symbols(sp) / K50_val) ** hc
184                    term *= (1 + a * ratio) / (1 + ratio)
185            else:
186                K50_val = proc.K50[0] if isinstance(proc.K50, tuple) else mean(proc.K50)
187                ratio = (symbols(proc.regulating_species) / K50_val) ** proc.nH
188                term = (1 + proc.alpha * ratio) / (1 + ratio)
189
190            if isinstance(proc, RegulatedMichaelisMentenProcess):
191                Km_val = proc.Km[0] if isinstance(proc.Km, tuple) else mean(proc.Km)
192                c = symbols(proc.catalyst)
193                term *= c / (Km_val + sympify(proc.reacts_[0], locals=_clash))
194
195            return term
196        else:
197            return 1

Generate the multiplicative term (or terms) needed for generating the correct algebraic expressions for specialized processes (such as Michaelis-Menten and regulated processes).

def solve_ODEs(self): View Source

199    def solve_ODEs(self):
200        """
201        Solve system of ordinary differential equations (ODEs).
202
203        Notes
204        -----
205        Using the solver `scipy.integrate.solve_ivp`,
206        whose method can be one of the following:
207
208        - RK45 : Explicit Runge-Kutta method of order 5(4).
209        - RK23 : Explicit Runge-Kutta method of order 3(2).
210        - DOP853 : Explicit Runge-Kutta method of order 8.
211        - Radau : Implicit Runge-Kutta method of the Radau IIA family of order 5
212        - BDF : Implicit multistep variable-order (1 to 5) method based on a
213                 backward differentiation formula for the derivative approximation.
214        - LSODA : Adams/BDF method with automatic stiffness detection and switching.
215
216        Documentation
217        -------------
218        https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
219        """
220        # Converting ODE expressions: sympy -> scipy/numpy
221        _odes_lambdified = lambdify([sympify(sp, locals=_clash) for sp in self.odes.keys()],
222                                    list(self.odes.values()))
223
224        # Converting initial population values to a list to ensure that the order of
225        # species-specific ODEs and initial values are correctly ordered.
226        p0_list = [self.p0[sp] for sp in self.odes.keys()]
227        self.odes_sol = solve_ivp(lambda t, S: _odes_lambdified(*S),
228                                  t_span=[self.t_min, self.t_max],
229                                  y0=p0_list,
230                                  method=self.ode_method,
231                                  t_eval=None,  # Specify points where the solution is desired
232                                  dense_output=True)  # Compute a continuous solution

Solve system of ordinary differential equations (ODEs).

Notes

Using the solver scipy.integrate.solve_ivp, whose method can be one of the following:

RK45 : Explicit Runge-Kutta method of order 5(4).
RK23 : Explicit Runge-Kutta method of order 3(2).
DOP853 : Explicit Runge-Kutta method of order 8.
Radau : Implicit Runge-Kutta method of the Radau IIA family of order 5
BDF : Implicit multistep variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation.
LSODA : Adams/BDF method with automatic stiffness detection and switching.

Documentation

https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html