# Stochastic modeling for biotechnologies

## Anaerobic model AM2b

F. Campillo<sup>a</sup> — M. Chebbi<sup>b</sup> — S. Toumi<sup>c</sup>

**RÉSUMÉ.** Le modèle AM2b est classiquement représenté par un système d'équations différentielles. Toutefois ce modèle n'est valide qu'en grande population et notre objectif est d'établir plusieurs modèles stochastiques à différentes échelles. À l'échelle microscopique, on propose un modèle stochastic de saut pur que l'on peut simuler de façon exacte. Mais dans la plupart des situations ce genre de simulation n'est pas réaliste, et nous proposons des méthodes de simulation approchées de type poissonnien ou de type diffusif. La méthode de simulation de type diffusif peut être vue comme une discrétisation d'une équation différentielle stochastic. Nous présentons enfin de façon informelle un résultat de type loi des grands nombres/théorème central limite fonctionnelle qui démontre la convergence de ses modèles stochastiques vers le modèle déterministe initial.

**ABSTRACT.** The model AM2b is conventionally represented by a system of differential equations. However, this model is valid only in a large population context and our objective is to establish several stochastic models at different scales. At a microscopic scale, we propose a pure jump stochastic model that can be simulated exactly. But in most situations this exact simulation is not feasible, and we propose approximate simulation methods of Poisson type and of diffusive type. The diffusive type simulation method can be seen as a discretization of a stochastic differential equation. Finally, we formally present a result of law of large numbers and of functional central limit theorem which demonstrates the convergence of these stochastic models towards the initial deterministic models.

**MOTS-CLÉS** : modèle AM2b, processus de saut pur, équation différentielle ordinaire, approximation diffusion, équation différentielle stochastic.

**KEYWORDS** : AM2b Model, pure jump process, ordinary differential equation, diffusion approximation, stochastic differential equation.

## 1. Introduction

Stochastic models recently gain more credibility and numerical efficiency in many domains like chemistry [10], biotechnology [12] or system biology [16] where deterministic models have been extensively used. In the present paper we explain how a stochastic modeling approach deepens the insights allowed by the deterministic classical models.

We consider a model of a biotechnological process used for wastewater treatment. The principle of classical wastewater treatment is to couple the degradation of organic matter via the action of microorganisms and a decantation process. An alternative technology to decantation, implemented in membrane bioreactors (MBR), is based on the sieve principle, that is to say the passage of the effluent through calibrated pores of a physical membrane. The main drawback of this technology lies in the membrane performance degradation due to the fouling phenomenon. AM2b is a mathemati-

<footnote><sup>a</sup> Inria-France</footnote>

<footnote>Fabien.Campillo@inria.fr</footnote>

<footnote><sup>b</sup> ENIT-Laboratoire LAMSIN-University Tunis el Manar-Tunisie</footnote>

<footnote>mohsenchebbi@ymail.com</footnote>

<footnote><sup>c</sup> INSAT-Carthage University-Tunisie</footnote>

<footnote>salwa.toumi@lamsin.rnu.tn</footnote>
cal model of anaerobic membrane bioreactors developed in [1, 2], it is a variant of AM2 model [3], 2-steps Acidogenesis-Methanogenesis model, that includes a *soluble microbial products* (SMP) dynamics. The production and the degradation of SMP play an important role in the membrane fouling phenomenon.

First, we present the original *ordinary differential equation* (ODE) model of the AM2b. Then we introduce a pure jump Markov model of the same device and the associated exact Monte Carlo simulation method. Next we develop more efficient approximated simulation methods and we propose a *stochastic differential equation* (SDE) model of the AM2b. Finally a rescaling technique bridges these stochastic models with the original ODE model.

## 2. The ODE model

Define:

$S_1$: the organic matter,

$B_1$: the acidogenic biomass,

$S_2$: the volatile fatty acids (VFA),

$B_2$: the methanogenic biomass,

$S$: the soluble microbial products (SMP).

The state of the AM2b model is $x = (s_1, b_1, s_2, b_2, s)$ where $s_1, b_1, s_2, b_2, s$ are the concentrations in $S_1, B_1, S_2, B_2, S$. The AM2b model describes the dynamics of biological and anaerobic wastewater treatment, where substrate $s_1$ is degraded by a bacterial ecosystem $b_1$ to produce substrate $s_2$ and SMP $s$. Substrate $s_2$ is transformed by a consortium of bacteria $b_2$ into SMP $s$. SMP $s$ is also produced by the mortality of biomasses $b_1$ and $b_2$. SMP $s$ is degraded by $b_1$ to produce $s_2$. Schematically, AM2b model is:

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

