On real-valued SDE and nonnegative-valued SDE population models with demographic variability

Abstract

Population dynamics with demographic variability is frequently studied using discrete random variables with continuous-time Markov chain (CTMC) models. An approximation of a CTMC model using continuous random variables can be derived in a straightforward manner by applying standard methods based on the reaction rates in the CTMC model. This leads to a system of Itô stochastic differential equations (SDEs) which generally have the form $d \mathbf {y} = \varvec{\mu } \, dt + G \, d \mathbf {W},$ where $\mathbf {y}$ is the population vector of random variables, $\varvec{\mu }$ is the drift vector, and G is the diffusion matrix. In some problems, the derived SDE model may not have real-valued or nonnegative solutions for all time. For such problems, the SDE model may be declared infeasible. In this investigation, new systems of SDEs are derived with real-valued solutions and with nonnegative solutions. To derive real-valued SDE models, reaction rates are assumed to be nonnegative for all time with negative reaction rates assigned probability zero. This biologically realistic assumption leads to the derivation of real-valued SDE population models. However, small but negative values may still arise for a real-valued SDE model. This is due to the magnitudes of certain problem-dependent diffusion coefficients when population sizes are near zero. A slight modification of the diffusion coefficients when population sizes are near zero ensures that a real-valued SDE model has a nonnegative solution, yet maintains the integrity of the SDE model when sizes are not near zero. Several population dynamic problems are examined to illustrate the methodology.

Appendices

Appendix 1: a nonnegative SDE model with smooth diffusion coefficients

The nonnegative SDE model described in Sect. 2 does not have Lipschitz continuous diffusion coefficients at zero. A nonnegative model with smooth diffusion coefficients is described in this appendix where the reaction rates are modified slightly differently than in the nonnegative SDE model of Sect. 2. To describe the model, the changes $\lambda _{k,i}$ are considered for $i=1,\ldots , N$ for each reaction rate $k= 1,2, \ldots , M$. Given small $\delta >0$, recall that

$$\begin{aligned}&L_k \subseteq \{1,2, \ldots , N\} \;\; \text {is the set of all indices} \;\; i \in \{1,2, \ldots , N\} \;\; \text {where} \;\; \lambda _{k,i} \ne 0. \end{aligned}$$

Then, for this new SDE model, the kth rate is modified as

$$\begin{aligned} {{\hat{r}}}_k^{\delta }(t, \mathbf {y}(t)) = {{\hat{r}}}_k(t, \mathbf {y}(t)) \prod _{i \in L_k} (\phi _i^{\delta } (\mathbf {y}(t)))^2 \end{aligned}$$

where $\phi _i^{\delta } (\mathbf {y}(t))$ is defined in Eq. (10). In particular, if $y_i(t) \ge \delta$ for all $i \in L_k$, then ${{\hat{r}}}_k^{\delta }(t, \mathbf {y}(t)) = {\hat{r}}_k(t, \mathbf {y}(t)).$ Similar to the earlier nonnegative SDE, the modified diffusion coefficients are given by

$$\begin{aligned} (G^{\delta }(t, \mathbf {y}(t)))_{i,k} = \lambda _{k,i} \sqrt{ {{\hat{r}}}_k^{\delta }(t,\mathbf {y}(t))} \;\; \text {for} \;\; 1 \le i \le N \end{aligned}$$

for $k=1, \ldots , M$. Furthermore, the conservation property of Proposition 1 holds for both nonnegative SDE models and, computationally, the two nonnegative SDE models give similar results. For example, in the Michaelis–Menten worst-case scenario exit-time problem (described in Sect. 3.5) with $K=4$, $y_1(0)=5$, and $\delta = 0.1$, the calculated exit times for the real-valued SDE, the nonnegative-valued SDE, the complex CLE method, and the CTMC models give, respectively, mean values of 2.80, 2.92, 2.90, 3.17 and standard deviations of 0.72, 0.77, 0.80, 0.84 for 5000 sample paths. The nonnegative-valued SDE model described in this appendix gives a mean exit time of 3.14 with a standard deviation of 0.85 which compares well with the other models.

Appendix 2: additional Michaelis–Menten computations for $K=4$

Continuing the example from Sect. 3.5, means and standard deviations are compared for the substrate and enzyme molecules using the real-valued and nonnegative-valued SDE models and CLE method for $K=4$. The number of substrate molecules, $y_1(t)$, and the number of free enzyme molecules, $y_2(t)$, are calculated when $y_1(0)=1$, $y_2(0)= 3$, $y_3(0)=0$, $k_2=1$, $k_3=1$, $k_4=5$, $K=4$ up to time $t=5$. The Euler–Maruyama method is used to approximate the real-valued and nonnegative-valued SDEs with a time step of 0.0005, 1000 sample paths, and $\delta =0.1$. The computational results for average population sizes and standard deviations are displayed in Fig. 9 along with the results of the complex CLE method. The CLE method and the results of the real-valued and nonnegative-valued SDE models generally agree well with the CTMC model for the substrate molecules, but there is disagreement with the CTMC model for free enzyme molecules when the mean is less than two. In particular, the calculated standard deviations for the real-valued and nonnegative-valued SDE models differ from CTMC model by as much as 20% with a closer agreement between the complex CLE method and the CTMC model ($\approx 10\%)$.

Fig. 9

Average of $y_1(t)$ (top left panel), average of $y_2(t)$ (top right panel), standard deviation in $y_1(t)$ (bottom left panel), and standard deviation in $y_2(t)$ (bottom right panel) for an enzyme kinetics problem. The solid curves represent the results for the CTMC model, stars for the complex CLE method, circles for the nonnegative-valued SDE, and pluses for the real-valued SDE