A Backward Technique for Demographic Noise in Biological Ordinary Differential Equation Models

: Physical systems described by deterministic differential equations represent idealized situations since they ignore stochastic effects. In the context of biomathematical modeling, we distinguish between environmental or extrinsic noise and demographic or intrinsic noise, for which it is assumed that the variation over time is due to demographic variation of two or more interacting populations (birth, death, immigration, and emigration). The modeling and simulation of demographic noise as a stochastic process affecting units of populations involved in the model is well known in the literature, resulting in discrete stochastic systems or, when the population sizes are large, in continuous stochastic ordinary differential equations and, if noise is ignored, in continuous ordinary differential equation models. The inverse process, i.e., inferring the effects of demographic noise on a natural system described by a set of ordinary differential equations, is still an issue to be addressed. With this paper, we provide a technique to model and simulate demographic noise going backward from a deterministic continuous differential system to its underlying discrete stochastic process, based on the framework of chemical kinetics, since demographic noise is nothing but the biological or ecological counterpart of intrinsic noise in genetic regulation. Our method can, thus, be applied to ordinary differential systems describing any kind of phenomena when intrinsic noise is of interest.


Introduction
Physical systems are usually modeled by differential equations, either ordinary (ODEs) or delay (DDEs) or partial (PDEs). These models represent idealized situations, as they ignore stochastic effects. By incorporating random elements in a differential equation model, a system of stochastic differential equations (SDEs), either stochastic ordinary (SODEs) or stochastic delay (SDDEs) or stochastic partial differential equations (SPDEs), respectively, arises. Stochastic models arise if we assume that a physical system operates in a noisy environment or if we wish to study the noisy behavior in the systems themselves, for example the intrinsic variability of interaction between individuals of two or more competing species. In the first case, deterministic differential equation models are enhanced by the introduction of some kind of noise processes modeling environmental fluctuations and lead to stochastic differential equation models. In the context of biomathematical modeling, we should distinguish between environmental or external noise and demographic or internal noise, for which it is assumed that the variation over time is due to demographic variation of two or more interacting populations (birth, death, immigration, and emigration) and not due to environmental variation. The populations could represent predator-prey interaction, competition, mutualism, or a single population divided into sub-populations of susceptible and infected individuals. On the other hand, external noise can arise from random fluctuations of one or more model parameters around some given mean value or from stochastic fluctuations of the system around any of the equilibrium points, in which cases the resulting stochastic systems have multiplicative noise. Fluctuations in experimental data, instead, result in additive noise stochastic systems.
An important point with stochastic models is also the understanding of the effects of randomness as drivers for stability. The stability theory for stochastic (either ordinary, delay, or partial) differential equations can be much richer than that for deterministic differential equations. The works by Mao [1][2][3] showed that environmental noise can have the effect of stabilization of an unstable system allowing to state a number of theorems on stochastic asymptotic stability where the Lyapunov functions may not be positive definite nor the generating operators of Lyapunov functions positive definite. The same author showed that another important feature of stochastic models in population dynamics is that environmental noise can suppress deterministic explosion.
Other than Biomathematics, randomness plays a crucial role when we try to model and simulate interactions among molecular species in single cells in Computational Cell Biology. In order to understand the functioning of an organism, the network of interactions between genes, mRNAs, miRNAs, proteins and other molecules needs to be addressed. Genetic (Gene) Regulation, relying on chemical kinetics theory, is the specific field that comes into play. In this case, cells are seen as biochemical reactors in which low reactant numbers can lead to significant statistical fluctuations in molecule numbers and reaction rates. These fluctuations are intrinsic: they are determined by the reaction parameters and connectivity of the underlying biochemical networks. In fact, in single cells, "molecular interactions produce stochastic fluctuations that may dominate cellular dynamics" [4]. Therefore, the modeling of molecular interactions in single cells requires the use of stochastic biochemical systems, which give a probabilistic description of cellular dynamics, and suitable simulation methods must be used.
Moreover, dynamical features such as oscillations or bi-stability cannot be observed as the effects of external (environmental) noise and the appropriate way to model them is to assume that intrinsic noise is affecting the system. The discrete stochastic methods we use do not follow individual molecules over time: they track only total molecular numbers. This leads to algorithms that describe the time evolution of a discrete nonlinear Markov process [5,6]. Gillespie first conceived a rigorous mathematical modeling framework through the Stochastic Simulation Algorithm [7]. Unfortunately, the original Gillespie algorithm was computationally very intensive; thus, refinements of it were introduced by many authors in order to reduce the computation time in simulations [8][9][10][11][12]. Gillespie [8] introduced the τ-leap and midpoint τ-leap methods where the sampling of most likely reactions was taken from a Poisson distribution, Tian and Burrage [9] proposed the sampling from a Binomial distribution and Cao et al. [10][11][12] introduced new strategies for the choice of the step size in τ-leap methods using Poisson τ-leaping.
A definite acceleration of the simulation process, though, can only be done by turning to the continuous stochastic system obtained by considering the ∆t → 0 limit of the discrete Markov process, leading to a system of stochastic ordinary differential equations (SODEs).
The modeling and simulation of demographic (intrinsic) noise in biological, ecological and biochemical systems have been addressed by many authors, among whom are Smith and Trevino [13], who studied demographic stochasticity in bacteriophage parasitism of a host bacteria, Mandal et al. [14], who addressed the problem of demographic noise in phytoplankton allelopathy, and Carletti et al. [15], who constructed a stochastic model for the regulation of gene PTEN implied in malignant brain tumours. All of these approaches have one feature in common, i.e., they start from considering the phenomenon from its intrinsic discrete nature, then go onward to the continuous stochastic setting, and finally result in a deterministic ODEs system. What is still lacking in the literature is an inverse method, i.e., a method to understand the impact of intrinsic noise when the only known dynamics is given by an ordinary differential system. With this paper, we intend to fill this gap by providing a technique that can reveal the effects of intrinsic/demographic noise starting from an ODEs model and going backwards to the discrete stochastic system. Our technique relies on the stochastic modeling and simulation of Gene Regulatory Networks (GRNs) by considering the involved biological species as they were reactants of a chemically reacting system.