On the right of this scheme is represented the membrane fouling model; the separation of matter is as follows: substrates $s_1$ and $s_2$ go through the membrane without retention (the size of their molecules is assumed to be smaller than pore diameter), biomasses $b_1$ and $b_2$ are retained by the membrane, and a fraction $\beta s$ of the SMP go through the membrane and leaves the reactor ((1 - $\beta$) $s$ will be considered as macromolecules).

The evolution of $x(t)$ is classically described by a system of ODEs:

$$ \dot{s}_1 = D_{in} [S_{1in} - s_1] - k_1 \mu_1(s_1) b_1, \quad (1a) $$

$$ \dot{b}_1 = [\mu_1(s_1) + \mu(s) - D_{dec} - D_{wit}] b_1, \quad (1b) $$

$$ \dot{s}_2 = D_{in} [S_{2in} - s_2] - k_2 \mu_2(s_2) b_2 + [c_{12} \mu_1(s_1) + c_{02} \mu(s)] b_1, \quad (1c) $$

$$ \dot{b}_2 = [\mu_2(s_2) - D_{dec} - D_{wit}] b_2, \quad (1d) $$

$$ \dot{s} = [c_{10} \mu_1(s_1) + D_{dec} - k_0 \mu(s)] b_1 + [c_{20} \mu_2(s_2) + D_{dec}] b_2 - Ms \quad (1e) $$

with:
<table>
  <tr>
    <td>k<sub>i</sub> degradation rate of s<sub>i</sub> by b<sub>i</sub>,</td>
    <td>M = &beta; D<sub>in</sub> + (1 &minus; &beta;) D<sub>wit</sub>,</td>
  </tr>
  <tr>
    <td>k<sub>0</sub> degradation rate of s by b<sub>1</sub>,</td>
    <td>&beta; fraction of SMP leaving the bioreactor,</td>
  </tr>
  <tr>
    <td>c<sub>12</sub> production rate of s<sub>2</sub> by b<sub>1</sub> from s<sub>1</sub>,</td>
    <td>D<sub>in</sub> dilution rate,</td>
  </tr>
  <tr>
    <td>c<sub>02</sub> production rate of s<sub>2</sub> by b<sub>1</sub> from s,</td>
    <td>D<sub>dec</sub> decay rate of biomass,</td>
  </tr>
  <tr>
    <td>c<sub>i0</sub> production rate of s by b<sub>i</sub> from s<sub>i</sub>,</td>
    <td>D<sub>wit</sub> withdrawal rate of biomass,</td>
  </tr>
  <tr>
    <td>S<sub>in</sub> input concentrations of s<sub>i</sub>,</td>
    <td>D<sub>out</sub> = D<sub>in</sub> &minus; D<sub>wit</sub> outflow rate of the bioreactor</td>
  </tr>
</table>

(i = 1, 2). The growth functions are :

$$
\mu_1(s_1) = m_1 \frac{s_1}{K_1 + s_1}, \quad \mu(s) = m \frac{s}{K+s}, \quad \mu_2(s_2) = m_2 \frac{s_2}{K_2 + s_2 + \frac{s_2^2}{K_3}}.
$$

$\mu_1$ and $\mu$ are of Monod type; $\mu_2$ is of Haldane type to model the phenomenon of inhibition of the
eventual accumulation of the VFA in the bioreactor during the methanogenesis, the major problem
of the anaerobic digestion.

System (1) could also be denoted :

$$
\dot{x}(t) = f(x(t)), \quad x(0) = x_0. \tag{2}
$$

Model (2) relies on the fact that the stochastic effects can be neglected, thanks to the law of large numbers, or at least can be averaged out. Although this level of description is sufficient for numerous applications of interest, it could be interesting to account for the stochastic nature of the process which could prominent at low biomasses concentrations or when small perturbations cumulate in the context of multi-species. Also, whereas the experimental results observed in well mastered laboratory conditions match closely the ODE theoretical behavior, a discrepancy may occur in operational conditions. In these cases, stochastic features may not be neglected.

3. Pure jump model

