# Graphical Conditions ensuring Equality between Differential and Mean Stochastic Dynamics

Hugo Buscemi, François Fages

► To cite this version:

Hugo Buscemi, François Fages. Graphical Conditions ensuring Equality between Differential and Mean Stochastic Dynamics. CMSB 2024 - 22nd International Conference on Computational Methods in Systems Biology,, Sep 2024, Pisa (IT), Italy. hal-04623675

HAL Id: hal-04623675

https://inria.hal.science/hal-04623675

Submitted on 25 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Graphical Conditions ensuring Equality between Differential and Mean Stochastic Dynamics

Hugo Buscemi¹,² and François Fages²[0000-0001-5650-8266]

¹ ENS Paris-Saclay, Univ. Paris-Saclay, France
Hugo.Buscemi@ens-paris-saclay.fr

² Inria Saclay, Palaiseau, France
Francois.Fages@inria.fr

**Abstract.** Complex systems can be advantageously modeled by formal reaction systems (RS), a.k.a. chemical reaction networks in chemistry. Reaction-based models can indeed be interpreted in a hierarchy of semantics, depending on the question at hand, most notably by Ordinary Differential Equations (ODEs), Continuous Time Markov Chains (CTMCs), discrete Petri nets and asynchronous Boolean transition systems. The last three semantics can be easily related in the framework of abstract interpretation. The first two are classically related by Kurtz's limit theorem which states that if reactions are density-dependent families, then, as the volume goes to infinity, the mean reactant concentrations of the CTMC tends towards the solution of the ODE. In the more realistic context of bounded volumes, it is easy to show, by moment closure, that the restriction to reactions with at most one reactant ensures similarly that the mean of the CTMC trajectories is equal to the solution of the ODE at all time points. In this paper, we generalize that result in presence of polyreactant reactions, by introducing the Stoichiometric Influence and Modification Graph (SIMG) of an RS, and by showing that the equality between the two interpretations holds for the variables that belong to distinct SIMG ancestors of polyreactant reactions. We illustrate this approach with several examples. Evaluation on BioModels reveals that the condition for all variables is satisfied on models with no polymolecular reaction only. However, our theorem can be applied selectively to certain variables of the model to provide insights into their behaviour within more complex systems. Interestingly, we also show that the equality holds for a basic oscillatory RS implementing the sine and cosine functions of time.

**Keywords:** chemical reaction networks, Kurtz's theorem, continuous time markov chains, ordinary differential equations, oscillators, graphical conditions.

## 1 Introduction

Inside cells, numerous chemical reactions can occur. It is valuable to understand the behaviour and state of each chemical species involved in these reactions.
# Graphical Conditions ensuring Equality between Differential and Mean Stochastic Dynamics

Hugo Buscemi¹,² and François Fages²[0000-0001-5650-8266]

¹ ENS Paris-Saclay, Univ. Paris-Saclay, France
Hugo.Buscemi@ens-paris-saclay.fr

² Inria Saclay, Palaiseau, France
Francois.Fages@inria.fr

**Abstract.** Complex systems can be advantageously modeled by formal reaction systems (RS), a.k.a. chemical reaction networks in chemistry. Reaction-based models can indeed be interpreted in a hierarchy of semantics, depending on the question at hand, most notably by Ordinary Differential Equations (ODEs), Continuous Time Markov Chains (CTMCs), discrete Petri nets and asynchronous Boolean transition systems. The last three semantics can be easily related in the framework of abstract interpretation. The first two are classically related by Kurtz's limit theorem which states that if reactions are density-dependent families, then, as the volume goes to infinity, the mean reactant concentrations of the CTMC tends towards the solution of the ODE. In the more realistic context of bounded volumes, it is easy to show, by moment closure, that the restriction to reactions with at most one reactant ensures similarly that the mean of the CTMC trajectories is equal to the solution of the ODE at all time points. In this paper, we generalize that result in presence of polyreactant reactions, by introducing the Stoichiometric Influence and Modification Graph (SIMG) of an RS, and by showing that the equality between the two interpretations holds for the variables that belong to distinct SIMG ancestors of polyreactant reactions. We illustrate this approach with several examples. Evaluation on BioModels reveals that the condition for all variables is satisfied on models with no polymolecular reaction only. However, our theorem can be applied selectively to certain variables of the model to provide insights into their behaviour within more complex systems. Interestingly, we also show that the equality holds for a basic oscillatory RS implementing the sine and cosine functions of time.

**Keywords:** chemical reaction networks, Kurtz's theorem, continuous time markov chains, ordinary differential equations, oscillators, graphical conditions.

## 1 Introduction

Inside cells, numerous chemical reactions can occur. It is valuable to understand the behaviour and state of each chemical species involved in these reactions.
Such reactions can be represented by a Chemical Reaction Network (CRN), which describes the species involved, their resulting outcomes, and the specific rates at which they react. An Ordinary Differential Equation (ODE) can be associated with such a CRN [7]. Given a fixed initial state, each ODE has a unique deterministic solution.

On the other hand, a CRN can also be interpreted as a Continuous Time Markov Chain (CTMC) and simulated by the Stochastic Simulation Algorithm (SSA), a.k.a. Gillespie's Algorithm [8]. Simulations starting from the same initial state exhibit a variety of non-deterministic trajectories and a mean trajectory. This stochastic mean behaviour may or may not align with the deterministic solution of the ODE starting from the same initial state. Kurtz's limit theorem [9] asserts that, under certain circumstances, the mean behaviour converges to the deterministic one when the volume tends to the infinity. These circumstances involve rates that are density-dependent, possibly apart from a limited number of mode-switching transient states [11].

In the actual context of a cell however, the volume is bounded, and there is also a finite, and possibly small, number of molecules for some species, e.g. as an extreme, only one for genes. It is in this real-world model that we aim to study the relationship between the two continuous-time interpretations of CRNs, by ODE and by CTMC. Regarding the exact correspondence between the mean stochastic CTMC behaviours and the deterministic ODE behaviour, it is well known that if all reactions in the CRN have at most one reactant, then the equality holds. This can be easily by the calculus of moments derived from the Chemical Master Equation[1].

Our contribution is an extension of the unimolecular case, while still ensuring equality. Our result allows for polymolecular reactions with uncorrelated species, by guaranteeing the absence of correlations through graph-theoretic conditions on the structure of the network.

For the sake of generality of these results for dynamical systems beyond the molecular interaction scale, outside of chemistry, we prefer to speak of formal reaction systems (RS) over mathematical variables [3], rather than CRNs over molecular species. Apart from this update of terminology, there is no difference between RSs and CRNs, which are equipped with the same hierarchy of differential, stochastic, Petri net and Boolean semantics [6].

The rest of the paper is organized as follows. The next section illustrates the discrepancy there is between the differential and stochastic interpretations of a simple RS modeling the classical Lotka-Volterra prey-predator population dynamics system. The following section provides the necessary preliminaries on reaction systems with mass action kinetics and their differential and mean stochastic semantics, based respectively on ODE and CTMC interpretations. In Sec. 4, we introduce the Stoichiometric Influence and Modification Graph (SIMG) which somehow combines the classical notions of reaction hypergraph and influence graph of a reaction system [5]. That notion of SIMG is used in Sec. 5 to express and prove our main theorem of strong equality between the differential ODE and mean stochastic CTMC interpretations of RSs satisfying
Such reactions can be represented by a Chemical Reaction Network (CRN), which describes the species involved, their resulting outcomes, and the specific rates at which they react. An Ordinary Differential Equation (ODE) can be associated with such a CRN [7]. Given a fixed initial state, each ODE has a unique deterministic solution.

On the other hand, a CRN can also be interpreted as a Continuous Time Markov Chain (CTMC) and simulated by the Stochastic Simulation Algorithm (SSA), a.k.a. Gillespie's Algorithm [8]. Simulations starting from the same initial state exhibit a variety of non-deterministic trajectories and a mean trajectory. This stochastic mean behaviour may or may not align with the deterministic solution of the ODE starting from the same initial state. Kurtz's limit theorem [9] asserts that, under certain circumstances, the mean behaviour converges to the deterministic one when the volume tends to the infinity. These circumstances involve rates that are density-dependent, possibly apart from a limited number of mode-switching transient states [11].

In the actual context of a cell however, the volume is bounded, and there is also a finite, and possibly small, number of molecules for some species, e.g. as an extreme, only one for genes. It is in this real-world model that we aim to study the relationship between the two continuous-time interpretations of CRNs, by ODE and by CTMC. Regarding the exact correspondence between the mean stochastic CTMC behaviours and the deterministic ODE behaviour, it is well known that if all reactions in the CRN have at most one reactant, then the equality holds. This can be easily by the calculus of moments derived from the Chemical Master Equation[1].

Our contribution is an extension of the unimolecular case, while still ensuring equality. Our result allows for polymolecular reactions with uncorrelated species, by guaranteeing the absence of correlations through graph-theoretic conditions on the structure of the network.

For the sake of generality of these results for dynamical systems beyond the molecular interaction scale, outside of chemistry, we prefer to speak of formal reaction systems (RS) over mathematical variables [3], rather than CRNs over molecular species. Apart from this update of terminology, there is no difference between RSs and CRNs, which are equipped with the same hierarchy of differential, stochastic, Petri net and Boolean semantics [6].