Materials and Methods
In order to elucidate the novel technique for finding the effects of intrinsic noise starting from an ordinary differential system, we consider a biological model for phage-bacteria interaction in an open marine environment. The choice is due to the fact that all the analytical properties of the model are well known from the paper by Beretta and Kuang [16] and that, when we address intrinsic noise in the form of demographic noise, we are interested in stochasticity in birth and death rates, as well as proliferation and migration of biological or ecological populations.
Our technique is based on the theory developed in the field of chemical kinetics, according to which an ODE model can be seen as the deterministic counterpart of the stochastic differential system arising from the ∆t → 0 limit of the discrete Markov process modeling the interaction between species, assuming them to be expressed by number of units instead of by concentration of populations. In our framework, this theory can be seen as the "direct" method to approach demographic noise in natural processes, assuming them to be inherently stochastic [5,[17][18][19]. The method we introduce in this paper, instead, is meant to give insights to demographic noise in an "inverse" or backward way, i.e., assuming that only the deterministic behavior of the process is known and modeled by ordinary differential equations.
The outline of this article is the following: • In Section 2.1 we provide a background to the modeling of gene regulation, pointing out the reason why it is an intrinsically stochastic mechanism and describing the most common simulation techniques. • In Section 2.2 we illustrate an example of modeling and simulation of a GRN in bacteriophage λ, based on the paper by Goutsias [20] and the simulation technique described in [17].

•
In Section 2.3 we propose the novel technique through which we can study the effects of demographic noise in the ODEs model by Beretta and Kuang [16], going backwards to the discrete stochastic process and the continuous stochastic system (characterised by SODEs) underlying the deterministic model.

Stochastic Modeling of Genetic Regulation
In molecular biology, the general transfers describe the normal flow of biological information: DNA can be copied to DNA (DNA replication), DNA information can be copied into mRNA (messenger RNA), (transcription), and proteins can be synthesized using the information in mRNA as a template (translation).
A Genetic Regulatory Network is a collection of DNA segments in a cell that interact with each other and with other substances in the cell, thereby governing the rates at which genes in the network are transcribed into mRNA. Regulation of gene expression, or gene regulation, refers to the cellular control of the amount and timing of changes to the appearance of the functional product of a gene.
The mechanisms that act on gene expression, i.e., the ability of a gene to produce a biologically active protein, are based on biochemical processes that are inherently stochastic and are much more complex in eukaryotes than in prokaryotes. A major difference is that in eukaryotes there is the presence of the nuclear membrane, which prevents the simultaneous transcription and translation that occurs in prokaryotes [21]. Whereas, in prokaryotes, control of transcriptional initiation is the major point of regulation, in eukaryotes the regulation of gene expression can be controlled from many different points. In eukaryotes, transcriptional regulation tends to involve combinatorial interactions between several transcription factors. By "regulation" of transcription, we mean when transcription occurs and how much RNA is created.
The simulation of a discrete chemical reacting system must be done through Stochastic Simulation Algorithms (SSA), in which we assume low and integer molecular numbers. The first SSA is due to Gillespie [7] and is a rigorous procedure for numerically simulating the evolution of a set of chemical reactions in a (well-stirred) homogeneous chemical reacting system, by taking into account the randomness inherent in such a system. If S 1 , S 2 , . . . , S N represent N interacting molecular species through M reaction channels, then the evolution of such a system is characterized by a discrete nonlinear Markov process in which a vector X(t) of dimension N, representing numbers (integer values) of the N molecular species at time t, is evolved through time. The state vector X(t) is a discrete jump Markov process, where the time evolution equation describing the probability that, given X(t 0 ) = x 0 , then X(t) = x, i.e., P(x, t|x 0 , t 0 ) is called the Chemical Master Equation (CME) and can be written as This discrete parabolic partial differential equation is generally too difficult to solve (either analytically or numerically) and different methods are needed to simulate its solution.
Any method to solve a system of chemical reactions is uniquely characterized by two sets of quantities: the stoichiometric vectors ν 1 , . . . , ν M which represent, in turn, the update of the numbers of molecules in the system if reactions R 1 , . . . , R M occur, respectively, and a set of M propensity functions a 1 (X(t)), . . . , a M (X(t), that describe the relative probabilities of reactions R 1 , . . . , R M occurring, respectively.
The SSA is an essentially exact variable time step procedure, where "exact" means that it generates the same probability distribution as the CME. At each step, two random variables are sampled from a uniform density function U [0,1] in order to choose which reaction will occur and what length the time step will have. The chemical species are then upgraded according to the corresponding stoichiometric vector, and the process is repeated until the simulation time range is covered. The crucial point with SSA is that the time step can become very small if the number of molecules of some reactant species is large, and it is difficult to guarantee that the simulation proceeds one reaction at a time. SSA is, thus, a computationally demanding approach that avoids its applicability when (large) reaction networks are needed to model biologically realistic phenomena.
The point with it is the multiscale nature of the underlying problem. First, some reactions are much faster than others and rapidly reach a stable state. In the SSA, the dynamics of the system are driven by the slow reactions. In a deterministic setting, this issue is called "stiffness". Second, if some chemical species are present in relatively small quantities, then modeling by a discrete stochastic process is appropriate, whereas, if other species appear in large quantities, then stochastic or deterministic differential equations are better tools to use. In the SSA, instead, all the species are considered as discrete stochastic processes [10].
Burrage et al. [18] addressed the issue of multiscale modeling by extending the approach of Haseltine and Rawlings in [22] in classifying slow, intermediate and fast reactions to slow, intermediate, and fast regimes, resulting in a completely general, adaptive, and partitioning approach for the stochastic simulation of multiscaled systems involving chemical reactions.
Finally, the issue of non-homogeneity, typical of in vivo conditions, has been considered, among others, in the papers by Berry [23], Schnell and Turner [24], and Nicolau and Burrage [25]. All of these authors use Monte Carlo simulations on a two-dimensional lattice. In particular, Nicolau and Burrage [25] introduced the Anomalous Stochastic Simulation Algorithm (ASSA) that proved to be an effective and efficient algorithm to simulate the anomalous diffusion when the environment is not well-mixed, due to the obstacles present within the cell, such as fixed proteins or cytoskelatal structures.
However, the SSA has been successfully applied to simulate, for example, the regulation of lambda phage [26], circadian rhythms [27,28], and many other kinds of GRNs.
When noise effects are still important but continuity arguments hold, by the application of the Central Limit Theorem and by matching the first two moments of the CME, we can write an SODE that describes the continuous time evolution of the Markov process X(t) involved in the SSA [29]. This equation is sometimes called the Chemical Langevin Equation (CLE) and takes the form where W(t) = (W 1 (t), . . . , W N (t)) is an N dimensional vector in which its individual elements are independent Wiener processes and B(X(t)) is a N × N matrix satisfying where ν = [ν 1 , . . . , ν M ] is the N × M stoichiometric matrix. Equation (2) is an Itô SODE representing the continuous stochastic evolution equation describing intrinsic noise.
Finally, when large numbers of molecules are taken into account so that we can talk about concentrations, the deterministic continuous setting holds. This averaged behavior is described by an ODE of initial value type of the form This is the standard regime of chemical kinetics where the Law of Mass Action applies.

An Example of Genetic Regulation in Bacteriophage λ
A known and simple model organism is that of bacteriophage λ, which is a virus that infects the bacterium E. coli and can take two forms: lysis or lysogen. In the lysogen state, the virus replicates passively whenever E. coli replicates. Under appropriate conditions, a state change from lysogen to lysis can be induced. Many mathematical models have been developed to describe the genetic regulatory behavior of λ. These include both deterministic [30] and stochastic [26,31] kinetics models.
Most of these models consider the behavior of the λ right operator control system and the action of regulatory proteins, CI repressor and Cro. The right operator region O R consists of three binding sites which are situated between the promoters P RM and P I for genes CI and Cro. The dimeric forms of repressor and Cro bind to these binding sites to regulate transcription [32]. The lysogenic state is mainly regulated by the P RM promoter and the CI gene, while lysis is mainly regulated by the P R promoter of the Cro gene. Shea and Ackers [33] have given a detailed model based on 40  The network can be described by the following reactions: The rate constant corresponding to reaction j is denoted by c j . According to the method used in [17] for modeling first order, second order, and homodimeric reactions, we use variables X i (t), i = 1, . . . , 6, to characterize the molecular state of the transcriptional regulatory system at time t ≥ 0, where X i (t) = (M, D, RN A, DN A, DN A1, DN A2). From the Goutsias model (5), we obtain the following propensity functions a j (X)), j = 1, . . . , 10, and the stoichiometric vectors ν j , j = 1, . . . , 10 a 1 (X) = c 1 X 5 a 2 (X) = c 2 X 3 a 3 (X) = c 3 X 1 a 4 (X) = c 4 X 3 a 5 (X) = c 5 X 2 1 a 6 (X) = c 6 X 2 a 7 (X) = c 7 X 2 X 4 a 8 (X) = c 8 X 5 a 9 (X) = c 9 X 2 X 5 a 10 (X) = c 10 X 6 We initialize the system with m ≥ 1 monomers and d ≥ 1 dimers, and we assume g ≥ 1 DNA templates per cell. In this case, we have the initial state vector X 0 = (m, d, 0, g, 0, 0).
For simulating the solution of the Goutsias discrete model (5), we used an enhancement of the Gillespie SSA developed by the Kevin Burrage's team in the years 2004-2006 (K.B, Tianhai Tian, Manuel Barrio, André Leier, Tatiana Marquez-Lago, and the first author of this paper) using a "sorting" algorithm, from now on addressed as the "sorting SSA", that allows the speeding up of the performance of the algorithm when choosing what reaction number should take place at every time step and uses a special time stepsize for monitoring and storing the updated concentrations of species [34].
We recall that, in a stochastic setting, we can provide either single simulations of just one strong solution over a given path, or a number of independent simulations in order to collect information about mean behaviors. Since individual solutions are, however, representative of the dynamics of the processes being modeled, in Figure 1 that follow we will show single simulations of one strong solution.

Demographic Noise in Bacteriophage Infections: The Backward Gene Regulation Approach
Demographic noise in population dynamics and epidemic models has been studied by many authors with different approaches [5,14,[35][36][37][38]. As already mentioned in the previous sections, in all of these approaches the starting system is discrete. Now, suppose we have an ODE model of different kinds (biological, ecological, biochemical, etc.), and we wish to address the issue of how demographic noise affects the underlying interacting species, in order to capture the intrinsic variability neglected by the deterministic modeling. We consider the ODE model describing the dynamics of epidemics induced by bacteriophages in marine bacteria populations, introduced by Beretta and Kuang in [16].
The model equations are where S(t) is the fraction of total bacteria population susceptible of being infected at time t, I(t) is the fraction of infected bacteria at time t, P(t) are the (virulent) phage at time t, α > 0 is the bacteria intrinsic birth rate constant, C > 0 is the bacteria carrying capacity, K > 0 is the effective per bacteria contact rate constant with viruses, λ > 0 is the lysis death rate constant, and µ > 0 is the viruses death rate constant. Parameter b > 1 is the virus replication factor. The initial condition for system (9) is any point in By rescaling the variables on the carrying capacity C and introducing the dimensionless time τ = KCt, we obtain the dimensionless model equations where are the dimensionless parameters. In the simplified dimensionless form (10) the main mathematical features of model (9) were analyzed in [16]. The parameter b ∈ (l, +∞) was chosen as the parameter on which the dynamics of the infection depends, due to the fact that viral infection induces bacterial lysis, thus releasing, on average, b viruses for each infected bacterium into the marine environment. Beretta and Kuang proved that there exists a threshold value for b beyond which the (positive) endemic equilibrium bifurcates from the free disease one. In turn, for large values of b, the endemic equilibrium bifurcates toward a periodic solution. Model (9) was a seminal study for a quantity of other models, either deterministic or stochastic. Authors have also taken into account the influence of time delay, either discrete or distributed delay, resulting in delay differential equations (DDEs) or stochastic delay differential equations (SDDEs) [39][40][41][42][43][44][45][46]. We now view the Beretta and Kuang model (9) as the deterministic counterpart of the CLE associated with the discrete stochastic phenomenon, that is the bacteriophage infection on "units" of marine bacteria, under the effect of demographic noise. In other words, we consider the ODE model as the ∆t → 0 limit of a jump Markov process X(t) = (X 1 (t), X 2 (t), X 3 (t)), where X 1 (t) = S(t), X 2 (t) = I(t), X 3 (t) = P(t), satisfying system (9). The initial number of units of the single interacting species (susceptible and infected bacteria and phages) is assumed to be large.

Results
A thorough discussion of the analytical properties of the Beretta and Kuang ODEs model (9), including existence and uniqueness of the solution and the dependence of the dynamics of the model on the phage replication number b, can be found in [16]. For our purposes, using a standard ODEs method for non-stiff differential equations, the simulations of the solution in continuous deterministic setting result in Figure 2. Figure 2a shows that, if b is less than the critical value b * ≈ 16, only the boundary equilibrium E f = (C, 0, 0) exists. Susceptible bacteria reach their maximum value, while infected bacteria and phage vanish. Figure 2b,c show that for b * < b < b c , where a second critical value b c ≈ 95, a globally stable positive equilibrium E + = (S * , I * , P * ) exists, thus providing coexistence of bacteria and phage. As a last case, Figure 2d shows that when b > b c , where b c is the bifurcation value, then a positive equilibrium exists but is unstable. If we apply the backward method described in Section 2.3 to the ODEs model (9), we obtain the discrete stochastic system (12) and (13) modeling demographic noise in the phage-bacteria interaction and simulate the solutions via the sorting SSA. We expect that the more parameter b increases, the more the dynamics of the discrete system (12) and (13) is modified, accordingly to the continuous case. Figure 3 shows that the solutions fluctuate much more randomly with respect to the continuous model, particularly in the wide range of parameter b providing coexistence of the two species, i.e., b * < b < b c ( Figure 3b) and with respect to the phage population. Moreover, we note that the order of the solutions, whatever b is, diminish by one when we turn from "concentrations" (continuous deterministic case) to "integer number of individuals" (discrete stochastic case). The simulations of system (12) and (13) with the SSA solver is, however, slow; thus, the integration time is much shorter than the integration time reached for the fast modeling regime via the ODE system (9). So, if we wish to know the long-time behavior of the stochastic biological system, we have to consider the continuous stochastic regime represented by the Itô SODE (14). Using the Euler-Maruyama method for SODEs [41], we obtained Figure 4 confirming that even in the long-time behavior parameter b is crucial: for low values of b, the solution still goes to the corresponding equilibrium state (Figure 4a), although more randomly than in the deterministic setting. For higher values of b the effects of demographic noise, especially for the phage population, are much more evident: the positive equilibrium E + = (S * , I * , P * ) becomes unstable for values of parameter b much smaller than the deterministic bifurcation value b c = 94.451 (Figure 4b,c).  We also note that the order of the solutions, for all values of parameter b, increases by one with respect to the discrete stochastic case, thus regaining the same order of the original continuous, though deterministic, setting.

Discussion
In this article, we provided a method to offer insights into deterministic differential models when intrinsic noise is of interest. If the ordinary differential system is of biological or ecological nature, intrinsic noise is the so-called demographic noise. Our framework is that of chemical kinetics used to model intrinsic noise involved in molecular biology, in particular, in genetic regulation using SSAs. We first gave an example of application of the direct method for modeling demographic noise to a gene regulation network in bacteriophage λ, using a particular type of SSA developed by Burrage et al. in [17,34,47]. Next, we inferred the effects of demographic noise on an ODEs model of phage-bacteria interaction [16] going backward from the deterministic differential system to its discrete stochastic counterpart, thus giving insights to the effects of randomly fluctuating birth and death rates, immigration and/or emigration of the involved molecular species. We observe that our method is quite general and can be applied to ODE models of different nature, not only in life sciences. Unfortunately, at the moment its applicability is limited to ODEs where the underlying chemical reaction channels are either first order, heterodimeric (second and higher order), homodimeric, and Hill-type reactions for which a mathematical formulation can be provided, according to [17,19].
We are not able to compare our results with other authors' results since, to our knowledge, no equivalent methods starting from the deterministic setting can be found in the literature. Nevertheless, it would be interesting to compare the outcomes of our backward method with those obtained by applying the direct method to the same biological process in such a way to evaluate the robustness of the backward technique. We expect that the two approaches (direct and inverse) lead to very similar outcomes. Finally, a higher level of biological realism could be reached if, jointly with demographic noise, we consider time delay, for example, for modeling incubation time or latency period. In such a case, for the direct method, the Delay Stochastic Simulation Algorithm (DSSA) and the Delay Chemical Master Equation [15,34,48] would be the appropriate simulation tools to use. As with the backward (delay) method, there are good models to start with, for example, the Beretta and Kuang or the Beretta et al. [39,49] phage infection models with delay. All of these issues are being addressed in a new research article in preparation.