Based on the reactions featured in the AM2b model, we propose a AM2b model as a pure jump Markov process $(X(t))_{t \ge 0}$ which takes values in $\mathbb{R}_+^d$. We follow the approach proposed in [5]. Let $X = (S_1, B_1, S_2, B_2, S)$ be the concentrations in $S_1, B_1, S_2, B_2, S$ respectively. This process will encompass $J = 15$ reactions : each reaction $j$ is characterized by its rate function $\lambda_j(x)$ and its jump function $\nu_j(x)$. Then $X(t)$ is a time-homogeneous markov process and can be described as : being at $X(t) = x$ :

$$
X(t + \Delta t) = \begin{cases} x + \nu_j(x), & \text{with probability } \lambda_j(x) \Delta t + o(\Delta t) \text{ for } j = 1, \dots, 15, \\ x, & \text{with probability } 1 - \lambda_j(x) \Delta t + o(\Delta t) \end{cases} \quad (3)
$$

where :

$$
\lambda(x) \stackrel{\text{def}}{=} \sum_{j=1}^{15} \lambda_j(x), \qquad p_j(x) \stackrel{\text{def}}{=} \lambda_j(x) / \lambda(x).
$$

Reaction network

The reaction scheme of AM2b model is based on a modification of reaction scheme of the AM2
model [2, Chapter 2] by introducing the variable s for the SMP component. The production of SMP
is related to the degradation of s₁, s₂, and biomass decay b₁ and b₂. They are also degraded and
converted into s₂ and carbon dioxide CO₂ by the consortium b₁.
We define a network of J = 15 reactions classified in three sets :
1) Biological reactions :

$$
\begin{array}{l@{\quad}l@{\quad}l}
\text{Acidogenesis} + \text{SMP production} & k_1 S_1 + B_1 \xrightarrow{\lambda_1} 2 B_1 + c_{12} S_2 + c_{10} S + \kappa_1 \text{CO}_2, & \\
\text{Methanogenesis} + \text{SMP production} & k_2 S_2 + B_2 \xrightarrow{\lambda_2} 2 B_2 + c_{20} S + \kappa_2 \text{CH}_4, & \\
\text{SMP degradation} & k_0 S + B_1 \xrightarrow{\lambda_3} 2 B_1 + c_{02} S_2 + \kappa_3 \text{CO}_2, & \\
\text{SMP production from biomass decay} & B_1 \xrightarrow{\lambda_4} S, \quad B_2 \xrightarrow{\lambda_5} S &
\end{array}
$$

2) *Substrate inflow and substrate outflows* represented as reactions :

$$
\begin{array}{l@{\quad}l@{\quad}l}
\text{Substrate inflow} & \emptyset \xrightarrow{\lambda_6} S_1, & \emptyset \xrightarrow{\lambda_7} S_2, \\
\text{Substrate outflows} & S_1 \xrightarrow{\lambda_8} \emptyset, & S_2 \xrightarrow{\lambda_9} \emptyset, & S \xrightarrow{\lambda_{10}} \emptyset
\end{array}
$$

3) *Biomass and substrate withdrawal* described as reactions :

$$
\begin{array}{l@{\quad}l@{\quad}l@{\quad}l}
\text{Withdrawal} & S_1 \xrightarrow{\lambda_{11}} \emptyset, & B_1 \xrightarrow{\lambda_{12}} \emptyset, & S_2 \xrightarrow{\lambda_{13}} \emptyset, \\
& B_2 \xrightarrow{\lambda_{14}} \emptyset, & S \xrightarrow{\lambda_{15}} \emptyset &
\end{array}
$$

where the rate functions $\lambda_j$ are to be specified.

The second and third set of reactions are not biochemical reactions they just describe the inflows and outflows in the AM2b process. In reaction 10 only a proportion $\beta$ of the SMP goes through the membrane, and in reaction 13 a proportion $1 - \beta$ of the SMP is withdrawn.

## Pure jump Markov process

The AM2b reaction network described in previous section is translated into the pure jump Markov process $X(t)$ defined by (3) thanks to the stochastic law of mass action [16] and mass conservation : each reaction $j$ is described as an instantaneous jump $X(t) \rightarrow X(t) + \nu_j(X(t))$ occurring with intensity $\lambda_j(X(t))$ defined respectively by :

$$
\lambda_j(x) \stackrel{\text{def}}{=} N_j \tilde{\lambda}_j(x), \qquad \nu_j(x) \stackrel{\text{def}}{=} \left[x + \frac{1}{N_j} \tilde{\nu}_j\right]_+ - x \quad (4)
$$

where $[x]_+$ the orthogonal projection of $x$ onto $\mathbb{R}_+^5$, and :