The rest of the paper is organized as follows. The next section illustrates the discrepancy there is between the differential and stochastic interpretations of a simple RS modeling the classical Lotka-Volterra prey-predator population dynamics system. The following section provides the necessary preliminaries on reaction systems with mass action kinetics and their differential and mean stochastic semantics, based respectively on ODE and CTMC interpretations. In Sec. 4, we introduce the Stoichiometric Influence and Modification Graph (SIMG) which somehow combines the classical notions of reaction hypergraph and influence graph of a reaction system [5]. That notion of SIMG is used in Sec. 5 to express and prove our main theorem of strong equality between the differential ODE and mean stochastic CTMC interpretations of RSs satisfying
some structural conditions, independently of the kinetics and any limit consid-
erations. In Sec. 6, we provide a low complexity algorithm to test our SIMG
conditions. Then in the following sections, we provide several examples to il-
lustrate these results. Interestingly, we also show that strong equality property
for an oscillatory RS which implements the sine and cosine functions of time in
both differential and stochastic interpretations. Then we report on the evalua-
tion of our condition in the repository of RS models Biomodels [10]. Finally, we
conclude on some perspectives for future work with three open questions.

## 2 Motivating Example

<figure><img src="image_1.png" /><figcaption>Fig. 1: Deterministic ODE simulation of Lotka-Volterra reaction system.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2: Stochastic CTMC simulations of Lotka-Volterra reaction system.</figcaption></figure>
some structural conditions, independently of the kinetics and any limit consid-
erations. In Sec. 6, we provide a low complexity algorithm to test our SIMG
conditions. Then in the following sections, we provide several examples to il-
lustrate these results. Interestingly, we also show that strong equality property
for an oscillatory RS which implements the sine and cosine functions of time in
both differential and stochastic interpretations. Then we report on the evalua-
tion of our condition in the repository of RS models Biomodels [10]. Finally, we
conclude on some perspectives for future work with three open questions.

## 2 Motivating Example

<figure><img src="image_2.png" /><figcaption>Fig. 1: Deterministic ODE simulation of Lotka-Volterra reaction system.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Stochastic CTMC simulations of Lotka-Volterra reaction system.</figcaption></figure>
The Lotka-Volterra reaction system, $F + R \to 2.F$, $R \to 2.R$, $F \to \emptyset$, is a classical model of interaction between two populations of preys (rabbits $R$) and predators (foxes $F$). When interpreted using ODEs, the system exhibits oscillatory behaviour, with the populations of rabbits and foxes fluctuating over time (see Fig. 1). However, when interpreted by a CTMC, the system leads almost surely to the extinction of the predator or both populations (see Fig. 2a). This is due to the presence in the stochastic semantics of an absorbing state 0 that is neither reached in the continuous semantics, nor by Kurtz's limit mean CTMC trajectory obtained by making the volume tends to infinity [9]. This is one reason why Kurtz's theorem should not be applied to this kind of counter-example.

On the other hand, one can remark that the addition of immigration and emigration reactions prevents the system from reaching an absorbing state where one population is completely extinct. By allowing for the possibility of individuals entering and leaving the system, there is no absorbing state anymore and the populations exhibit sustained stochastic oscillations (see Fig. 2b).

### 3 Preliminaries on Reaction Systems

We recall here the standard notions of Chemical Reaction Networks (CRN) over molecular species [7], by using a more general terminology of Reaction Systems (RS) over mathematical variables, in order to apply them to population dynamics or dynamical processes that have nothing to do with chemistry. Apart from this change of terminology to reflect the generality of that formalism of reaction systems, there is no difference between both formalisms of RSs and CRNs.

**Definition 1.** A reaction $\rho$ over a set of variables $X$ is a triple $(R, P, \Lambda)$, also noted

$$ \rho : \sum_{i} \ell_{i} R_{i} \Rightarrow^{A} \sum_{j} r_{j} P_{j} $$

where $R, P$ are multisets of variables, $R_i, P_j \in X$ are the elements of $R, P$ with non-zero coefficient multiplicity $\ell_i, r_j \in \mathbb{N}^*$, and $\Lambda : X \to \mathbb{R}^+$ is a rate function. A Reaction System (RS) over $X$ is a finite set of reactions over $X$.

The $R_i$'s are called reactants and the $P_i$'s products. A reaction with no reactant is called a synthesis reaction. A reaction with only one reactant and coefficient equal to 1 is called a unireactant reaction. Other reactions are polyreactant reactions. A polyreactant reaction with no product is called an annihilation reaction.

**Notation 1** For a reaction $\rho$ and a variable $x \in X$, we write $\zeta_{\rho,x} = r_{\rho,x} - \ell_{\rho,x}$ the difference of coefficients between the production of product $x$ and the consumption of reactant $x$ in $\rho$.

In this paper, we shall consider mass action law rate functions only.
The Lotka-Volterra reaction system, $F + R \to 2.F$, $R \to 2.R$, $F \to \emptyset$, is a classical model of interaction between two populations of preys (rabbits $R$) and predators (foxes $F$). When interpreted using ODEs, the system exhibits oscillatory behaviour, with the populations of rabbits and foxes fluctuating over time (see Fig. 1). However, when interpreted by a CTMC, the system leads almost surely to the extinction of the predator or both populations (see Fig. 2a). This is due to the presence in the stochastic semantics of an absorbing state 0 that is neither reached in the continuous semantics, nor by Kurtz's limit mean CTMC trajectory obtained by making the volume tends to infinity [9]. This is one reason why Kurtz's theorem should not be applied to this kind of counter-example.

On the other hand, one can remark that the addition of immigration and emigration reactions prevents the system from reaching an absorbing state where one population is completely extinct. By allowing for the possibility of individuals entering and leaving the system, there is no absorbing state anymore and the populations exhibit sustained stochastic oscillations (see Fig. 2b).

### 3 Preliminaries on Reaction Systems

We recall here the standard notions of Chemical Reaction Networks (CRN) over molecular species [7], by using a more general terminology of Reaction Systems (RS) over mathematical variables, in order to apply them to population dynamics or dynamical processes that have nothing to do with chemistry. Apart from this change of terminology to reflect the generality of that formalism of reaction systems, there is no difference between both formalisms of RSs and CRNs.

**Definition 1.** A reaction $\rho$ over a set of variables $X$ is a triple $(R, P, \Lambda)$, also noted

$$ \rho : \sum_{i} \ell_{i} R_{i} \Rightarrow^{\Lambda} \sum_{j} r_{j} P_{j} $$

where $R, P$ are multisets of variables, $R_i, P_j \in X$ are the elements of $R, P$ with non-zero coefficient multiplicity $\ell_i, r_j \in \mathbb{N}^*$, and $\Lambda : X \to \mathbb{R}^+$ is a rate function. A Reaction System (RS) over $X$ is a finite set of reactions over $X$.

The $R_i$'s are called reactants and the $P_i$'s products. A reaction with no reactant is called a synthesis reaction. A reaction with only one reactant and coefficient equal to 1 is called a unireactant reaction. Other reactions are polyreactant reactions. A polyreactant reaction with no product is called an annihilation reaction.

**Notation 1** For a reaction $\rho$ and a variable $x \in X$, we write $\zeta_{\rho,x} = r_{\rho,x} - \ell_{\rho,x}$ the difference of coefficients between the production of product $x$ and the consumption of reactant $x$ in $\rho$.

In this paper, we shall consider mass action law rate functions only.
**Definition 2 (Mass-Action Law).** The mass action law rate of a reaction $(R, P, \Lambda)$ is of the form $\Lambda = \lambda \prod_i \prod_{k=0}^{l_i-1} (R_i - k)$ for some rate constant $\lambda \in \mathbb{R}$.

*Example 1.* Mass action RS are written with just their rate constant, which can also be omitted if equal to 1, as in the reactions below:

$$ 
\begin{aligned}
A &\Rightarrow^2 A + B \\
- &\Rightarrow B \\
B + C &\Rightarrow^6 B + D \\
D &\Rightarrow -
\end{aligned}
 $$

**Definition 3 (Associated ODE).** The Ordinary Differential Equation (ODE) of an RS associates to each variable $x \in X$ the differential equation

$$ \frac{dx}{dt} = \sum_{\rho} \Lambda_{\rho}(X) \zeta_{\rho, x} $$

We denote by $\sigma$ the unique solution of the system starting at the initial state $\sigma(0)$, and $\sigma_x$ its projection on the variable $x$.

*Example 2.* The ODE associated with the RS of Ex. 1 is

$$ 
\begin{aligned}
\frac{dA}{dt} &= 0 & \frac{dB}{dt} &= 2A + 1 \\
\frac{dC}{dt} &= -6BC & \frac{dD}{dt} &= 6BC - D
\end{aligned}
 $$

**Definition 4 (Associated Continuous Time Markov Chain).** The Continuous Time Markov Chain (CTMC) associated to an RS is defined by the following transition probability

$$ \mathbb{P}\{X(t + \Delta t) = j | X(t) = i\} \sim_{\Delta t \to 0} q_{ij} \Delta t $$

where

$$ q_{X(t),X(t)+\zeta} = \sum_{\rho: r_{\rho} - l_{\rho} = \zeta} \Lambda_{\rho}(X(t)) $$

We shall make use of the following properties shown in [9] of the CTMC associated to an RS.

**Proposition 1.** The mean of $X(t)$ at each point of time follows the equation

$$ \frac{d\mathbb{E}[X(s)]}{dt} = \sum_{\rho} \mathbb{E}[\Lambda_{\rho}(X(s))] \zeta_{\rho} $$

**Definition 5 (Generator).** The generator of the CTMC is the function $A$ defined for each function $f$ by

$$ Af(x) = \sum_{\rho} \Lambda_{\rho}(x)(f(x + \zeta_{\rho}) - f(x)) $$
**Definition 2 (Mass-Action Law).** The mass action law rate of a reaction $(R, P, \Lambda)$ is of the form $\Lambda = \lambda \prod_i \prod_{k=0}^{l_i-1} (R_i - k)$ for some rate constant $\lambda \in \mathbb{R}$.

*Example 1.* Mass action RS are written with just their rate constant, which can also be omitted if equal to 1, as in the reactions below:

$$ 
\begin{aligned}
A &\Rightarrow^2 A + B \\
- &\Rightarrow B \\
B + C &\Rightarrow^6 B + D \\
D &\Rightarrow -
\end{aligned}
 $$

**Definition 3 (Associated ODE).** The Ordinary Differential Equation (ODE) of an RS associates to each variable $x \in X$ the differential equation

$$ \frac{dx}{dt} = \sum_{\rho} \Lambda_{\rho}(X) \zeta_{\rho, x} $$

We denote by $\sigma$ the unique solution of the system starting at the initial state $\sigma(0)$, and $\sigma_x$ its projection on the variable $x$

*Example 2.* The ODE associated with the RS of Ex. 1 is

$$ 
\begin{aligned}
\frac{dA}{dt} &= 0 & \frac{dB}{dt} &= 2A + 1 \\
\frac{dC}{dt} &= -6BC & \frac{dD}{dt} &= 6BC - D
\end{aligned}
 $$

**Definition 4 (Associated Continuous Time Markov Chain).** The Continuous Time Markov Chain (CTMC) associated to an RS is defined by the following transition probability

$$ \mathbb{P}\{X(t + \Delta t) = j | X(t) = i\} \sim_{\Delta t \to 0} q_{ij} \Delta t $$

where

$$ q_{X(t),X(t)+\zeta} = \sum_{\rho: r_{\rho} - l_{\rho} = \zeta} \Lambda_{\rho}(X(t)) $$

We shall make use of the following properties shown in [9] of the CTMC associated to an RS.

**Proposition 1.** The mean of $X(t)$ at each point of time follows the equation

$$ \frac{d\mathbb{E}[X(s)]}{dt} = \sum_{\rho} \mathbb{E}[\Lambda_{\rho}(X(s))] \zeta_{\rho} $$

**Definition 5 (Generator).** The generator of the CTMC is the function $A$ defined for each function $f$ by

$$ Af(x) = \sum_{\rho} \Lambda_{\rho}(x)(f(x + \zeta_{\rho}) - f(x)) $$
**Proposition 2.** $f(X(t)) - f(X(0)) - \int_0^t Af(X(s))ds$ is a martingale and so
$$ \mathbb{E}[f(X(t))] = \mathbb{E}[f(X(0))] + \int_0^t \mathbb{E}[Af(X(s))] ds. \quad \frac{d\mathbb{E}[f(X(s))]}{dt} = \mathbb{E}[Af(X(t))] $$

When $\zeta_{\rho,x} = 0$ the reaction $\rho$ has no direct influence on the variable $x$, and one can focus on reactions where $\zeta_{\rho,x} \neq 0$ only.

# 4 Stoichiometric Influence and Modification Graph

In this section, we introduce a special kind of influence graph associated to a reaction system, which will be used to define our structural conditions ensuring the equality between the mean stochastic and differential interpretations of a given RS.

**Definition 6 (SIMG).** The Stoichiometric Influence and Modification Graph (SIMG) of an RS is a bipartite multigraph with variables and reactions as vertices, with $l_{\rho,x}$ arcs from variable $x$ to reaction $\rho$, and one additional arc from $\rho$ to $x$ if $\zeta_{\rho,x} \neq 0$.

*Example 3.* The SIMG of the RS of Ex. 1 is the following:

<figure><img src="image_3.png" /></figure>

Once this directed graph is constructed, the directed acyclic graph (DAG) of its strongly connected components (SCCs) can be created. This DAG provides a hierarchy and depth for each node (reactions and variables) within the graph. In this graph, each node's ancestors are all the SCCs (including itself) higher in the hierarchy and leading to it.

**Definition 7 (Ancestors).** A variable $y \in X$ is called an ancestor of a variable $x \in X$ if there exists a path in the SIMG from $y$ to $x$. Especially, $x$ is its own ancestor.

**Definition 8 (Depth).** The depth of a SCC is defined recursively as follows:

- All SCCs with no arc leading to them are of depth 0.

- Any other SCC is at depth $n+1$ where $n$ is the maximum depth of the SCCs leading to it.
**Proposition 2.** $f(X(t)) - f(X(0)) - \int_0^t Af(X(s))ds$ is a martingale and so
$$ \mathbb{E}[f(X(t))] = \mathbb{E}[f(X(0))] + \int_0^t \mathbb{E}[Af(X(s))] ds. \quad \frac{d\mathbb{E}[f(X(s))]}{dt} = \mathbb{E}[Af(X(t))] $$

When $\zeta_{\rho,x} = 0$ the reaction $\rho$ has no direct influence on the variable $x$, and one can focus on reactions where $\zeta_{\rho,x} \neq 0$ only.

# 4 Stoichiometric Influence and Modification Graph

In this section, we introduce a special kind of influence graph associated to a reaction system, which will be used to define our structural conditions ensuring the equality between the mean stochastic and differential interpretations of a given RS.

**Definition 6 (SIMG).** The Stoichiometric Influence and Modification Graph (SIMG) of an RS is a bipartite multigraph with variables and reactions as vertices, with $\ell_{\rho,x}$ arcs from variable $x$ to reaction $\rho$, and one additional arc from $\rho$ to $x$ if $\zeta_{\rho,x} \neq 0$.

*Example 3.* The SIMG of the RS of Ex. 1 is the following:

<figure><img src="image_4.png" /></figure>

Once this directed graph is constructed, the directed acyclic graph (DAG) of its strongly connected components (SCCs) can be created. This DAG provides a hierarchy and depth for each node (reactions and variables) within the graph. In this graph, each node's ancestors are all the SCCs (including itself) higher in the hierarchy and leading to it.

**Definition 7 (Ancestors).** A variable $y \in X$ is called an ancestor of a variable $x \in X$ if there exists a path in the SIMG from $y$ to $x$. Especially, $x$ is its own ancestor.

**Definition 8 (Depth).** The depth of a SCC is defined recursively as follows:

- All SCCs with no arc leading to them are of depth 0.

- Any other SCC is at depth $n+1$ where $n$ is the maximum depth of the SCCs leading to it.
By abuse of notation, we call the depth of a vertex y, the depth of the SCC of y.

Example 4. We show in the next section that the SIMG of the following RS satisfies graphical conditions that entail the equality between the deterministic ODE solution and the mean CTMC behaviour. Simulation traces are shown in Fig. 3a. The actual difference between the mean CTMC trace and the deterministic ODE solution is shown in Fig. 3b, using a scale that shows that the difference is indeed neglectable.

$$A+B \Rightarrow A+B+C$$

$$A+B \Rightarrow A+B+D$$

$$C \Rightarrow D$$

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>(a) Stochastic mean trace and deterministic solution</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Difference between mean stochastic and deterministic solutions</figcaption></figure>

<caption>Fig. 3: Simulation traces for Ex. 4 based on 10000 stochastic traces.</caption>

# 5 Main Theorem

**Theorem 1.** For an RS with the mass-action law, if each reaction is

1. either at most unireactant,

2. or polyreactant and has all reactants with disjoint SIMG ancestors,

then the mean stochastic behaviour, starting from a fixed initial state, is equal to the differential solution with the same initial state.
By abuse of notation, we call the depth of a vertex y, the depth of the SCC of y.

Example 4. We show in the next section that the SIMG of the following RS satisfies graphical conditions that entail the equality between the deterministic ODE solution and the mean CTMC behaviour. Simulation traces are shown in Fig. 3a. The actual difference between the mean CTMC trace and the deterministic ODE solution is shown in Fig. 3b, using a scale that shows that the difference is indeed neglectable.

$$A + B \Rightarrow A + B + C$$

$$A + B \Rightarrow A + B + D$$

$$C \Rightarrow D$$

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>(a) Stochastic mean trace and deterministic solution</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Difference between mean stochastic and deterministic solutions</figcaption></figure>

<caption>Fig. 3: Simulation traces for Ex. 4 based on 10000 stochastic traces.</caption>

# 5 Main Theorem

**Theorem 1.** For an RS with the mass-action law, if each reaction is

1. either at most unireactant,

2. or polyreactant and has all reactants with disjoint SIMG ancestors,

then the mean stochastic behaviour, starting from a fixed initial state, is equal to the differential solution with the same initial state.
*Proof.* Let us examine the different types of reactions that can occur in the RS. As a reminder, they have the form:

$$ \rho : \sum_{i} \ell_{\rho, R_i} R_i \Rightarrow^{\Lambda_{\rho}} \sum_{j} r_{\rho, P_j} P_j $$

A first observation is that each variable is its own ancestor, which means it cannot have multiple arcs leading to the same reaction. Hence, for any reaction $\rho$ and variable $x$, we have $\ell_{\rho,x} \in \{0, 1\}$.

Therefore, for a specific variable $x$, reactions with $\zeta_{\rho,x} \neq 0$ can be categorized into the following five sets according to the reactants and their ancestors:

$$ 1. \pi_0(x) = \{ \rho : \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} > 0 \} $$

$$ 2. \pi_1(x) = \{ \rho : x \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} \neq 1 \} $$

$$ 3. \pi'_1(x) = \{ \rho : R \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} > 0, R \neq x \} $$

$$ 4. \pi_2(x) = \{ \rho : x + \sum_i R_i \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} \neq 1, R_i \neq x \} \text{ and } x \text{ and all } R_i \text{ have disjoint ancestors.} $$

$$ 5. \pi'_2(x) = \{ \rho : \sum_i R_i \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} > 0, R_i \neq x \} \text{ and } R_i \text{ have disjoint ancestors} $$

It is worth noting that if $\rho \in \pi_2(x)$, then reactant $R_0$ is an ancestor of both $x$ and $R_0$ itself. However, this violates the disjoint ancestors property condition. Therefore, we conclude that $\pi_2(x) = \emptyset$.

Now, we can express the set of differential equations for the mean value of each variable:

$$ \frac{d\mathbb{E}[x]}{dt} = \sum_{\rho \in \pi_0(x)} \lambda_{\rho} \zeta_{\rho,x} + \sum_{\rho \in \pi_1(x)} \lambda_{\rho} \mathbb{E}[x] \zeta_{\rho,x} + \sum_{\rho \in \pi'_1(x)} \lambda_{\rho} \mathbb{E}[R] \zeta_{\rho,x} + \sum_{\rho \in \pi'_2(x)} \lambda_{\rho} \mathbb{E}\left[\prod_i R_i\right] \zeta_{\rho,x} $$

Now, by applying Lemma 2 proven below, the differential equation becomes:

$$ \frac{d\mathbb{E}[x]}{dt} = \sum_{\rho \in \pi_0(x)} \lambda_{\rho} \zeta_{\rho,x} + \sum_{\rho \in \pi_1(x)} \lambda_{\rho} \mathbb{E}[x] \zeta_{\rho,x} + \sum_{\rho \in \pi'_1(x)} \lambda_{\rho} \mathbb{E}[R] \zeta_{\rho,x} + \sum_{\rho \in \pi'_2(x)} \lambda_{\rho} \prod_i \mathbb{E}[R_i] \zeta_{\rho,x} $$
*Proof.* Let us examine the different types of reactions that can occur in the RS. As a reminder, they have the form:

$$ \rho : \sum_{i} \ell_{\rho, R_i} R_i \Rightarrow^{\Lambda_{\rho}} \sum_{j} r_{\rho, P_j} P_j $$

A first observation is that each variable is its own ancestor, which means it cannot have multiple arcs leading to the same reaction. Hence, for any reaction $\rho$ and variable $x$, we have $\ell_{\rho,x} \in \{0, 1\}$.

Therefore, for a specific variable $x$, reactions with $\zeta_{\rho,x} \neq 0$ can be categorized into the following five sets according to the reactants and their ancestors:

$$ 1. \pi_0(x) = \{ \rho : \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} > 0 \} $$

$$ 2. \pi_1(x) = \{ \rho : x \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} \neq 1 \} $$

$$ 3. \pi'_1(x) = \{ \rho : R \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} > 0, R \neq x \} $$

$$ 4. \pi_2(x) = \{ \rho : x + \sum_i R_i \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} \neq 1, R_i \neq x \} \text{ and } x \text{ and all } R_i \text{ have disjoint ancestors.} $$

$$ 5. \pi'_2(x) = \{ \rho : \sum_i R_i \Rightarrow^{\Lambda_{\rho}} r_{\rho,x}x + \sum_j r_{\rho,P_j} P_j \mid r_{\rho,x} > 0, R_i \neq x \} \text{ and } R_i \text{ have disjoint ancestors} $$

It is worth noting that if $\rho \in \pi_2(x)$, then reactant $R_0$ is an ancestor of both $x$ and $R_0$ itself. However, this violates the disjoint ancestors property condition. Therefore, we conclude that $\pi_2(x) = \emptyset$.

Now, we can express the set of differential equations for the mean value of each variable:

$$ \frac{d\mathbb{E}[x]}{dt} = \sum_{\rho \in \pi_0(x)} \lambda_{\rho} \zeta_{\rho,x} + \sum_{\rho \in \pi_1(x)} \lambda_{\rho} \mathbb{E}[x] \zeta_{\rho,x} + \sum_{\rho \in \pi'_1(x)} \lambda_{\rho} \mathbb{E}[R] \zeta_{\rho,x} + \sum_{\rho \in \pi'_2(x)} \lambda_{\rho} \mathbb{E}\left[\prod_i R_i\right] \zeta_{\rho,x} $$

Now, by applying Lemma 2 proven below, the differential equation becomes:

$$ \frac{d\mathbb{E}[x]}{dt} = \sum_{\rho \in \pi_0(x)} \lambda_{\rho} \zeta_{\rho,x} + \sum_{\rho \in \pi_1(x)} \lambda_{\rho} \mathbb{E}[x] \zeta_{\rho,x} + \sum_{\rho \in \pi'_1(x)} \lambda_{\rho} \mathbb{E}[R] \zeta_{\rho,x} + \sum_{\rho \in \pi'_2(x)} \lambda_{\rho} \prod_i \mathbb{E}[R_i] \zeta_{\rho,x} $$
This is exactly the same set of differential equations satisfied by the differential solution $\sigma$. Furthermore, we have $\mathbb{E}[X](0) = \sigma(0)$. Therefore, we conclude that $\mathbb{E}[X] = \sigma$.

**Lemma 1.** Under the hypotheses of the theorem, a reaction is either alone in its SCC, or all the reactions of its SCC are unireatant reactions.

*Proof.* By simple inspection of the different cases.

**Lemma 2.** Under the hypotheses of the theorem, for all variables $x_i$'s with no common ancestors, the following equality holds:

$$ \mathbb{E} \left[ \prod_i x_i \right] = \prod_i \mathbb{E} [x_i] $$

*Proof.* Since the variables $x_i$'s have no common ancestors, they are not modified by a same reaction. Therefore, for any reaction $\rho$, we have $\zeta_{\rho,x_i}\zeta_{\rho,x_j} = 0$ for any $x_i$ and $x_j$ with no common ancestor and no $x_i$ and $x_j$ are in the same SCC. We prove the result by induction on the maximum depth of the variables occurring in the product.

*Base case:* At depth 0, for any $i$, the sets $\pi_0(x_i)$ and $\pi_2'(x_i)$ are empty otherwise the SIMG would have one of this patterns $\rho \to x_i$ or $R_1, R_2 \Rightarrow \rho \to x_i$ with $\rho$ in its own SCC, using Lemma 1, and so $x_i$ would have a depth > 0.

By using the equation from Proposition 2, we have:

$$ \frac{d\mathbb{E}[\prod_{i \in I} x_i]}{dt} = \sum_{j \in I} \sum_{\rho \in \pi_1(x_j)} \lambda_{\rho} \mathbb{E} \left[ \prod_{i \in I} x_i \right] \zeta_{\rho, x_j} + \sum_{\rho \in \pi_1'(x_j)} \lambda_{\rho} \mathbb{E} \left[ R \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} $$

For the product of means, we have:

$$ \frac{d \prod_{i \in I} \mathbb{E}[x_i]}{dt} = \sum_{j \in I} \sum_{\rho \in \pi_1(x_j)} \lambda_{\rho} \prod_{i \in I} \mathbb{E}[x_i] \zeta_{\rho, x_j} + \sum_{\rho \in \pi_1'(x_j)} \lambda_{\rho} \mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i] \zeta_{\rho, x_j} $$

All variables appearing in those equations are of depth at most 0. We can rewrite those equation for any variable in each SCC of the $x_i$'s. This gives us a whole system of equations with one equation for each possible product of variables. Since only those products appear in the system, it has a unique solution at a fixed initial state. The system of equations for the mean of the product and the product of means have the same structure. Since they have the same initial conditions, they are equivalent. Therefore, $\mathbb{E}[\prod_i x_i] = \prod_i \mathbb{E}[x_i]$ holds for any product of variables of depth 0.

*Inductive step:* Let us assume $\Delta$ that the result holds for any product of variables with a maximum depth of $n$, and prove it for products of variables with
This is exactly the same set of differential equations satisfied by the differential solution $\sigma$. Furthermore, we have $\mathbb{E}[X](0) = \sigma(0)$. Therefore, we conclude that $\mathbb{E}[X] = \sigma$.

**Lemma 1.** Under the hypotheses of the theorem, a reaction is either alone in its SCC, or all the reactions of its SCC are unireatant reactions.

*Proof.* By simple inspection of the different cases.

**Lemma 2.** Under the hypotheses of the theorem, for all variables $x_i$'s with no common ancestors, the following equality holds:

$$ \mathbb{E}\left[\prod_{i} x_{i}\right] = \prod_{i} \mathbb{E}[x_{i}] $$

*Proof.* Since the variables $x_i$'s have no common ancestors, they are not modified by a same reaction. Therefore, for any reaction $\rho$, we have $\zeta_{\rho,x_i}\zeta_{\rho,x_j} = 0$ for any $x_i$ and $x_j$ with no common ancestor and no $x_i$ and $x_j$ are in the same SCC. We prove the result by induction on the maximum depth of the variables occurring in the product.

*Base case:* At depth 0, for any $i$, the sets $\pi_0(x_i)$ and $\pi_2'(x_i)$ are empty otherwise the SIMG would have one of this patterns $\rho \to x_i$ or $R_1, R_2 \Rightarrow \rho \to x_i$ with $\rho$ in its own SCC, using Lemma 1, and so $x_i$ would have a depth > 0.

By using the equation from Proposition 2, we have:

$$ \frac{d\mathbb{E}[\prod_{i \in I} x_i]}{dt} = \sum_{j \in I} \sum_{\rho \in \pi_1(x_j)} \lambda_{\rho} \mathbb{E}\left[\prod_{i \in I} x_i\right] \zeta_{\rho, x_j} + \sum_{\rho \in \pi_1'(x_j)} \lambda_{\rho} \mathbb{E}\left[R \prod_{i \neq j} x_i\right] \zeta_{\rho, x_j} $$

For the product of means, we have:

$$ \frac{d \prod_{i \in I} \mathbb{E}[x_i]}{dt} = \sum_{j \in I} \sum_{\rho \in \pi_1(x_j)} \lambda_{\rho} \prod_{i \in I} \mathbb{E}[x_i] \zeta_{\rho, x_j} + \sum_{\rho \in \pi_1'(x_j)} \lambda_{\rho} \mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i] \zeta_{\rho, x_j} $$

All variables appearing in those equations are of depth at most 0. We can rewrite those equation for any variable in each SCC of the $x_i$'s. This gives us a whole system of equations with one equation for each possible product of variables. Since only those products appear in the system, it has a unique solution at a fixed initial state. The system of equations for the mean of the product and the product of means have the same structure. Since they have the same initial conditions, they are equivalent. Therefore, $\mathbb{E}[\prod_i x_i] = \prod_i \mathbb{E}[x_i]$ holds for any product of variables of depth 0.

*Inductive step:* Let us assume $\Delta$ that the result holds for any product of variables with a maximum depth of $n$, and prove it for products of variables with
a maximum depth of $n + 1$. We use induction on the $n + 2$-tuple $(d_{n+1}, \dots, d_0)$ where $d_i$ represents the number of variables of the product at depth $i$, with lexicographic order.

★ Base case: If there are no variable at depth $n + 1$, then the induction hypothesis $\Delta$ gives the result for products of variables with a maximum depth of $n$. Thus, the result holds.

★ Induction step: Assume that the result holds for any product with at most $k$ variables of depth exactly $n + 1$ (and all others at most $n$). Considering a product with exactly $k + 1$ variables of depth exactly $n + 1$. We denote its $n + 2$-tuple by $\tau = (k + 1, d_n, \dots, d_0)$. Consider the differential equation system for the product $\prod_{i \in I} x_i$:

$$ \frac{d\mathbb{E}[\prod_{i \in I} x_i]}{dt} = \sum_{j \in I} \sum_{\rho \in \pi_0(x_j)} \lambda_\rho \mathbb{E} \left[ \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} \\ + \sum_{\rho \in \pi_1(x_j)} \lambda_\rho \mathbb{E} \left[ \prod_{i \in I} x_i \right] \zeta_{\rho, x_j} \\ + \sum_{\rho \in \pi'_1(x_j)} \lambda_\rho \mathbb{E} \left[ R \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} \\ + \sum_{\rho \in \pi'_2(x_j)} \lambda_\rho \mathbb{E} \left[ \prod_{\ell} R_\ell \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} $$

For $\rho \in \pi_0(x_j)$, the tuple $\tau_{j,\rho} = \tau[d(x_j) \leftarrow \tau(d(x_j)) - 1]$ is such that $\tau_{j,\rho} < \tau$. Then, the induction hypothesis applies to the product $\mathbb{E}[\prod_{i \neq j} x_i] = \prod_{i \neq j} \mathbb{E}[x_i]$.

Since for $\rho \in \pi'_2(x_j)$, using Lemma 1, the $R_\ell$ are higher in the hierarchy than $x_j$, and they have no common ancestors with $x_i$ for $i \neq j$. So we have $\tau_{j,\rho} < \tau$ and thus $\mathbb{E}[\prod_\ell R_\ell \prod_{i \neq j} x_i] = \prod_\ell \mathbb{E}[R_\ell] \prod_{i \neq j} \mathbb{E}[x_i]$. The same resoning applies to $\rho \in \pi'_1(x_j)$ such that $R$ is not from the SCC of $x_j$. In this case, $R$ is higher and independant from the other variables, and $\mathbb{E}[R \prod_{i \neq j} x_i] = \mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i]$.

It remains to consider $\rho \in \pi'_1(x_j)$ such that $R$ is part of the SCC of $x_j$, and also $\rho \in \pi_1(x_j)$ which is the same but with $R = x_j$. In this case, we can write the system of differential equations for $R \prod_{i \neq j} x_i$ for all $R$ of the SCC. Since we have the same number of equations as unknown functions and the initial values are fixed, we conclude that $\mathbb{E}[R \prod_{i \neq j} x_i](0) = (\mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i])(0)$, which implies $\mathbb{E}[R \prod_{i \neq j} x_i] = \mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i]$ for the considered product.
a maximum depth of $n + 1$. We use induction on the $n + 2$-tuple $(d_{n+1}, \dots, d_0)$ where $d_i$ represents the number of variables of the product at depth $i$, with lexicographic order.

★ Base case: If there are no variable at depth $n + 1$, then the induction hypothesis $\Delta$ gives the result for products of variables with a maximum depth of $n$. Thus, the result holds.

★ Induction step: Assume that the result holds for any product with at most $k$ variables of depth exactly $n + 1$ (and all others at most $n$). Considering a product with exactly $k + 1$ variables of depth exactly $n + 1$. We denote its $n + 2$-tuple by $\tau = (k + 1, d_n, \dots, d_0)$. Consider the differential equation system for the product $\prod_{i \in I} x_i$:

$$ \frac{d\mathbb{E}[\prod_{i \in I} x_i]}{dt} = \sum_{j \in I} \sum_{\rho \in \pi_0(x_j)} \lambda_\rho \mathbb{E} \left[ \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} \\ + \sum_{\rho \in \pi_1(x_j)} \lambda_\rho \mathbb{E} \left[ \prod_{i \in I} x_i \right] \zeta_{\rho, x_j} \\ + \sum_{\rho \in \pi'_1(x_j)} \lambda_\rho \mathbb{E} \left[ R \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} \\ + \sum_{\rho \in \pi'_2(x_j)} \lambda_\rho \mathbb{E} \left[ \prod_{\ell} R_\ell \prod_{i \neq j} x_i \right] \zeta_{\rho, x_j} $$

For $\rho \in \pi_0(x_j)$, the tuple $\tau_{j,\rho} = \tau[d(x_j) \leftarrow \tau(d(x_j)) - 1]$ is such that $\tau_{j,\rho} < \tau$. Then, the induction hypothesis applies to the product $\mathbb{E}[\prod_{i \neq j} x_i] = \prod_{i \neq j} \mathbb{E}[x_i]$.

Since for $\rho \in \pi'_2(x_j)$, using Lemma 1, the $R_\ell$ are higher in the hierarchy than $x_j$, and they have no common ancestors with $x_i$ for $i \neq j$. So we have $\tau_{j,\rho} < \tau$ and thus $\mathbb{E}[\prod_\ell R_\ell \prod_{i \neq j} x_i] = \prod_\ell \mathbb{E}[R_\ell] \prod_{i \neq j} \mathbb{E}[x_i]$. The same resoning applies to $\rho \in \pi'_1(x_j)$ such that $R$ is not from the SCC of $x_j$. In this case, $R$ is higher and independant from the other variables, and $\mathbb{E}[R \prod_{i \neq j} x_i] = \mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i]$.

It remains to consider $\rho \in \pi'_1(x_j)$ such that $R$ is part of the SCC of $x_j$, and also $\rho \in \pi_1(x_j)$ which is the same but with $R = x_j$. In this case, we can write the system of differential equations for $R \prod_{i \neq j} x_i$ for all $R$ of the SCC. Since we have the same number of equations as unknown functions and the initial values are fixed, we conclude that $\mathbb{E}[R \prod_{i \neq j} x_i](0) = (\mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i])(0)$, which implies $\mathbb{E}[R \prod_{i \neq j} x_i] = \mathbb{E}[R] \prod_{i \neq j} \mathbb{E}[x_i]$ for the considered product.
## 6 SIMG Ancestor Condition Checking Algorithm

**input**: SIMG of an RS of size *n* (number of variables and reactions)
**output**: Is the graphical condition of Thm. 1 verified?

1 `sccs()` calculus; // in O(n)
2 Sort sccs() with depth(scc 1) ≤ depth(scc 2) ≤ ...;
3 `foreach` scc ∈ sccs do
4 | Check scc is unireactant; // in O(|scc|)
5 | ancestors ← scc;
6 | `foreach` variable x ∈ scc do
7 | | `foreach` reaction ρ with x as SIMG successor do
8 | | | match reactants of ρ with :
9 | | | | ∅ : do nothing;
10 | | | | {R} : ancestors.add(R.ancestors);
11 | | | | {R₁...Rₖ} :
12 | | | | if ∀i < j, Rᵢ.ancestors ∩ Rⱼ.ancestors = ∅ then
13 | | | | ancestors.add(∪<sub>i=1</sub><sup>k</sup> Rᵢ.ancestors)
14 | | | | else
15 | | | | return False
16 | | | | end
17 | end
18 | end
19 | `foreach` variable x ∈ scc do
20 | | x.ancestors ← ancestors;
21 | | end
22 | end
23 | return True

<caption>**Algorithm 1:** Checking of the graphical condition</caption>

The construction of the SIMG of an RS can be easily computed from the reaction graph and is of linear size. Now, Alg. 1 provides a means to verify whether the graphical condition is fulfilled given the SIMG of an RS. By simple inspection of the loops, one can show:

**Proposition 3.** (*termination and complexity*) Alg. 1 terminates in $O(n^2)$ steps where *n* is the number of variables and reactions and *d* is the maximal number of reactants for a reaction.

**Proposition 4.** (*correctness*) Alg. 1 checks the condition of Thm. 1.

*Proof.* Let us prove its correctness by first showing that for each variable *x*, *x.ancestors* are indeed the real ancestors of *x*. Then we'll see why this implies the correctness of the algorithm.

First, it is clear that all members of a SCC have the same ancestors.
## 6 SIMG Ancestor Condition Checking Algorithm

```c
input : SIMG of an RS of size n (number of variables and reactions)
output: Is the graphical condition of Thm. 1 verified?

1 sccs() calculus ; // in O(n)
2 Sort sccs() with depth(scc 1) ≤ depth(scc 2) ≤ ...;
3 foreach scc ∈ sccs do
4 | Check scc is unireactant ; // in O(|scc|)
5 | ancestors ← scc;
6 | foreach variable x ∈ scc do
7 | | foreach reaction ρ with x as SIMG successor do
8 | | | match reactants of ρ with :
9 | | | | ∅ : do nothing;
10 | | | | {R} : ancestors.add(R.ancestors);
11 | | | | {R₁...Rₖ} :
12 | | | | if ∀i < j, Rᵢ.ancestors ∩ Rⱼ.ancestors = ∅ then
13 | | | | ancestors.add(∪<k>ᵢ=₁ Rᵢ.ancestors)
14 | | | | else
15 | | | | return False
16 | | end
17 | end
18 | end
19 | foreach variable x ∈ scc do
20 | | x.ancestors ← ancestors;
21 | end
22 end
23 | return True
```

<caption>**Algorithm 1:** Checking of the graphical condition</caption>

The construction of the SIMG of an RS can be easily computed from the reaction graph and is of linear size. Now, Alg. 1 provides a means to verify whether the graphical condition is fulfilled given the SIMG of an RS. By simple inspection of the loops, one can show:

**Proposition 3.** (*termination and complexity*) *Alg. 1 terminates in O(n<sup>2</sup>) steps where n is the number of variables and reactions and d is the maximal number of reactants for a reaction.*

**Proposition 4.** (*correctness*) *Alg. 1 checks the condition of Thm. 1.*

*Proof.* Let us prove its correctness by first showing that for each variable *x*, *x.ancestors* are indeed the real ancestors of *x*. Then we'll see why this implies the correctness of the algorithm.

First, it is clear that all members of a SCC have the same ancestors.
Now, let us prove by induction on the depth of SCCs that $x.ancestors$ is actually the set of ancestors of $x$.

- At depth 0:

* $\supseteq$ The only ancestors of any $x$ at depth 0 are the variables from its SCC. They are added in Alg. line 5

* $\subseteq$ The only reactions leading to any $x$ (in the SIMG) are the ones from its SCC. Then only Alg. line 10 adds ancestors but it is nothing else than the SCC itself.

So ancestors at depth 0 are well computed.

- At depth $n$:

* $\supseteq$ If $A$ is an ancestor of $x$

* Either $A$ is in the same SCC than $x$ and is added in Alg. line 5

* Either it is not and there exists a path in the SIMG such that $A \to \dots \to R \to \Box \to x'$, with some $R$ such that $depth(R) < n$ and $x'$ is in the same SCC than $x$. Then by induction hypothesis, $A.ancestors$ is well formed. So $R \in A.ancestors$ and $A.ancestors$ is adding to $x'.ancestors = x.ancestors$ when considering $x'$ in Alg. line 6.

Then all ancestors of $x$ are in $x.ancestors$.

* $\subseteq$ If $A$ is in $x.ancestors$.

* Either it is added with Alg. line 13 and was in some $R_i.ancestors$ such that $R_i \to \Box \to x'$ with $x'$ in the same SCC than $x$. So $depth(R_i) < n$ and then by induction hypothesis $A$ is an ancestor of $R_i$ and so of $x'$ and of $x$.

* Either it is added with Alg. line 10 and so $R \to \Box \to x'$ with $x'$ in the SCC of $x$ and $A$ in $R.ancestors$. The only way to be adde is then with $R.ancestors \neq \emptyset$ which can't happen with $R$ in the same SCC. So $depth(R) < n$ and the same reasoning as just before with induction hypothesis and transitivity of ancestors leads to $A$ as an ancestor of $x$.

So all variables in $x.ancestors$ are true ancestors of $x$.

Therefore $x.ancestors$ is well computed for each $x$ at depth $n$.

By the induction principle, it is the case for any depth. Finally, checking the condition of Thm. 1 is done by Alg. line 4 and line 12.

# 7 Examples

*Example 5.* One can check that the following RS satisfies the condition of Thm. 1 on all variables. Fig. 4 shows the equality between the ODE trace and the mean CTMC trace, together with their high variance.
Now, let us prove by induction on the depth of SCCs that $x.ancestors$ is actually the set of ancestors of $x$.

- At depth 0:

* $\supseteq$ The only ancestors of any $x$ at depth 0 are the variables from its SCC. They are added in Alg. line 5

* $\subseteq$ The only reactions leading to any $x$ (in the SIMG) are the ones from its SCC. Then only Alg. line 10 adds ancestors but it is nothing else than the SCC itself.

So ancestors at depth 0 are well computed.

- At depth $n$:

* $\supseteq$ If $A$ is an ancestor of $x$

* Either $A$ is in the same SCC than $x$ and is added in Alg. line 5

* Either it is not and there exists a path in the SIMG such that $A \to \dots \to R \to \square \to x'$, with some $R$ such that $depth(R) < n$ and $x'$ is in the same SCC than $x$. Then by induction hypothesis, $A.ancestors$ is well formed. So $R \in A.ancestors$ and $A.ancestors$ is adding to $x'.ancestors = x.ancestors$ when considering $x'$ in Alg. line 6.

Then all ancestors of $x$ are in $x.ancestors$.

* $\subseteq$ If $A$ is in $x.ancestors$.

* Either it is added with Alg. line 13 and was in some $R_i.ancestors$ such that $R_i \to \square \to x'$ with $x'$ in the same SCC than $x$. So $depth(R_i) < n$ and then by induction hypothesis $A$ is an ancestor of $R_i$ and so of $x'$ and of $x$.

* Either it is added with Alg. line 10 and so $R \to \square \to x'$ with $x'$ in the SCC of $x$ and $A$ in $R.ancestors$. The only way to be adde is then with $R.ancestors \neq \emptyset$ which can't happen with $R$ in the same SCC. So $depth(R) < n$ and the same reasoning as just before with induction hypothesis and transitivity of ancestors leads to $A$ as an ancestor of $x$.

So all variables in $x.ancestors$ are true ancestors of $x$.

Therefore $x.ancestors$ is well computed for each $x$ at depth $n$.

By the induction principle, it is the case for any depth. Finally, checking the condition of Thm. 1 is done by Alg. line 4 and line 12.

# 7 Examples

*Example 5.* One can check that the following RS satisfies the condition of Thm. 1 on all variables. Fig. 4 shows the equality between the ODE trace and the mean CTMC trace, together with their high variance.
$$A+B \Rightarrow A+B+C$$

$$C \Rightarrow D$$

$$D \Rightarrow \_$$

$$D \Rightarrow D+E$$

$$E \Rightarrow E+F+G$$

$$E \Rightarrow \_$$

$$F \Rightarrow \_$$

$$F \Rightarrow F+C$$

$$G \Rightarrow H$$

$$E+I \Rightarrow E+I+H$$

$$H \Rightarrow \_$$

$$D+J \Rightarrow D+J+K$$

$$K \Rightarrow \_$$

$$K+L \Rightarrow K+L+M$$

$$\_ \Rightarrow M$$

$$M \Rightarrow \_$$

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>(a) Stochastic mean and deterministic solutions.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Difference between mean stochastic and deterministic solutions.</figcaption></figure>

<caption>Fig. 4: ODE trace and 100 CTMC traces with mean in Ex. 5.</caption>

*Example 6.* The following RS implements arithmetic multiplication [4]. Indeed, with A and B constant, in the ODE semantics we get at stable state, $\mathbb{E}[C] = AB$. Thm. 1 applies on all variables and thus entails the correctness of the result in the mean stochastic interpretation as well:
$$A+B \Rightarrow A+B+C$$

$$C \Rightarrow D$$

$$D \Rightarrow \_$$

$$D \Rightarrow D+E$$

$$E \Rightarrow E+F+G$$

$$E \Rightarrow \_$$

$$F \Rightarrow \_$$

$$F \Rightarrow F+C$$

$$G \Rightarrow H$$

$$E+I \Rightarrow E+I+H$$

$$H \Rightarrow \_$$

$$D+J \Rightarrow D+J+K$$

$$K \Rightarrow \_$$

$$K+L \Rightarrow K+L+M$$

$$\_ \Rightarrow M$$

$$M \Rightarrow \_$$

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>(a) Stochastic mean and deterministic solutions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Difference between mean stochastic and deterministic solutions.</figcaption></figure>

<caption>Fig. 4: ODE trace and 100 CTMC traces with mean in Ex. 5.</caption>

*Example 6.* The following RS implements arithmetic multiplication [4]. Indeed, with $A$ and $B$ constant, in the ODE semantics we get at stable state, $\mathbb{E}[C] = AB$. Thm. 1 applies on all variables and thus entails the correctness of the result in the mean stochastic interpretation as well:
$$A + B \Rightarrow A + B + C \\ C \Rightarrow -$$

<figure><img src="image_10.png" /></figure>

*Example 7.* Let us consider the following RS with SIMG violating the conditions:

$$A + B \Rightarrow A + B + C \\ C \Rightarrow - \\ C \Rightarrow C + D \\ D \Rightarrow - \\ C + D \Rightarrow C + D + E \\ E \Rightarrow -$$

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 5: Difference between the mean over 100 iterations and the deterministic solution for Ex. 7 with $A(0) = 5$ and $B(0) = 11$</figcaption></figure>

The difference of behaviours between the stochastic and deterministic interpretations is plotted in Fig. 5 for the different variables. Interestingly, one can observe no difference on all variables except E. This can be proven on the SIMG as follows. Because the sub-graph until D verifies the graphical condition, we have $\mathbb{E}[A] = \sigma_A = a$, $\mathbb{E}[B] = \sigma_B = b$ and $\mathbb{E}[C] = \sigma_C$ and $\mathbb{E}[D] = \sigma_D$.
$$A + B \Rightarrow A + B + C$$
$$C \Rightarrow \_$$

<figure><img src="image_11.png" /></figure>

*Example 7.* Let us consider the following RS with SIMG violating the conditions:

$$A + B \Rightarrow A + B + C$$
$$C \Rightarrow \_$$
$$C \Rightarrow C + D$$
$$D \Rightarrow \_$$
$$C + D \Rightarrow C + D + E$$
$$E \Rightarrow \_$$

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 5: Difference between the mean over 100 iterations and the deterministic solution for Ex. 7 with $A(0) = 5$ and $B(0) = 11$</figcaption></figure>

The difference of behaviours between the stochastic and deterministic interpretations is plotted in Fig. 5 for the different variables. Interestingly, one can observe no difference on all variables except E. This can be proven on the SIMG as follows. Because the sub-graph until D verifies the graphical condition, we have $E[A] = \sigma_A = a$, $E[B] = \sigma_B = b$ and $E[C] = \sigma_C$ and $E[D] = \sigma_D$.
Then $\frac{d\sigma_C}{dt} = ab + \sigma_C$ so, with $\sigma_C(0) = 0$ we have $\sigma_C = ab(1 - e^{-t})$. Thus

$$ 
\begin{aligned}
\frac{d\xi}{dt}(t) &= \frac{d(\mathbb{E}[C^2] - \sigma_C^2)}{dt}(t) = \mathbb{E}[ab((C+1)^2 - C^2) + C((c-1)^2 - C^2)](t) \\
&\quad - (2(ab - \sigma_C)\sigma_C)(t) \\
&= (ab + \sigma_C - 2\xi)(t)
\end{aligned}
$$

and $\xi(t) = ab(1 - e^{-t})$.

The same reasoning leads for $\psi = \mathbb{E}[CD] - \sigma_C\sigma_D$ to $\psi' + 2\psi = \xi$ and so $\psi(t) = \frac{ab}{2}(1 - 2e^{-t} + e^{-2t})$.

Then $\phi = \mathbb{E}[E] - \sigma_E$ verifies $\phi' + \phi = \psi$ and so $\phi = \frac{ab}{2}(1 - 2te^{-t} - e^{-2t}) \neq 0$.

## 8 Equality Property for Cosine Oscillatory RS

The sine and cosine functions of time constitute a fundamental example of oscillator. A CRN implementing $sin(t)$ and $cos(t)$ with 4 molecular species was given in [4] for illustrating the Turing-completeness of continuous CRNs, and the compiler of computable mathematical functions in CRNs implemented in BIOCHAM<fnref n="1" /> [2]. That automatically generated RS, called cosRS, uses the dual-rail encoding of negative values as the difference between the concentrations of two reactants with a fast annihilation reaction between both with rate $f$:

$$ 
\begin{aligned}
A &\Rightarrow A + B \\
B &\Rightarrow B + C \\
C &\Rightarrow C + D \\
D &\Rightarrow D + A \\
A + C &\Rightarrow^f \_ \\
B + D &\Rightarrow^f \_ \\
\end{aligned}
$$

With initial conditions $A(0) = 1$ and $B(0) = C(0) = D(0) = 0$, we have $A - C = \cos(t)$ and $D - B = \sin(t)$ in the ODE semantics of cosRS. We show that the same holds for $\mathbb{E}[A - C]$ and $\mathbb{E}[D - B]$ in the stochastic CTMC semantics, and does generalize to any initial conditions.

**Theorem 2.** For any initial conditions, the mean stochastic traces of $A - C$ and $D - B$ are equal to their ODE traces in the ODE semantics of cosRS.

*Proof.* We just have to show that the mean stochastic CTMC traces of $D - B$ and $A - C$ obey the differential equations of the ODE semantics of cosRS.

Variables $A$ and $C$ represent the positive and negative parts of the same function, $A - C$, while $B$ and $D$ represent similarly another function, $D - B$.

By linearity, we have $\frac{d\mathbb{E}[A - C]}{dt} = \mathbb{E}[D - B]$ and $\frac{d\mathbb{E}[D - B]}{dt} = \mathbb{E}[C - A]$.

Therefore, we get $\frac{d^2\mathbb{E}[A - C]}{dt^2} = -\mathbb{E}[A - C]$ and $\frac{d^2\mathbb{E}[D - B]}{dt^2} = -\mathbb{E}[D - B]$.

<footnote n="1">¹ http://lifeware.inria.fr/biocham</footnote>
Then $\frac{d\sigma_C}{dt} = ab + \sigma_C$ so, with $\sigma_C(0) = 0$ we have $\sigma_C = ab(1 - e^{-t})$. Thus

$$ 
\begin{aligned}
\frac{d\xi}{dt}(t) &= \frac{d(\mathbb{E}[C^2] - \sigma_C^2)}{dt}(t) = \mathbb{E}[ab((C+1)^2 - C^2) + C((c-1)^2 - C^2)](t) \\
&\quad - (2(ab - \sigma_C)\sigma_C)(t) \\
&= (ab + \sigma_C - 2\xi)(t)
\end{aligned}
$$

and $\xi(t) = ab(1 - e^{-t})$.

The same reasoning leads for $\psi = \mathbb{E}[CD] - \sigma_C\sigma_D$ to $\psi' + 2\psi = \xi$ and so $\psi(t) = \frac{ab}{2}(1 - 2e^{-t} + e^{-2t})$.

Then $\phi = \mathbb{E}[E] - \sigma_E$ verifies $\phi' + \phi = \psi$ and so $\phi = \frac{ab}{2}(1 - 2te^{-t} - e^{-2t}) \neq 0$.

## 8 Equality Property for Cosine Oscillatory RS

The sine and cosine functions of time constitute a fundamental example of oscillator. A CRN implementing $sin(t)$ and $cos(t)$ with 4 molecular species was given in [4] for illustrating the Turing-completeness of continuous CRNs, and the compiler of computable mathematical functions in CRNs implemented in BIOCHAM<fnref n="1" /> [2]. That automatically generated RS, called cosRS, uses the dual-rail encoding of negative values as the difference between the concentrations of two reactants with a fast annihilation reaction between both with rate $f$:

$$ 
\begin{aligned}
A &\Rightarrow A + B \\
B &\Rightarrow B + C \\
C &\Rightarrow C + D \\
D &\Rightarrow D + A \\
A + C &\Rightarrow^f \_ \\
B + D &\Rightarrow^f \_ \\
\end{aligned}
$$

With initial conditions $A(0) = 1$ and $B(0) = C(0) = D(0) = 0$, we have $A - C = \cos(t)$ and $D - B = \sin(t)$ in the ODE semantics of cosRS. We show that the same holds for $\mathbb{E}[A - C]$ and $\mathbb{E}[D - B]$ in the stochastic CTMC semantics, and does generalize to any initial conditions.

**Theorem 2.** For any initial conditions, the mean stochastic traces of $A - C$ and $D - B$ are equal to their ODE traces in the ODE semantics of cosRS.

*Proof.* We just have to show that the mean stochastic CTMC traces of $D - B$ and $A - C$ obey the differential equations of the ODE semantics of cosRS.

Variables $A$ and $C$ represent the positive and negative parts of the same function, $A - C$, while $B$ and $D$ represent similarly another function, $D - B$.

By linearity, we have $\frac{d\mathbb{E}[A - C]}{dt} = \mathbb{E}[D - B]$ and $\frac{d\mathbb{E}[D - B]}{dt} = \mathbb{E}[C - A]$.

Therefore, we get $\frac{d^2\mathbb{E}[A - C]}{dt^2} = -\mathbb{E}[A - C]$ and $\frac{d^2\mathbb{E}[D - B]}{dt^2} = -\mathbb{E}[D - B]$.

<footnote n="1">¹ http://lifeware.inria.fr/biocham</footnote>
# 9 Evaluation in Biomodels and Partial Approximations

Alg. 1 has been applied to the first 2000 models of the database BioModels [10]. The only models satisfying the condition did not contain polyomolecular reactions. As mentioned in the introduction, the property for at most unimolecular reaction models is well-known and our graphical condition for all variables does not apply to more models in Biomodels repository. Nevertheless, even when the conditions of Thm. 1 are not satisfied on the whole SIMG of the model, they can be used to identify a subset of reactant variables for which the ODE solution is indeed equal to the mean CTMC behaviour.

The key observation is that the condition is additive with respect to the variables ordering based on their depth in the SIMG. If variables do not violate the condition, then it is true for all previous variables considered in the algorithm. Therefore, the algorithm only needs to focus on the partial portion of the graph where all ancestors satisfy the condition. The result then holds for those variables and means nothing for the other variables.

A partial application of Thm. 1 can thus be used to identify RS input variables which may be non constant, but with provably equal differential and stochastic mean dynamics. Independently of the more complex behaviours of the other variables, the theorem guarantees that the ODE semantics provides the exact mean concentration of these input variables over time.

This is illustrated here in Ex. 7 where the partial condition applies to variables A, B, C and D. This explains why in Fig. 5, variables A, B, C and D follow the deterministic ODE solution in mean, unlike variable E. In Biomodels, the statistics are however not very significant since many reactions have no mass-action law kinetics, or have non well-formed kinetics [3].

# 10 Conclusion and Perspectives

Kurtz's theorem is the state-of-the-art result to relate the mean stochastic CTMC and deterministic ODE interpretations of a reaction system. As a limit theorem however, it does not apply to bounded systems and leads to wrong results if applied for instance to population dynamics models in finite space.

We have shown that the strong equality holds between the mean CTMC and deterministic ODE interpretations independently of the kinetics and without any limit approximation, provided some graphical conditions are satisfied on the structure the RS. The interest of our result is that it is not a limit theorem and can thus be applied to bounded dynamical systems with low values for the variables, such as models of gene regulatory networks if we think that one gene is in single copy in a cell.

The evaluation of our graphical conditions in the reaction models in Biomodels shows, on the one hand, that they are satisfied for all species, only in models containing no polymolecular reactions, but, on the other hand, that they can be used to identify in each model a subset of molecular species for which the property holds. Interestingly, we have also shown that the basic oscillatory RS
# 9 Evaluation in Biomodels and Partial Approximations

Alg. 1 has been applied to the first 2000 models of the database BioModels [10]. The only models satisfying the condition did not contain polyomolecular reactions. As mentioned in the introduction, the property for at most unimolecular reaction models is well-known and our graphical condition for all variables does not apply to more models in Biomodels repository. Nevertheless, even when the conditions of Thm. 1 are not satisfied on the whole SIMG of the model, they can be used to identify a subset of reactant variables for which the ODE solution is indeed equal to the mean CTMC behaviour.

The key observation is that the condition is additive with respect to the variables ordering based on their depth in the SIMG. If variables do not violate the condition, then it is true for all previous variables considered in the algorithm. Therefore, the algorithm only needs to focus on the partial portion of the graph where all ancestors satisfy the condition. The result then holds for those variables and means nothing for the other variables.

A partial application of Thm. 1 can thus be used to identify RS input variables which may be non constant, but with provably equal differential and stochastic mean dynamics. Independently of the more complex behaviours of the other variables, the theorem guarantees that the ODE semantics provides the exact mean concentration of these input variables over time.

This is illustrated here in Ex. 7 where the partial condition applies to variables A, B, C and D. This explains why in Fig. 5, variables A, B, C and D follow the deterministic ODE solution in mean, unlike variable E. In Biomodels, the statistics are however not very significant since many reactions have no mass-action law kinetics, or have non well-formed kinetics [3].

# 10 Conclusion and Perspectives

Kurtz's theorem is the state-of-the-art result to relate the mean stochastic CTMC and deterministic ODE interpretations of a reaction system. As a limit theorem however, it does not apply to bounded systems and leads to wrong results if applied for instance to population dynamics models in finite space.

We have shown that the strong equality holds between the mean CTMC and deterministic ODE interpretations independently of the kinetics and without any limit approximation, provided some graphical conditions are satisfied on the structure the RS. The interest of our result is that it is not a limit theorem and can thus be applied to bounded dynamical systems with low values for the variables, such as models of gene regulatory networks if we think that one gene is in single copy in a cell.

The evaluation of our graphical conditions in the reaction models in Biomodels shows, on the one hand, that they are satisfied for all species, only in models containing no polymolecular reactions, but, on the other hand, that they can be used to identify in each model a subset of molecular species for which the property holds. Interestingly, we have also shown that the basic oscillatory RS
implementing the sine and cosine functions enjoys the strong equality property for all species.

These results also leave some new questions open. First, back to our first example, can we determine the mean stochastic behaviour of the Lotka-Volterra system extended with immigration and emigration reactions for which we have shown the existence of sustained stochastic oscillations ? Second, can we generalize the reasoning done for the sine-cosine RS, to prove the property for general RSs similarly involving annihilation reactions as sole degradation reactions ? Third, our assumptions require no correlations between reactants, while some variables may be balanced and yield the same result, they do not meet our condition. These questions should deserve further investigations.

## Acknowledgements

We are grateful to the reviewers of CMSB for their comments, which where useful to improve the presentation of our results, and to Mathieu Hemery for discussions, and for improving the parsing of SBML models in BIOCHAM.

## References

1. Anderson, D.F., Kurtz, T.G.: Stochastic analysis of biochemical systems, vol. 674. Springer (2015)

2. Calzone, L., Fages, F., Soliman, S.: BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge. Bioinformatics 22(14), 1805-1807 (2006). https://doi.org/10.1093/bioinformatics/btl172

3. Fages, F., Gay, S., Soliman, S.: Inferring reaction systems from ordinary differential equations. Theoretical Computer Science 599, 64-78 (Sep 2015). https://doi.org/10.1016/j.tcs.2014.07.032

4. Fages, F., Le Guludec, G., Bournez, O., Pouly, A.: Strong Turing Completeness of Continuous Chemical Reaction Networks and Compilation of Mixed Analog-Digital Programs. In: CMSB'17: Proceedings of the fiveteen international conference on Computational Methods in Systems Biology. LNCS, vol. 10545, pp. 108-127. Springer-Verlag (Sep 2017). https://doi.org/10.1007/978-3-319-67471-1_7

5. Fages, F., Martinez, T., Rosenblueth, D., Soliman, S.: Influence networks compared with reaction networks: Semantics, expressivity and attractors. IEEE/ACM Transactions on Computational Biology and Bioinformatics (Feb 2018). https://doi.org/10.1109/TCBB.2018.2805686

6. Fages, F., Soliman, S.: Abstract interpretation and types for systems biology. Theoretical Computer Science 403(1), 52-70 (2008). https://doi.org/10.1016/j.tcs.2008.04.024

7. Feinberg, M.: Mathematical aspects of mass action kinetics. In: Lapidus, L., Amundson, N.R. (eds.) Chemical Reactor Theory: A Review, chap. 1, pp. 1-78. Prentice-Hall (1977)

8. Gillespie, D.T.: Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry 81(25), 2340-2361 (1977)

9. Kurtz, T.G.: Solutions of ordinary differential equations as limits of pure jump markov processes. Journal of applied Probability 7(1), 49-58 (1970)
implementing the sine and cosine functions enjoys the strong equality property for all species.

These results also leave some new questions open. First, back to our first example, can we determine the mean stochastic behaviour of the Lotka-Volterra system extended with immigration and emigration reactions for which we have shown the existence of sustained stochastic oscillations ? Second, can we generalize the reasoning done for the sine-cosine RS, to prove the property for general RSs similarly involving annihilation reactions as sole degradation reactions ? Third, our assumptions require no correlations between reactants, while some variables may be balanced and yield the same result, they do not meet our condition. These questions should deserve further investigations.

## Acknowledgements

We are grateful to the reviewers of CMSB for their comments, which where useful to improve the presentation of our results, and to Mathieu Hemery for discussions, and for improving the parsing of SBML models in BIOCHAM.

## References

1. Anderson, D.F., Kurtz, T.G.: Stochastic analysis of biochemical systems, vol. 674. Springer (2015)

2. Calzone, L., Fages, F., Soliman, S.: BIOCHAM: An environment for modeling biological systems and formalizing experimental knowledge. Bioinformatics **22**(14), 1805-1807 (2006). https://doi.org/10.1093/bioinformatics/btl172

3. Fages, F., Gay, S., Soliman, S.: Inferring reaction systems from ordinary differential equations. Theoretical Computer Science **599**, 64-78 (Sep 2015). https://doi.org/10.1016/j.tcs.2014.07.032

4. Fages, F., Le Guludec, G., Bournez, O., Pouly, A.: Strong Turing Completeness of Continuous Chemical Reaction Networks and Compilation of Mixed Analog-Digital Programs. In: CMSB'17: Proceedings of the fiveteen international conference on Computational Methods in Systems Biology. LNCS, vol. 10545, pp. 108-127. Springer-Verlag (Sep 2017). https://doi.org/10.1007/978-3-319-67471-1_7

5. Fages, F., Martinez, T., Rosenblueth, D., Soliman, S.: Influence networks compared with reaction networks: Semantics, expressivity and attractors. IEEE/ACM Transactions on Computational Biology and Bioinformatics (Feb 2018). https://doi.org/10.1109/TCBB.2018.2805686

6. Fages, F., Soliman, S.: Abstract interpretation and types for systems biology. Theoretical Computer Science **403**(1), 52-70 (2008). https://doi.org/10.1016/j.tcs.2008.04.024

7. Feinberg, M.: Mathematical aspects of mass action kinetics. In: Lapidus, L., Amundson, N.R. (eds.) Chemical Reactor Theory: A Review, chap. 1, pp. 1-78. Prentice-Hall (1977)

8. Gillespie, D.T.: Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry **81**(25), 2340-2361 (1977)

9. Kurtz, T.G.: Solutions of ordinary differential equations as limits of pure jump markov processes. Journal of applied Probability **7**(1), 49-58 (1970)
10. Malik-Sheriff, R.S., Glont, M., Nguyen, T.V.N., Tiwari, K., Roberts, M.G., Xavier, A., Vu, M.T., Men, J., Maire, M., Kananathan, S., Fairbanks, E.L., Meyer, J.P., Arankalle, C., Varusai, T.M., Knight-Schrijver, V., Li, L., Dueñas-Roca, C., Dass, G., Keating, S.M., Park, Y.M., Buso, N., Rodriguez, N., Hucka, M., Hermjakob, H.: BioModels — 15 years of sharing computational models in life science. Nucleic Acids Research 48(D1), D407–D415 (1 2020). https://doi.org/10.1093/nar/gkz1055, https://doi.org/10.1093/nar/gkz1055, gkz1055

11. Piho, P., Hillston, J.: Fluid approximation-based analysis for mode-switching population dynamics. ACM Trans. Model. Comput. Simul. 31(2) (feb 2021). https://doi.org/10.1145/3441680, https://doi.org/10.1145/3441680
10. Malik-Sheriff, R.S., Glont, M., Nguyen, T.V.N., Tiwari, K., Roberts, M.G., Xavier, A., Vu, M.T., Men, J., Maire, M., Kananathan, S., Fairbanks, E.L., Meyer, J.P., Arankalle, C., Varusai, T.M., Knight-Schrijver, V., Li, L., Dueñas-Roca, C., Dass, G., Keating, S.M., Park, Y.M., Buso, N., Rodriguez, N., Hucka, M., Hermjakob, H.: BioModels — 15 years of sharing computational models in life science. Nucleic Acids Research 48(D1), D407–D415 (1 2020). https://doi.org/10.1093/nar/gkz1055, https://doi.org/10.1093/nar/gkz1055, gkz1055

11. Piho, P., Hillston, J.: Fluid approximation-based analysis for mode-switching population dynamics. ACM Trans. Model. Comput. Simul. 31(2) (feb 2021). https://doi.org/10.1145/3441680, https://doi.org/10.1145/3441680