Intrinsic noise in systems with switching environments

We study individual-based dynamics in finite populations, subject to randomly switching environmental conditions. These are inspired by models in which genes transition between on and off states, regulating underlying protein dynamics. Similarly switches between environmental states are relevant in bacterial populations and in models of epidemic spread. Existing piecewise-deterministic Markov process (PDMP) approaches focus on the deterministic limit of the population dynamics while retaining the randomness of the switching. Here we go beyond this approximation and explicitly include effects of intrinsic stochasticity at the level of the linear-noise approximation. Specifically we derive the stationary distributions of a number of model systems, in good agreement with simulations. This improves existing approaches which are limited to the regimes of fast and slow switching.


I. INTRODUCTION
There is now a broad consensus that noise plays a crucial role in most dynamical systems in biology, chemistry, and in the social sciences. The theory with which to describe these stochastic processes is well established and has its roots in statistical physics. Modelling tools such as master equations, Fokker-Planck equations, and Langevin dynamics are standard and can be found in a number of textbooks [1][2][3]. Much of this work focuses on processes between discrete interacting individuals, which can be members of a population in epidemiology [4][5][6], atoms or molecules in chemical reaction systems [7,8], or proteins in the context of gene regulatory networks [9,10]. Many of these models are Markovian and their natural description is in terms of a master equation which describes the time evolution of the underlying probability distribution of microstates.
Solving the master equation analytically is often a difficult task: Only very simple linear model dynamics allow further treatment [1][2][3]11,12]. It is therefore common to employ approximation techniques, most notably ones built around the assumption that the size of the population is large but finite. The inverse system size (or its square root) then serves as a small parameter in which an expansion can be performed. The lowest order in this expansion corresponds to the limit of infinite systems and provides a deterministic description devoid of stochasticity. The second-order term in the expansion introduces some stochastic effects of noise in the population but approximates the individual-level dynamics by a simpler Gaussian process on a continuum domain [13]. These techniques have been very successful in capturing elements of noise-induced phenomena, for example, the weak selection of competing species [14][15][16][17][18], cyclic behavior, patterns, and waves [19][20][21]. The main techniques used to characterize these effects are system-size expansion methods, most notably the Kramers-Moyal and van Kampen expansions. The latter is also known as the linear-noise approximation (LNA) [1,2].
These methods are now used routinely for the analysis of individual-based models with intrinsic noise in the weaknoise limit. Most applications so far focus on systems in which the reaction rates are set by constant model parameters, and the only time dependence is in the evolution of the population of individuals itself. Recently, exceptions have gained attention [6,[22][23][24][25][26]. In these models, reaction rates vary continuously and deterministically in time, for example, to capture periodically varying infection rates to model seasonal variation in epidemic spreading. Crucially, no additional stochasticity is introduced in these dynamics by the environment, and the only discreteness in the dynamics is in the evolution of a finite population of individuals. Other authors have considered models with an environment which varies stochastically and continuously [27][28][29].
For many model systems it is more realistic to assume that model parameters switch between different discrete states. This includes phases of antibiotic treatment in the context of bacteria [30,31], genetic switches [9,[32][33][34][35][36][37][38], evolutionary game theory [39], and predator-prey models in switching environmental conditions [40,41]. Such models describe two types of discreteness: that of the state of the environment and that of the population of interacting individuals. In principle, the environmental switching can occur stochastically or follow a deterministic pattern (e.g., prescribed periods of antibiotic treatment, regularly interspersed with periods of no antibiosis).
In the present work we focus on stochastically switching environments. Assuming again a large but finite population, the demographic noise in the population can be approximated using the above expansion techniques. The noise relating to switches in the environmental state, however, cannot be dealt with in this way: There is no large parameter to expand in when the number of environmental states is small.
In models with switching environments the lowest-order expansion in the strength of the intrinsic noise leads to a "piecewise-deterministic Markov process" (PDMP) [42,43]. In this approximation the dynamics of the population of individuals is described by deterministic rate equations between stochastic switches of the environment. This approach neglects all intrinsic stochasticity from the reaction dynamics within the population-the population scale is taken to be infinite. The only type of randomness retained is that of the switching of the environmental states. The application of PDMPs has recently gained attention in the description of genetic networks [32,44,45].
The PDMP approximation has been surprisingly effective in modeling systems with very large populations [32,45,46]; however, it does not produce accurate results outside this limit. Recently, an alternative approach has incorporated some effects of demographic noise, but it is only valid if there is a very large separation between the time scales of environmental switching and that of the population dynamics [37,38]. Here, we develop the theory further and construct a systematic expansion in the noise strength about the PDMP.
The remainder of this paper is organized as follows: In Sec. II we give a more formal introduction to the problem using a relatively simple linear model. We show how the system-size expansion can be applied in the presence of switching environments and we analytically derive the resulting stationary distribution for the linear model. In Sec. III we construct a more general theory and describe how the method can be applied to a larger set of model dynamics. Section IV contains applications to a number of model systems with nonlinear reaction rates and/or other additional features. In Sec. V we summarize our findings and give an outlook on future work. The Appendix contains further details of the relevant calculations.

A. Model definition
We first focus on a simple example with linear reaction rates. We consider a population of individuals of type A, and we write n for the number of individuals in the population at any given time. Individuals can be created and they can decay, so the model describes a birth-death process. The death rate per individual is assumed to be a constant d. The creation rate is taken to depend on the state of the external environment; individuals are born with rate b σ , where σ represents the state of the environment. This can be summarized as follows: for times at which the environment is in state σ . The parameter has been introduced as per normal convention to set a typical scale of the population size [1,2]. If the environment were to be fixed at σ , then the number of individuals in the system would fluctuate around the value b σ /d in the long run. These fluctuations can be expected to be of order 1/2 and reflect the demographic stochasticity.
At any given time, the state of the full system is completely described by the state σ of the environment and the number of individuals in the population n. We restrict ourselves to cases in which the environment has two states, σ ∈ {0,1}. The switching between these states is assumed to be independent of the state of the population n and it occurs with constant rates. We write λ + for the rate of switching from state 0 to state 1 and λ − for the rate of switching from state 1 to state 0.
This stylized model can be interpreted in the context of genetic networks as follows [9]. The two states of the environment, often labeled as G 0 and G 1 , correspond to regimes in which a certain promoter-a region on DNA which initiates transcription-is either inactive (G 0 ) or activated (G 1 ). A protein (A) is produced with rates b 0 and b 1 in the corresponding state. Independently of the state of the gene, proteins degrade at a constant rate d. This simple model has been studied, for example, in Refs. [9,37,38,45,47], and it can be written in the form The reaction rates of this model are linear in n. It is possible to develop exact solutions to linear models of this type using a generating-function approach [11,12,[47][48][49][50]. However, such solutions are often limited to simple model systems and frequently they only provide limited insight into the actual physical dynamics. We use this linear model to introduce our approximation method. In later sections we will then apply this approach to cases with nonlinear reaction rates, where an exact solution is no longer feasible.

B. Simulation of results and general behavior
To illustrate the general behavior of the model, we first show the outcome of a set of characteristic simulations in Fig. 1. Figures 1(a), 1(c), and 1(e) depict individual simulation runs. In each of these panels a trajectory of the population density x(t) = n(t)/ is shown as a solid line, and the switching of the environmental states is indicated by the shading of the background. Figure 1(a) shows an example of relatively slow switching. In each environmental state, the population tends to a fixed point, φ * σ = b σ /d, specific to the environment. It then fluctuates about this fixed point. The stationary distribution [ Fig. 1(b)] is bimodal. As the switching rate is increased, Figs. 1(c) and 1(d), the stochastic dynamics spends more time in between the two fixed points and the bimodality of the stationary distribution is lost; we observe a nearly flat distribution between the two fixed points. At very fast switching, Figs. 1(e) and 1(f), the stationary distribution becomes unimodal, peaked at a value between the two fixed points. The system spends most of its time fluctuating about a point in the interior of phase space, away from either of the two fixed points.
We set out to characterize this behavior analytically, and our aim is to approximate the stationary outcome of the dynamics. We will develop a general approximation technique, which is applicable to a wider class of models, in the next section. Before we do this, it is useful to outline our approach and to describe the main steps of the analysis in the simpler model defined in Eqs. (2).

C. Master equation
We write n(t) for the number of individuals at time t and σ (t) for the state of the environment. These follow random-jump Markov processes [2]. We write P (n,σ,t) for the probability to find the system and environment in state (n,σ ) at time t. We will frequently suppress the explicit time dependence to keep the notation compact.   , and (f) on the right show the stationary distribution from simulations, averaged over multiple runs (histogram). We also plot the theoretical prediction for the stationary distribution * (x) obtained from Eqs. (16), shown as solid lines for the three scenarios. The trajectories and stationary distributions have been obtained by application of the Gillespie algorithm [51,52]. Model parameters are b 0 = 1/3, b 1 = 5/3, d = 1, and = 150.
These equations consist of two components. The operators L 0 and L 1 characterize the creation and removal of individuals assuming a fixed state of the environment. They are given by where E is the shift operator [1]: Ef (n) = f (n + 1). It is important to note that operators, such as E or L σ , act on everything that follows to their right throughout our paper, e.g., Enf (n) = (n + 1)f (n + 1). The latter two terms in the master equation describe the switching between the two environmental states. The master equation (3) describes the time evolution of the full process, and in our analysis, we wish to calculate the joint stationary distribution of the number of individuals n and environmental state σ , P * (n,σ ).

D. Approximation of the master equation
We now proceed to approximate the above master equation. To this end, it is useful to define the population density x(t) = n(t)/ . Assuming that the typical system size is large but finite, the operators, L 0 and L 1 , can be approximated by a Taylor expansion with respect to −1 , This is the Kramers-Moyal expansion [1,2] in terms of the variable x = n/ , while keeping the state of the environment, σ , discrete. We write L σ for the operators obtained from this expansion, retaining only leading and subleading orders in −1 . The state of the system is now expressed in terms of x and σ , and we will write (x,σ ) for the resulting probability density. The explicit time dependence is again suppressed in this notation.
Substituting this into the master equation (3) allows us to approximate the process by This expansion differs from the standard Kramers-Moyal expansion in that the environmental states are not included in the expansion. Equation (6) accordingly is not a standard Fokker-Planck equation; it retains the discrete switching terms, akin to terms in the master equation of a conventional telegraph process. Equations (6) describe a diffusion process with a Markovian switching in between two sets of drift and diffusion [53,54].
Since Eqs. (6) contain multiplicative noise, it is difficult to solve these equations directly. In the following sections, we propose an approximation method to analyze this dynamics. Our approach is similar to the conventional linear-noise approximation, and describes the effects of intrinsic noise to subleading order.

E. Leading-order approximation: Piecewise-deterministic Markov process
In the above expansion we have retained leading and subleading terms in −1 . To proceed, it is useful to first consider the leading-order terms only, i.e., to study the limit → ∞. We obtain where we have written φ = lim →∞ n/ to indicate that we have taken the limit of infinite populations. In this limit the L σ are Liouville operators describing deterministic flow of φ. The functional form of this flow at any one time is entirely determined by the state of the environment, σ . The process φ(t) is a piecewise deterministic Markov process [42,43,53,55], a random process composed of deterministic motion in between the discrete environmental transitions. The PDMP is a description of the system which accounts for the stochasticity of the switching only; the demographic noise, on the other hand, is neglected in the above truncation after the leading-order term.
The long-term behavior of the system can be determined from inspection of the Liouville operators. In each state σ the process tends towards a stable fixed point, φ * σ = b σ /d. In the following, we assume that b 0 < b 1 so φ * 0 < φ * 1 . The dynamics can be illustrated by the flow diagram in Fig. 2. When the environment is in state 0 the trajectory of the PDMP moves towards φ * 0 , and when the environment switches the direction is reversed towards φ * 1 . In the long run, the PDMP will always take values in the interval between the fixed points φ * 0 and φ * 1 . The stationary distribution of the PDMP, * (φ,σ ), can be found by setting the derivatives with respect to time in Eqs. (7) to zero, followed by integration and rearrangement. Further details are discussed in a more general setting in the next This result is consistent with those reported by other authors [45,55,56].

F. Comparison against simulations
In Fig. 3 we illustrate the behavior of the PDMP.  3(d), and 3(f) depict the corresponding stationary distributions of the PDMP, as obtained from Eqs. (8). For slow switching rates, the marginal stationary distribution * (φ) = * (φ,0) + * (φ,1) is bimodal, with singularities at the end points. In this regime, the PDMP typically spends enough time in each environment between switches to come close to the corresponding fixed point. For fast switching * (φ) is unimodal. In this situation, the system typically does not have sufficient time to reach the vicinity of the fixed points. An intermediate case is shown in Figs. 3(c) and 3(d). For this particular choice of parameters, the resulting stationary distribution is seen to be flat between the two fixed points, φ * 0 and φ * 1 . Comparison of the stationary distributions of the PDMP with those of the process in finite populations (Fig. 1) shows that the PDMP approximation recovers some of the qualitative features of the full system but not all. The transition from a bimodal to a unimodal shape is successfully reproduced. On the other hand, the singularities for slow switching rates seen in the stationary state of the PDMP are not observed in the full process. The PDMP is confined to the interval of the stationary distribution of the model with intrinsic noise includes concentrations below φ * 0 and above φ * 1 . These results stress the significance of intrinsic noise for the dynamics of the system. In order to approximate the stationary behavior to a better accuracy, it is necessary to include higher-order terms in the above Kramers-Moyal expansion. The corresponding formalism is well established for systems without random switches of the environment, and it takes the form of an expansion about the deterministic path of the infinite system [26]. In our case, the leading-order behavior (the PDMP) is a stochastic process itself due to the randomness of the environmental switching. The subleading description we discuss below is hence an expansion about this random process.

G. Subleading order: Linear-noise approximation
Before we present a more detailed account of the expansion to subleading order (Sec. III), we briefly outline the general idea. For a fixed realization of the environmental switching process, σ (t), we decompose the dynamics of the population as follows: This is along the lines of the system-size expansion in systems with time-dependent rates [26,57]. For a given path of the environment, σ (t), the trajectory φ(t) of the resulting PDMP is given by the solution oḟ where the birth rate b σ (t) at time t is determined by the state of the environment at that time, σ (t). The dynamics of Eq. (11) is the equivalent of the usual rate equations for systems without environmental switching.
To subleading order, systems of this form, but without environmental switching, are described by stochastic differential equations of the forṁ where we have suppressed the obvious time dependence of x on the right-hand side. The dependence of the amplitude w(x) on x indicates multiplicative noise, and η(t) is Gaussian white noise of unit amplitude, i.e., In a birth-death process with fixed birth rate b and death rate d The equivalent of Eq. (12) for the model with environmental switching is given bẏ Equation (13) is similar to the outcome of the procedure described in Refs. [26,57], where the system-size expansion was carried out for systems with periodically driven rates. The main difference is the fact that σ (t) is now a stochastic process itself, whereas the external driving was deterministic in Ref. [26]. In the limit → ∞ the noise η(t) does not contribute, and we recover the PDMP dynamics of Eq. (11). For the case of the linear model described above we have v σ (x) = b σ − dx and w σ (x) = b σ + dx. These are the drift term and noise amplitude one would obtain from a standard Kramers-Moyal expansion at fixed environmental state σ .
In the spirit of the usual linear-noise approximation (LNA) we next write x(t) = φ(t) + ξ (t)/ √ and obtaiṅ We have written v σ = dv σ (φ)/dφ, and we have suppressed the time dependence of φ on the right-hand side.
From the LNA it is possible to proceed and to approximate the stationary distribution of the process, * (ξ,φ,σ ) = * (ξ |φ,σ ) * (φ,σ ). The distribution * (φ,σ ) is the stationary outcome of the PDMP, and it can be computed exactly, see Eqs. (8) for the linear birth-death model. The general case is discussed in the following section.
In the stationary regime the distribution of ξ will, in principle, depend on the state σ of the environment and on the variable φ. In order to proceed we now assume that the dependence on σ can be neglected, so we write * (ξ |φ,σ ) ≈ * (ξ |φ). This is an approximation, but it turns out to work well for all models we have tested. Making this assumption we write * (ξ,φ,σ ) ≈ * (ξ |φ) * (φ,σ ).
The distribution * (ξ |φ) is Gaussian with mean zero, and with a variance which depends on φ and which can be obtained analytically as described in the next section. The resulting picture is illustrated in Fig. 4. A given realization of the environmental process generates a realization of the PDMP. The state of a finite population, subject to the same path of environmental states, will fluctuate about the PDMP trajectory, as indicated by the shading in Fig. 4. From these approximations the stationary distribution of x = φ + −1/2 ξ then can be estimated as Returning to Fig. 1, we compare the stationary distribution of the linear model as obtained from Eq. (16) against the results of numerical simulations. The results from the theory are shown as solid lines, and the simulation data as histograms. We find that the approximation reproduces the numerical results to a good accuracy.

A. Definition and master equation
The linear model discussed so far was deliberately simple and the main purpose of studying it was to develop a general intuition. In order to extend the method beyond the linear case, we now consider a more general model. This will introduce several new aspects to the problem.
As before, we restrict our discussion to the case of a single species and two environmental states. We write n for the number of individuals in the population, and σ ∈ {0,1} for the state of the environment. We use the notation λ + (n) for the rate with which the environment switches from state 0 to state 1, and λ − (n) for the rate of switches in the opposite direction. Unlike in the previous sections, we now allow for an explicit dependence of these rates on the state of the population, λ ± = λ ± (n).
We assume that there are M possible reactions in the population, labeled m = 1, . . . ,M. Each reaction m occurs with rate a m,σ (n) dependent on the current state of the environment and population. Any occurrence of a reaction of type m is taken to change the number of individuals in the population by S m . These are the underlying stoichiometric coefficients. The set of propensity functions, a m,σ (n), together with the stoichiometric coefficients completely define the dynamics of the population.
The operators L 0 and L 1 are given by Anticipating the system-size expansion, we write the transition rates in the form a m,σ (n) = r m,σ (x), where x = n/ and where r m,σ carries no explicit dependence. In slight abuse of notation we will write λ ± (x) instead of λ ± ( x). The master equation can be expressed in terms of (x,σ,t), the probability density of finding the random processes at x,σ at time t. The dynamics of the joint probability distribution can be approximated by where the operators L 0 and L 1 are given by As before, we have truncated the expansion of the operators in powers of −1 after the subleading terms. In explicit form we have where we have introduced Again, the process described by Eqs. (19) contains multiplicative noise and it is difficult to solve these equations in general.
In the next section, we begin by analyzing the dynamics in the limit of an infinite population.

B. Leading-order approximation: Piecewise-deterministic Markov process
We first analyze the process in the limit of infinite system size. Here we outline the main steps of the analysis and give the central results. The details of the calculation can be found in Appendix A. As before we write φ instead of x in the limit of an infinite system, and we find the PDMP dynamics This indicates a flow of the forṁ in between switches of the environment. These switches in turn occur with rates λ ± (φ). We now proceed by assuming that the deterministic flow in each of the environments is towards a unique fixed point, i.e., that there are points φ * σ for σ ∈ {0,1} such that v σ (φ * σ ) = 0.
For the time being we assume that there is only one such fixed point per state; for a more general case see Sec. IV C. We assume φ * 0 < φ * 1 without loss of generality. After a potential transient the PDMP will eventually be confined to the interval (φ * 0 ,φ * 1 ). Our assumption above (only one fixed point in each environmental state) implies that the v σ do not change sign on this interval, v 0 (φ) < 0 and v 1 (φ) > 0 for φ ∈ (φ * 0 ,φ * 1 ): the flow is always towards fixed point φ * 0 in state σ = 0, and towards φ * 1 in state σ = 1, as also illustrated in Fig. 2. To analyze the PDMP further, it is useful to introduce currents The physical interpretation of these currents is illustrated in Fig. 5. We focus on a domain (φ * 0 ,φ), where φ * 0 φ φ * 1 . The first term of the right-hand side of Eqs. (26) accounts for the probability flowing out of such a domain at location φ due to the Liouville flow in environmental state σ . The term containing the integral describes the net flow of probability out of the domain due to switching between the environmental states. Thus, the quantity J σ (φ + φ) − J σ (φ) represents the total amount of probability leaving the interval (φ,φ + φ) per unit time.

C. Subleading order: Linear-noise approximation
We now proceed by including contributions of intrinsic noise to subleading order. We focus on a time interval between switches of the environment, i.e., we assume σ is constant. During such time intervals the environmental noise has no effect, and the problem reduces to that of a conventional individual-based system with a fixed environment. Following established procedures we write n/ = φ(t) + −1/2 ξ (t). The deterministic dynamics is given byφ = v σ (φ). The outcome of a standard LNA can be expressed as a linear Langevin equation, for fluctuations about φ(t), where These relations describe the evolution of the population (within the LNA) between switches of the environment. When a change of the environment occurs the variable σ (t) changes at a discrete point in time, and the next such interval of constant σ begins. Within the LNA the evolution of the probability to observe state φ,ξ,σ is described by the following set of equations: While the expressions on the right-hand side look complicated, the different terms have a clear physical meaning. The ξ (φ,ξ,σ )/2 capture the evolution of ξ within the LNA of Eq. (31). Finally, the terms proportional to λ ± (φ) describe switching of the environment. Consistent with the expansion in the system-size we have replaced λ ± (x) by λ ± (φ), i.e., any dependence of the switching rates on the state of the population is taken to be on the state of the PDMP φ.

D. Stationary state within the linear-noise approximation
At stationarity the time derivatives on the left-hand side of Eqs. (33) vanish. Using asterisks as before to denote stationary distributions, and writing * (φ,ξ,σ ) = * (φ,σ ) * (ξ |φ,σ ), we find from summing the two equations (33) at stationarity. At this point we introduce a further approximation. We assume * (ξ |φ,0) ≈ * (ξ |φ,1) and simply write * (ξ |φ) for either of these. This will be justified below. Making this assumption leads to where we have used the relation v 0 (φ) * (φ,0) + v 1 (φ) * (φ,1) = 0, valid at stationarity, to show that the terms containing derivatives with respect to φ cancel out. Solving the stationary equation (35) is standard [3]; we find a Gaussian distribution, * (ξ |φ), with zero mean and with variance Using Eq. (30), this can be simplified to give We conclude this section with a brief comment on the approximation * (ξ |φ,σ ) ≈ * (ξ |φ). In the fast-switching limit the approximation is plausible, the system switches between environmental states too fast for the variable ξ to equilibrate to a distribution specific to the environmental state. For this case, the approximation reproduces the result of Ref. [37]. In the limit of slow switching, on the other hand, the system spends most of its time near the fixed points φ * 0 and φ * 1 . Our approach then also recovers appropriate distributions of ξ . For example, when φ ≈ φ * 0 , we have v 0 (φ = φ * 0 ) = 0 and Eq. (37) hence reduces to s 2 (φ * 0 ) = −[w 0 (φ * 0 )/(2v 0 (φ * 0 ))]. This is the stationary variance of the processξ = −v 0 (φ * 0 )ξ + w 0 (φ * 0 )η. A similar argument applies when φ ≈ φ * 1 . Outside the limits of fast and slow switching, the quality of the approximation can be verified numerically. In Fig. 6 we illustrate this for a model with nonlinear dynamics (discussed in more detail in the next section). To test the validity of the above assumption we have implemented the following measurement protocol: We first generate a combined realisation of the environmental process and the PDMP, i.e., a realization of the process defined by Eqs. (23). Subsequently, we generate a realization of the birth-death dynamics in a finite population for the same path of the environment. In the stationary state we then measure the variance of ξ conditioned on the state φ of the PDMP and on the environmental state σ . We denote this variance by s 2 σ (φ). This is then averaged over multiple realizations and compared against the approximation of Eq. (37). The data in Fig. 6 show that the approximation works well, with only slight deviations in the region away from the two fixed points φ * 0 and φ * 1 . For the linear model, numerical results and the theoretical prediction are virtually indistinguishable.

A. Nonlinear reactions rates
In order to demonstrate the generality of our approach, we now proceed to a model system with nonlinear reaction rates. The system is identical to the one described in Eq. (2), except that the last reaction (removal of individuals) is replaced by The above notation indicates that this last reaction occurs with rate dn(n − 1)/(2 ), where n is the number of individuals in the system. The reaction can be interpreted as a dimerization process, in which two particles of type A form a complex which is chemically inert and hence effectively removed. In addition to the switches between G 0 and G 1 , this system is described by two reactions, with stoichiometric coefficients and reaction rates We have written x = n/ as before. Using the notation of the preceding sections, we have the following drift and the diffusion terms In the limit → ∞, we obtain a PDMP with fixed points φ * σ = √ b σ /d for σ = 0,1. As before we assume b 0 < b 1 . The stationary distribution of the PDMP is found from Eq. (28), * (φ,0) = N for φ ∈ (φ * 0 ,φ * 1 ), and where N is the usual normalization constant.
From Eq. (37) finally we find s 2 (φ) = 3φ/4, i.e., * (ξ |φ) The stationary distribution * (x) is then obtained by numerically evaluating Eq. (16). Results are compared against simulation in Fig. 7, and we find convincing agreement between theoretical predictions and simulations. This confirms the validity of the assumptions and approximations made during the course of the analytical calculation.

B. System-dependent environmental transition rates
We now turn to another variation of the original linear model [Eqs. (2)]. For the dynamics within the population we use the same reactions and rates as in Eqs. (2), but we consider the 052119-9 case in which the rates with which the environment switches between states depends on the state of the population, i.e., where we choose the linear form λ ± (n) = α ± + β ± (n/ ). The coefficients α ± and β ± are constants, chosen such that λ ± remain non-negative.
The drift and diffusion terms of the dynamics within the population are as before. The stationary distribution of the PDMP in the limit of infinite populations is obtained from Eqs. (28) as where the exponents κ ± are given by The dynamics within the population is the same as in the linear model above, so evaluating Eq. (37) again leads to s 2 (φ) = φ. The theoretical estimate of the stationary distribution * (x), finally, is again found from numerical integration of Eq. (16). Figure 8 shows a sample path of the dynamics and the stationary distribution for the choice λ + (x) = λ − (x) = x (i.e., α ± = 0 and β ± = 1). In this case, transitions between states are more likely at higher concentrations of the system. Again, we find good agreement between the predicted distribution and the simulation results. It is important to note that the parameters we used in the figure are not special in any way: We have tested other choices of the parameter, and we find an agreement between simulations and theory of a similar quality.

C. Multiple fixed points
In the final example, we consider more complicated dynamics within the population such that there are multiple fixed points of the flows v σ (x). The switching between the two environmental states is taken to occur with constant rates λ ± . Specifically, the population dynamics is now modelled by the following reactions: with constant parameters c i,σ > 0. The corresponding stoichiometric coefficients for the four reactions are and the propensities r m,σ (x) read Additional details of the model and the numerical values of the parameters we used for our analysis can be found in Appendix B. Again, we first consider the PDMP, i.e., the limit → ∞. For the parameters chosen for our analysis one finds three fixed points of the dynamicsφ = v 0 (φ), at φ * ≈ 0.20,0.90, and 1.6. We label these φ (1) * ,φ (3) * , and φ (5) * . These fixed points are linearly stable, unstable, and stable, respectively. In the environmental state σ = 1 we have fixed points φ * ≈ 0.4,1.1, and 1.8, labeled φ (2) * ,φ (4) * ,and φ (6) * (again stable, unstable, and stable, respectively). This arrangement of fixed points is illustrated in Fig. 9. The dynamics of the PDMP depends on the initial condition. If started between the fixed points φ (1) * and φ (2) * , then the system will be confined between these two values and follow a dynamics similar to that of the system in Sec. IV A. Similarly, if the initial condition is between the two fixed points φ (5) * and φ (6) * , then the PDMP will operate in the interval between these two points. We will refer to these as the two "stable modes" of the PDMP dynamics. For other initial conditions, the PDMP will eventually reach one of these two stable modes as well, and which one this is will depend on the starting point and on the exact path the environment takes.
If the system is finite, then the dynamics are subject to intrinsic noise. Two representative trajectories are shown Fig. 10. Depending on initial conditions, the system either remains close to the interval between φ (1) * and φ (2) * or near the interval between φ (5) * and φ (6) * . Intrinsic noise will allow for excursions outside these intervals, similar to what was observed in Sec. IV A. The finite system can traverse the region between φ (2) * and φ (5) * , at least in principle, and move from one of the dynamic modes to the other. This is similar to escaping from a basin of attraction in systems with constant environment. We generally expect that the rate with which this happens is exponentially suppressed in the noise strength. We will assume that such events are sufficiently rare so they can be ignored for the purposes of our analysis.
The general approach we have developed can then be applied to the system in either of the two dynamic modes. Due to the complexity of the flow fields v σ (x), the relevant equations have to be evaluated numerically. In Fig. 10 we show the resulting predictions for the stationary distribution in either mode. As in the previous examples, comparison with data from numerical simulations shows very good agreement.

V. CONCLUSION
In summary, we have constructed a systematic approach with which to investigate the effects of demographic noise in systems subject to sudden random switches of reaction rates. We have focused on relatively simple birth-death processes of one single species and in which birth and death rates depend on both the state of the population and on the state of an external environment. The states of the environment follow a telegraph process, with transition rates which may depend on the state of the population. Previously existing approaches either disregard intrinsic fluctuations and only account for environmental noise or they focus on cases in which there is a clear separation of time scales between the population dynamics and the dynamics of the environment. The former approach in particular leads to the well-established picture of so-called piecewise deterministic Markov processes. Our work systematically improves on this view; we take into account demographic fluctuations to leading order and carry out a system-size expansion and linear-noise approximation about the PDMP dynamics.
Using the linear-noise approximation, and retaining the discreteness of the environmental process, we then approximate the resulting stationary distribution of the population dynamics. We have tested the resulting theory on a number of model systems, both with linear and nonlinear reaction rates, situations in which environmental switching depends on the state of the population, and covering systems with single dynamic modes and cases where there are multiple attractors. In all cases we have tested, the approximation is in good agreement with results from simulations. In particular, no separation of time scales is required.
The technique we provide makes our understanding of processes involving both intrinsic and extrinsic noise more complete. While the existing PDMP description has been shown to be successful in many instances, it disregards intrinsic noise entirely. We are now in a position to describe the effects of demographic noise in systems with switching environments in the spirit of the van Kampen expansion. This allows us to investigate models subject to combinations of intrinsic and extrinsic noise and, in particular, systems in which some degrees of freedom can be treated within a linear-noise approximation, while other variables remain fundamentally discrete. Such systems can be of relevance, for example, in the context of genetic switches, bacterial populations subject to varying external conditions, or to predator-prey dynamics. In order to make the method more applicable, several extensions can be considered in future work. This includes the case of multiple environmental states, systems with more than one species, or, indeed, environmental dynamics beyond simple telegraph processes. Work along these lines is currently in progress.
Equation (A3) and its analog for ρ 1 (φ) can directly be integrated and we find * (φ,0) where and where N 0 and N 1 are normalization constants. Using again the relation v 0 (φ) * (φ,0) + v 1 (φ) * (φ,1) = 0, derived above, we conclude N 0 = −N 1 ≡ −N , and so we have Recalling that v 0 < 0 and v 1 > 0 throughout the domain of the PDMP, N is a positive constant, to be determined from the normalization condition The accuracy of our approach can be characterized by computing the Kullback-Leibler divergence [58] between the true stationary distribution and our approximation. Figure 11 system size Ω shows how the Kullback-Leibler divergence changes with system size for the linear model. Each line represents one of the parameter regimes displayed in Fig. 1. Birth rates, death rates, and switching rates have been chosen such that is the average population of the system. The divergence decreases as the system size is increased, with relatively small system sizes having appreciably small divergences. Similar results are observed for the nonlinear and feedback models.