<table>
  <thead>
    <tr>
      <th>j</th>
      <th>λ̃<sub>j</sub>(x)</th>
      <th>ν̃<sub>j</sub><sup>*</sup></th>
      <th>j</th>
      <th>λ̃<sub>j</sub>(x)</th>
      <th>ν̃<sub>j</sub><sup>*</sup></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>1</td>
      <td>μ<sub>1</sub>(s<sub>1</sub>) b<sub>1</sub></td>
      <td>(-k<sub>1</sub>, 1, c<sub>12</sub>, 0, c<sub>10</sub>)</td>
      <td>9</td>
      <td>D<sub>out</sub> s<sub>2</sub></td>
      <td>(0, 0, -1, 0, 0)</td>
    </tr>
    <tr>
      <td>2</td>
      <td>μ<sub>2</sub>(s<sub>2</sub>) b<sub>2</sub></td>
      <td>(0, 0, -k<sub>2</sub>, 1, c<sub>20</sub>)</td>
      <td>10</td>
      <td>β D<sub>out</sub> s</td>
      <td>(0, 0, 0, 0, -1)</td>
    </tr>
    <tr>
      <td>3</td>
      <td>μ(s) b<sub>1</sub></td>
      <td>(0, 1, c<sub>02</sub>, 0, -k<sub>0</sub>)</td>
      <td>11</td>
      <td>D<sub>wit</sub> s<sub>1</sub></td>
      <td>(-1, 0, 0, 0, 0)</td>
    </tr>
    <tr>
      <td>4</td>
      <td>D<sub>dec</sub> b<sub>1</sub></td>
      <td>(0, -1, 0, 0, 1)</td>
      <td>12</td>
      <td>D<sub>wit</sub> b<sub>1</sub></td>
      <td>(0, -1, 0, 0, 0)</td>
    </tr>
    <tr>
      <td>5</td>
      <td>D<sub>dec</sub> b<sub>2</sub></td>
      <td>(0, 0, 0, -1, 1)</td>
      <td>13</td>
      <td>D<sub>wit</sub> s<sub>2</sub></td>
      <td>(0, 0, -1, 0, 0)</td>
    </tr>
    <tr>
      <td>6</td>
      <td>D<sub>in</sub>S<sub>1in</sub></td>
      <td>(1, 0, 0, 0, 0)</td>
      <td>14</td>
      <td>D<sub>wit</sub> b<sub>2</sub></td>
      <td>(0, 0, 0, -1, 0)</td>
    </tr>
    <tr>
      <td>7</td>
      <td>D<sub>in</sub>S<sub>2in</sub></td>
      <td>(0, 0, 1, 0, 0)</td>
      <td>15</td>
      <td>(1 - β) D<sub>wit</sub> s</td>
      <td>(0, 0, 0, 0, -1)</td>
    </tr>
    <tr>
      <td>8</td>
      <td>D<sub>out</sub> s<sub>1</sub></td>
      <td>(-1, 0, 0, 0, 0)</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
  </tbody>
</table>

The coefficients $N_j$ are rescaling parameters which characterize the size of the jump in the reaction $j$. Large $N_j$ corresponding to frequent and small jumps. These scale parameters $N_j$ do not act on the mean values of the increments but on their variances, large $N_j$ will correspond to small variances. We can assume that the $N_j$'s range from $10^4$ to $10^9$. When a reaction involves only substrate molecules the corresponding $N_j$'s range from $10^7$ to $10^9$. When a reaction involves only bacteria the corresponding $N_j$'s range from $10^4$ to $10^6$. A standard case could be :

$$
N_j = \begin{cases} 10^6 & \text{for } j = 4, 5, 12, 14, \\ 10^5 & \text{for } j = 1, 2, 3, \\ 10^4 & \text{for } 6 \le j \le 15. \end{cases} \qquad (5)
$$
We define :

$$
F(x) \stackrel{\text{def}}{=} \sum_{j=1}^{15} \lambda_j(x) \nu_j(x) \qquad (6)
$$

note that $F(x)$ depends on the $N_j$'s but $F(x) = \sum_{j=1}^{15} \frac{\lambda_j(x)}{N_j} \times N_j \nu_j(x)$ so that, for large $N_j$'s :

$$
F(x) \approx \sum_{j=1}^{15} \tilde{\lambda}_j(x) \tilde{\nu}_j \qquad (7)
$$

and :

$$
f(x) \stackrel{\text{def}}{=} \sum_{j=1}^{15} \tilde{\lambda}_j(x) \tilde{\nu}_j \qquad (8)
$$

is exactly the second hand term of Equation (2).

Note that :

$$
\nu_j(x) \stackrel{\text{def}}{=} [x + \frac{1}{N_j} \tilde{\nu}_j]_+ - x
$$

We have then for each component : $\forall i = 1, \dots, 5$ :

$$
- \text{ if } x_i + \frac{1}{N_j} \{\tilde{\nu}_j\}_i \ge 0 \text{ then } \{\nu_j(x)\}_i = \frac{1}{N_j} \{\tilde{\nu}_j\}_i;
$$

$$
- \text{ if } x_i + \frac{1}{N_j} \{\tilde{\nu}_j\}_i \le 0 \text{ then } \{\nu_j(x)\}_i = -x_i.
$$

In this last case $|\{\nu_j(x)\}_i| = x_i \le \frac{1}{N_j}|\{\tilde{\nu}_j\}_i|$ and then $\exists C$ such that for all $j$ :

$$
|N_j \nu_j(x)| \le C |\tilde{\nu}_j|, \quad \forall x \in \mathbb{R}_+^5. \qquad (9)
$$

Representation of the process $X(t)$

Trajectories of the process $X(t)$ can be exactly simulated according to the so-called *Stochastic Simulation Algorithm* (SSA) popularized by Gillespie [7]. Starting for an initial condition $X(0)$, we simulate the next time of jump $T_1$ from an exponential distribution of intensity $\lambda(X(0))$, then we draw the index $j_0$ of the reaction from the distribution $p_j(X(0))$, $j = 1, \dots, J$. Then $X(T_1) = X(0) + \nu_{j_0}(X(0))$, and so on. Hence the Markov process $X(t)$ can be represented as :

$$
X(t) = X(0) + \sum_{j=1}^{15} \int_{0}^{t} \int_{\mathbb{R}_{+}} \nu_j(X(s^{-})) 1_{[0, \lambda_j(X(s^{-}))]}(v) N_j(ds, dv) \quad (10)
$$

where the $N_j(ds, dv)$ are independent Poisson random measures (also independent from $X(0)$) of
intensity measure $ds \times dv$, the Lebesgue measure on $\mathbb{R}_+^2$.

Note that (10) is expressed in law and is a generalization of the expression obtained in the case of
constant jumps (see [17, theorem 4.1 chapter 6]), namely we used :

$$
\sum_{j=1}^{J} \int_{0}^{t} \mathcal{P}_{j}(ds \times \lambda_{j}(X_{s^{-}})) \nu_{j}(X_{s^{-}}) \stackrel{\mathcal{L}}{=} \sum_{j=1}^{J} \int_{0}^{t} \int_{0}^{\infty} 1_{[0, \lambda_{j}(X_{s^{-}})]}(u) \nu_{j}(X_{s^{-}}) \mathcal{N}_{j}(ds, du),
$$

for more justification, see [18].

The law of the process $X(t)$ is characterized by its infinitesimal generator define by :

$$
\mathcal{A}_{\phi}(x) = \lim_{t \to 0} \frac{\mathbb{E}[\phi(X(t)) - \phi(X(0)) | X(0) = x]}{t},
$$
we can easily check that :

$$ \mathcal{A}_{\phi}(x) = \lambda(x) \int_{\mathbb{R}_{+}^{5}} [\phi(y) - \phi(x)] \rho(x, dy) = \sum_{j=1}^{15} \lambda_{j}(x) [\phi(x + \nu_{j}(x)) - \phi(x)] $$

for all continuous with compact support function $\phi$, where :

$$ \rho(x, dy) \stackrel{\text{def}}{=} \sum_{j=1}^{15} p_j(x) \delta_{x+\nu_j(x)}(dy) $$

We introduced the compensated Poisson measure :

$$ \tilde{\mathcal{N}}_j(ds, dv) \stackrel{\text{def}}{=} \mathcal{N}_j(ds, dv) - ds \, dv. $$

From (10) we get :

$$ X(t) = X(0) + \int_{0}^{t} F(X(s)) \, ds + M_{t} \quad (11) $$

where $F(x)$ is defined by (6) and :

$$ M_t = \sum_{j=1}^{15} M_t^j, \text{ with } M_t^j \stackrel{\text{def}}{=} \int_0^t \int_{\mathbb{R}_+} \nu_j(X(s^{-})) 1_{[0, \lambda_j(X(s^{-}))]}(v) \tilde{\mathcal{N}}_j(ds, dv) $$

$M_t^j$ is a square integrable martingale with predictable quadratic variation :

$$ \langle M^j \rangle_t = \int_0^t \lambda_j(X(s)) |\nu_j(X(s))|^2 ds \quad (12) $$

In Equation (11) the process $X(t)$ is written as the sum of :

- a drift term $\int_0^t F(X(s)) ds$ which describes the deterministic dynamics of $X(t)$ and, for large $N_j$'s, it coincides with the model AM2b (2);

- a martingale term $M_t$ with quadratic variation (12), this quadratic variation characterizes the variance of the difference between the stochastic model and the deterministic model. In particular, this term converges to 0 when $N_j \to \infty$.

# 4. Discrete time approximations

The Stochastic Simulation Algorithmic (SSA) simulates each reaction of the ecosystem asynchronously in time. In many situations this detailed simulation is too cumbersome, this is why synchronous discrete time approximations have been proposed. Let :

$$ t_m \stackrel{\text{def}}{=} m \Delta t $$

for $\Delta t > 0$ fixed.

## Poisson approximation

We construct an approximation $(\tilde{X}(t_m))_{m \ge 1}$. On the interval $[t_m, t_{m+1})$ suppose that the different rate functions are constant with :

$$ \lambda_j(\tilde{X}(t) \approx \lambda_j(\tilde{X}(t_m)) $$
so that the $J = 15$ reactions are independent and occur at constant rates $\lambda_j(\tilde{X}(t_m))$, that is the occurrence of reaction of type $j$ is a Poisson process of intensity $\lambda_j(\tilde{X}(t_m))$. Hence, on the time interval $[t_m, t_{m+1})$ the number of reactions of type $j$ follows a Poisson distribution of parameter $\Delta t \lambda_j(\tilde{X}(t_m))$. We obtain the following approximation called $\tau$-leaping :

$$ \tilde{X}(t_{m+1}) = [\tilde{X}(t_m) + \sum_{j=1}^{15} \nu_j(\tilde{X}(t_m)) P_{j,m}]_+ \quad (13) $$

where for each $m$, the $P_{j,m}$ are independent Poisson distributed variables with parameter $\Delta t \lambda_j(\tilde{X}(t_m))$.

## Diffusion approximation

The Poisson distribution with parameter $\Delta t \lambda_j(\tilde{X}(t_m))$ can be approximated by a normal distribution of mean $\Delta t \lambda_j(\tilde{X}(t_m))$ and variance $\Delta t \lambda_j(\tilde{X}(t_m))$. From (13) we get :

$$ \tilde{\xi}(t_{m+1}) = \left[ \tilde{\xi}(t_m) + \sum_{j=1}^{15} \nu_j(\tilde{\xi}(t_m)) \left( \Delta t \lambda_j(\tilde{\xi}(t_m)) + \sqrt{\Delta t \lambda_j(\tilde{\xi}(t_m))} w_{j,m} \right) \right]_+ $$

where $w_{j,m}$ are independent $N(0, 1)$ random variables. This last equation can be rewritten :

$$ \tilde{\xi}(t_{m+1}) = \left[ \tilde{\xi}(t_m) + F(\tilde{\xi}(t_m)) \Delta t + \sum_{j=1}^{15} \frac{1}{\sqrt{N_j}} g_j(\tilde{\xi}(t_m)) [W_j(t_{m+1}) - W_j(t_m)] \right]_+ \quad (14) $$

where $W_j(t)$ are independent standard Brownian motions so that $W_j(t_{m+1}) - W_j(t_m)$ are independent and $N(0, \Delta t)$, and :

$$ g_j(x) \stackrel{\text{def}}{=} \sqrt{N_j \lambda_j(x)} \nu_j(x) = \sqrt{\tilde{\lambda}_j(x)} N_j \nu_j(x). $$

and according to (9) :

$$ |g_j(x)| \le C \sqrt{\tilde{\lambda}_j(x)} \tilde{\nu}_j. $$

## Stochastic differential equation

Equation (14), is an Euler-Maruyama approximation of the following SDE :

$$ d\xi(t) = F(\xi(t)) dt + \sum_{j=1}^{15} \frac{1}{\sqrt{N_j}} g_j(\xi(t)) dW_j(t). \quad (15) $$

# 5. Scales and asymptotics

## A first analysis

Suppose that $N_j = N$ for all $j$. We can prove rigorously [13, 14] that :

$$ \sup_{0 \le t \le T} |X(t) - x(t)| \xrightarrow{N \to \infty} 0 \quad (16) $$

in $L^2(\Omega)$ or in probability. In (16) we can also replace $X(t)$ by $\xi(t)$. So under specific conditions, when the population sizes are large and so the number of reactions is, i.e. $N$ large, the ODE model (2) is valid.
At an intermediate scale, a functional central limit theorem states that the process $\sqrt{N}(X(t) - x(t))$ can be approximated in law by $\sum_{j=1}^{15} \int_{0}^{t} g_j(x(s))d\tilde{W}_j(s)$ where the $\tilde{W}_j(s)$ are independent standard Brownian motions, that is formally:

$$X(t) \approx x(t) + \frac{1}{\sqrt{N}} \sum_{j=1}^{15} \int_{0}^{t} g_j(x(s)) d\tilde{W}_j(s).$$

This proves that the SDE model (15) is valid for $N$ large but also for $N$ moderately large.

## A second analysis : hybrid model

In many situation ODE and SDE models are not valid. This is the case when one of the bacterial population is present in low concentration but still affects the global dynamic of the process. This so-called molecular randomness may influence the global dynamic even when the population sizes are not so small [6]. In this case we may adopt an hybrid approach.

We present an example where we separate the dynamics of the substrates from the dynamics of the biomasses, i.e. separate substrate type reactions from biomass type reactions, then to describe the first ones as a system of ODE's and to describe the second ones as a pure jump Markov process. We obtain a system of ODE's describing the continuous evolution of the substrates and the SMP concentrations:

$$\dot{s}_1 = D(S_{1\text{in}} - s_1) - k_1 \mu_1(s_1) B_1,$$

$$\dot{s}_2 = D(S_{2\text{in}} - s_2) - k_2 \mu_2(s_2) B_2 + (c_{12} \mu_1(s_1) + c_{02} \mu(s)) B_1,$$

$$\dot{s} = (c_{10} \mu_1(s_1) + D_{\text{dec}} - k_0 \mu(s)) B_1 + (c_{20} \mu_2(s_2) + D_{\text{dec}}) B_2 - Ms$$

coupled to a 2-dimensional pure jump process describing the discrete evolution of the biomasses concentrations:

<table><thead><tr><th colspan="2">jump</th><th colspan="2">rate</th><th colspan="2">jump</th><th colspan="2">rate</th></tr></thead><tbody><tr><td>B<sub>1</sub> → B<sub>1</sub> + δ<sub>1</sub></td><td>μ<sub>1</sub>(s<sub>1</sub>)</td><td>B<sub>1</sub>/δ<sub>1</sub></td><td></td><td>B<sub>2</sub> → B<sub>2</sub> + δ<sub>2</sub></td><td>μ<sub>2</sub>(s<sub>2</sub>)</td><td>B<sub>2</sub>/δ<sub>2</sub></td><td></td></tr><tr><td>B<sub>1</sub> → B<sub>1</sub> + δ<sub>1</sub></td><td>μ(s)</td><td>B<sub>1</sub>/δ<sub>1</sub></td><td></td><td>B<sub>2</sub> → B<sub>2</sub> - δ<sub>2</sub></td><td>D<sub>dec</sub></td><td>B<sub>2</sub>/δ<sub>2</sub></td><td></td></tr><tr><td>B<sub>1</sub> → B<sub>1</sub> - δ<sub>1</sub></td><td>D<sub>dec</sub></td><td>B<sub>1</sub>/δ<sub>1</sub></td><td></td><td>B<sub>2</sub> → B<sub>2</sub> - δ<sub>2</sub></td><td>D<sub>1</sub></td><td>B<sub>2</sub>/δ<sub>2</sub></td><td></td></tr><tr><td>B<sub>1</sub> → B<sub>1</sub> - δ<sub>1</sub></td><td>D<sub>1</sub></td><td>B<sub>1</sub>/δ<sub>1</sub></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

# 6. Simulation

We now present some simulations of the ODE and SDE models using the values for the parameters of the AM2b model defined in [1, Table 1 p. 6].
<figure><img src="image_2.png" /><figcaption>FIGURE 1: Simulation of diffusion approximation (14) with the parameters (5). Time evolution of the organic matter concentration (left), time evolution of the acidogenic biomass concentration (right). Phase portrait acidogenic biomass/organic matter concentrations (below).</figcaption></figure>

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

# 7. Conclusion

We show that the ODE model (1) contains all the ingredients that can be used to establish a pure jump Markov model. The ODE model is valid only in large population scales. The SDE model
<figure><img src="image_4.png" /><figcaption>FIGURE 2 : Simulation of diffusion approximation (14) with the parameters (5). Time evolution of the AGV concentration (left), time evolution of the methanogenic biomass concentration (right). Phase portrait methanogenic biomass/AGV concentrations (below).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3 : Simulation of diffusion approximation (14) with the parameters (5). Time evolution of the SMP concentration (left). Phase portrait acidogenic biomass/SMP concentrations (right).</figcaption></figure>

is valid in large population scales but also in “moderately” large population scales. For small population scales, only the pure jump Markov model is valid. The most promising approach is to hybridize an ODE model for substrate dynamics coupled with a stochastic model for the dynamic of the biomasses.

## Acknowledgments

This work was partially funded by the NuWat LIRIMA Inria project.

## 8. Bibliographie

[1] B. Benyahia, T. Sari, B. Cherki, and J. Harmand, « Anaerobic membrane bioreactor modeling in the presence of Soluble Microbial Products (SMP) – the Anaerobic Model AM2b », *Chemical Engineering Journal*, 228(0):1011-1022, 2013.

[2] B. Benyahia, « Modélisation et observation des bioprocédés à membranes : application à la digestion anaérobie », *PhD thesis*, Université de Tlemcen et Université de Montpellier 2, 2012.
[3] O. Bernard, Z. Hadj-Sadock, D. Dochain, A. Genovesi, and J.-P. Steyer, « Dynamical model development and parameter identification for an anaerobic wastewater treatment process », *Biotechnology and Bioengineering*, 75:424-438, 2001.

[4] F. Campillo, M. Chebbi, and S. Toumi, « Stochastic modeling of the anaerobic model AM2b-Models at different scales », *CARI'13, 11-14 October, 2016, Hammamet, Tunisie*.

[5] F. Campillo, M. Joannides, and I. Larramendy-Valverde « Stochastic modeling of the chemo-stat », *Ecological Modelling*, 222(15):2676-2689, 2011.

[6] F. Campillo and C. Lobry, « Effect population size in a Predator-Prey model », *Ecological Modelling*, 246:1-10, 2012.

[7] D. T. Gillespie, « Exact stochastic simulation of coupled chemical reactions », *Journal of Physical Chemistry*, 81(25):2340-2361, 1977.

[8] D. T. Gillespie, « Approximate accelerated stochastic simulation of chemically reacting systems », *Journal of Physical Physics*, 115(4):1716-1733, 2001.

[9] D. T. Gillespie, and L. R. Petzold, « Improved leap-size selection for accelerated stochastic simulation », *Journal of Physical Physics*, 119:8229-8234, October 2003.

[10] D. J. Higham, « Modeling and simulating chemical reactions », *SIAM Review*, 50(2):347-368, 2008.

[11] G. Kersting, and F. C. Klebaner, « Sharp conditions for nonexplosions and explosions in Markov jump processes », *The Annals of Probability*, 23(1):268-272, 1995.

[12] S. Kops, K. Gernaey, O. Thas, and P. A. Vanrolleghem, « The modelling of noise processes in stochastic differential equations : Application to biotechnological processes », *Proceedings 7th IFAC Conference on Computer Applications in Biotechnology CAB7. Osaka, Japan, May 31-June 4*, pages 67-72, 1998.

[13] T. G. Kurtz, « Solutions of ordinary differential equations as limits of pure jump Markov processes », *Journal of Applied Probability*, 7(1):49-58, 1970.

[14] T. G. Kurtz, « Limit theorems for sequences of jump Markov processes approximating ordinary differential processes », *Journal of Applied Probability*, 8:344-356, 1971.

[15] M. Ullah, and O. Wolkenhauer, « Stochastic Approaches in Systems Biology », *Springer*, 2011.

[16] D. J. Wilkinson, « Stochastic Modelling for Systems Biology », *Chapman & Hall/CRC*, 2006.

[17] Ethier, Stewart N. and Kurtz, Thomas G., « Markov processes, characterization and convergence », *New York : John Wiley Sons Inc.*, 1986.

[18] Grigelionis B, « On the representation of integer-valued random measures by means of stochastic integrals with respect to the Poisson measure », *Lit. Math.J*, 11.1 (1971) : 93-108.