STOCHASTIC MODELS IN BIOLOGY AND THE INVARIANCE PROBLEM

. Invariance is a crucial property for many mathematical models of biological or biomedical systems, meaning that the solutions necessarily take values in a given range. This property reﬂects physical or biological constraints of the system and is independent of the model under consideration. While most classical deterministic models respect invariance, many recent stochastic extensions violate this fundamental property. Based on an invariance criterion for systems of stochastic diﬀerential equations we discuss several stochastic models exhibiting this behavior and propose classes of modiﬁed, admissible models as possible resolutions. Numerical simulations are presented to illustrate the model behavior.

1. Introduction. Random effects are encountered in many biological systems and it is often essential to take them into account in the modeling approach. Adding noise terms is one possibility to include stochasticity in a given deterministic ordinary differential equation (ODE) model and leads to stochastic differential equations (SDEs). In recent years, this problem has been discussed in a large variety of applications, and stochastic versions of well-known and well-understood deterministic models have been formulated. We mention several models within the biological and biomedical context: • Hodgkin-Huxley type models [6,10,11,18,27,28]; where f : R m → R m is a vector-valued function, and x 0 ∈ R m denotes the initial data. If parameters or coefficients in the equations are subject to random fluctuations, noise terms can be included and lead to systems of SDEs of the form dX = f (X)dt + g(X)dW, X = (X 1 , . . . , X m ), where X 0 ∈ R m , f is the deterministic part and g : R m → R m the stochastic perturbation. Moreover, W denotes a standard Wiener process and dW the corresponding Itô differential. Almost all cited stochastic models have been derived heuristically. Starting from a deterministic ODE model different approaches have been used to construct an SDE extension, and determining the "correct" stochastic perturbation is, in general, not an easy task. One method to derive SDEs is described in [1], Section 9.3., however, it can lead to models producing undesired values (e.g., see [33,37] and Section 3).
Since the paths of the Brownian motion are nowhere differentiable and of unbounded variation, (1) is formal notation and can be rigorously defined through the integral equation f (X(s))ds + t 0 g(X(s))dW (s).

STOCHASTIC MODELS AND THE INVARIANCE PROBLEM 2147
The last term cannot be explained as Riemann-Stieltjes integral. Among different possible interpretations Itô's and Stratonovich's definition of stochastic integrals are the most commonly used notions. However, the solutions generally depend on the choice of interpretation. They possibly exhibit different qualitative behavior that can lead to contradicting biological implications ( [3]). A long debate about this controversy has been going on, and depending on the specific application, different arguments have been given to support Itô's or Stratonovich's interpretation. Which one should be used is typically difficult to decide. For detailed discussions and possible resolutions we refer to [3,4,34,35] and the references therein. Now, assuming a stochastic model extension has been derived and the appropriate stochastic calculus been chosen, the next step is to analyze whether basic constraints of the underlying deterministic model, such as the existence and uniqueness of solutions and the invariance property, remain valid for the stochastic version. While for several of the cited SDE models the well-posedness was discussed, the invariance property was rarely addressed, even if the authors often recognized in simulations that solutions attain undesired values, which is a fundamental issue for the models' viability. We had already observed this modeling problem for stochastic Hodgkin-Huxley type models in our previous work [6], but it is, in fact, not rare. With few exceptions all cited stochastic models do not satisfy the invariance property and are, indeed, not realistic. This observation is in contrast to the fact that invariance criteria for systems of SDEs have been established. Some are formulated in an abstract mathematical framework, but others are explicit and can be easily used to validate stochastic models.
Most of the mentioned SDE models were studied by numerical experiments and, based on the simulations, biological implications were deduced. As we pointed out, in several cases undesired values of the solutions have been obtained. Typical resolution methods are, e.g., to include absorbing barriers or reflecting boundaries in the numerical code, where model variables exceeding the admissible range are set to the limiting value or reflected back into the allowed range (see [9], p. 2071, [15], p. 189, or [33], p. 454). However, since such methods modify the original model, this raises delicate questions that need to be discussed and analyzed: What is the relation between the simulations and the original model? Can we rely on the numerical simulations to validate the model behavior and to deduce biological interpretations?
Thus, several difficulties can arise when stochastic extensions of deterministic ODE models are developed: We need to determine the form of the stochastic perturbation and to decide which interpretation is appropriate for the model under consideration. Then, once the well-posedness is shown, we should analyze whether the invariance property of the underlying deterministic model is preserved, i.e., if solutions take values within the admissible range. If the model behavior is illustrated in numerical simulations, we should question whether they reflect the original model, in particular, if methods like absorbing boundaries were used.
Our aim is to provide an answer to the invariance problem for SDE models. To this end we recall a general invariance criterion for stochastic systems and show how it can be applied to classify admissible models. The structure of our article is as follows: In Section 2 we formulate explicit necessary and sufficient conditions for the invariance property for SDE systems similar to (1). They can be easily checked in applications and provide an important tool for model validation. Moreover, the criterion is valid for Itô's as well as Stratonovich's interpretation, i.e., with respect to invariance the qualitative behavior of solutions does not depend on the choice of interpretation. In the subsequent sections we discuss several previous SDE models and apply our invariance criterion for their validation. It allows to explicitly characterize the class of admissible stochastic perturbations preserving the invariance property of the underlying deterministic system. It turns out that nearly all cited stochastic models are, in fact, unrealistic since they do not respect invariance. We suggest modified SDE systems as possible resolutions and show the behavior of solutions in numerical simulations. Whether the choices of stochastic perturbations we provide are the correct and biologically relevant ones remains to be justified and analyzed in each particular situation.
2. Invariance criteria for stochastic differential equations: A reminder. We recall a general theorem for systems of SDEs that yields explicit necessary and sufficient conditions for the invariance of rectangular subsets. It was shown in [6] and extends a previous result by A. Milian ([21]). For the proof and further statements about stochastic invariance we refer to [6].
Let (Ω, F, P ) be a probability space with a right-continuous increasing family F = (F t ) t≥0 of sub-σ-fields of F each containing all sets of P -measure zero. We consider systems of Itô SDEs of the form where 1≤j≤r : [0, ∞[×R m → R m×r a Borel-measurable mapping into the set of all R m×r -matrices. Furthermore, W : [0, ∞[×Ω → R r denotes an r-dimensional F -adapted Wiener process, dW the corresponding Itô differential, the initial time t 0 is non-negative and X 0 ∈ R m is the given initial data.
Since the paths of the Brownian motion are nowhere differentiable and of unbounded variation, (2) is a formal notation. The initial value problem is rigorously defined through the integral equation where the last term denotes the Itô integral ( [8,22]). Besides Itô's notion Stratonovich's definition of stochastic integrals is frequently used. In applications it is generally difficult to decide which calculus is appropriate, and the qualitative behavior of solutions can depend on the chosen interpretation ( [3,4,34]). However, there is an explicit transformation formula (see [6,8,22]), and we will see that the invariance property for SDEs is, in fact, independent of the interpretation.
In the sequel, we denote by (f, g) stochastic initial value problems of the form (2). The function g represents the stochastic perturbation and f the deterministic part. In the particular case that g ≡ 0, we obtain an unperturbed deterministic ODE system, which we denote by (f, 0). Since our aim is not to establish the wellposedness, but to study the qualitative behavior of solutions, we assume that for any initial time t 0 ≥ 0 and initial data X 0 ∈ R m there exists a unique solution of the stochastic problem (f, g). Definition 2.1. We call the subset K ⊂ R m invariant for the stochastic system (f, g) if for every initial data X 0 ∈ K and initial time t 0 ≥ 0 the corresponding solution X(t), t ≥ t 0 , satisfies i.e., the solution almost surely attains values within the set K.
Theorem 2.2. Let I ⊂ {1, . . . , m} be a non-empty subset and a i , b i ∈ R ∪ {∞} such that b i > a i , i ∈ I. Then, the set for all t ≥ 0 and i ∈ I. This result applies independent of Itô's or Stratonovich's interpretation of SDEs.
One particular and important case for applications is the non-negativity of solutions, i.e., the invariance of the positive cone.
for all t ≥ 0 and i ∈ I. This result is valid independent of Itô's or Stratonovich's interpretation.
If we apply Theorem 2.2 to unperturbed deterministic systems, we recover the well-known tangential condition for systems of ODEs, which is necessary and sufficient for the flow invariance of subsets in R m (see [36] or [23]). In particular, for the deterministic system (f, 0) the set K is invariant if and only if the function f satisfies the conditions in Theorem 2.2.
A "sign condition" is not sufficient for g; if one model variable approaches the limiting values of the admissible range, the stochastic perturbations in the corresponding equation need to vanish. Otherwise, the random fluctuations can cause unrealistic values of the solution.
Example. We consider the stochastic system Corollary 1 implies that for arbitrary initial data X 0 , Y 0 , Z 0 ≥ 0 the solutions remain non-negative if and only if f and g satisfy for all X, Y, Z ≥ 0.
If we additionally require upper bounds, e.g., for the variable X, further conditions are imposed: By Theorem 2.2, for arbitrary initial data Y 0 , Z 0 ≥ 0 and X 0 ∈ [0, 1] the corresponding solutions Y and Z remain non-negative and X attains values within the interval [0, 1], if and only if the conditions on f and g above hold for all Y, Z ≥ 0 and X ∈ [0, 1] and, moreover, 3. Stochastic models for HIV-1 population dynamics. Our first example is the stochastic model for HIV-1 population dynamics proposed by H.C. Tuckwell and E. Le Corfec in [33]. Its solutions can attain unrealistic negative values. Using the stochastic invariance criterion we propose a modified, admissible model and present numerical simulations to illustrate its behavior. Finally, we discuss related SDE models exhibiting similar problems and indicate possible resolutions.

3.1.
A model for early HIV-1 population dynamics. A stochastic model for early HIV-1 population dynamics was proposed and numerically studied by H.C.
Tuckwell and E. Le Corfec in [33]. The model is formulated as system of SDEs with four components: The model variable T denotes the uninfected, activated CD4 + T cells, V the free virions, L the latently infected CD4 + T cells, which do not immediately participate in viral reproduction, and I the actively infected CD4 + T cells, which are involved in ongoing viral emission. The underlying deterministic model is a slight modification of the previous ODE models [24] and [25]. Noise terms were added in order to take stochastic effects into account which arise by virtue of nature in the infection process when virions bind to receptors and interact with uninfected cells. Moreover, randomness is also assumed in the transition of uninfected cells into latently infected, or actively infected, cells and described by the transition probabilities p and 1 − p, respectively.
The model [33] is given by the following system of SDEs where the W 1 , . . . , W 4 are independent standard Wiener processes, and dW i denotes the corresponding Itô differential. The constants a, c, k 1 , k 2 , α, γ, λ and µ are positive and p ∈ (0, 1). For a complete description of the model and parameter values we refer to [33]. The model variables T, L, I and V represent cell densities, which are necessarily non-negative. Non-negativity of solutions is preserved by the deterministic model, however, Corollary 1 implies that the stochastic extension does not possess this property. Indeed, solutions emanating from non-negative initial data take negative values with positive probability. Proposition 1. The stochastic model for HIV-1 population dynamics (3) has the following properties: (i) Solutions of the underlying deterministic system (i.e., if g 1 = · · · = g 4 ≡ 0) remain non-negative. (ii) For any constant k 1 = 0 the stochastic model does not preserve non-negativity, independent of Itô's or Stratonovich's interpretation of SDEs.
Proof. To shorten notations we introduce the functions We directly verify that which implies that solutions attain negative values with positive probability. The stochastic perturbations in the governing equations for L and I violate the necessary conditions for non-negativity, while the stochastic terms in the equations for T and V are admissible. However, such a partial result does not imply nonnegativity for the variables T and V , since the dynamics of V depends on I, and once I becomes negative it may cause negative values of V.
The simulations in Figure 1 illustrate such a situation. We plot here the real part of one sample path for the solutions, which is complex-valued due to the negativity of the term T V . Indeed, the stochastic perturbations are well-defined, only if T V is non-negative, which is not ensured by the SDE model (3).
The possibility of solutions to attain negative values was, in fact, observed by H.C. Tuckwell and E. Le Corfec (see [33], p. 454, Section 3): "In the simulations we found it expeditious to insert reflecting barriers at very small uninfected cells and virions numbers [. . . ]. This prevented the random fluctuation terms from taking these variables to unphysical negative values [. . . ]." However, the resulting numerical simulations do not correspond to the original model, and applying such resolution methods in the numerical schemes raises delicate questions: • Can we rely on the numerical simulations to validate the model behavior and to deduce biological interpretations? • What is the relation between the simulations and the original model? Could we formulate a modified stochastic model that reflects the absorbing barriers used in the simulations? In order to avoid such difficulties, in the following section we apply the invariance criterion to deduce a new class of admissible models. There are different options, and we consider one particular model that is as close as possible to the original model [33], but preserves the non-negativity of solutions.

3.2.
A class of admissible models. We suggest a class of admissible SDE models that slightly modifies the original system (3). Replacing the functions g 2 and g 3 in the equations for L and I we consider the system where we use the same notations as in (3), and h 1 and h 2 are two functions satisfying h 1 (0) = h 2 (0) = 0. By Theorem 2.2 the solutions of model (4) remain non-negative, independent of Itô's or Stratonovich's interpretation. One possible choice is where a and b are positive constants. The intensities of the stochastic perturbations in the governing equations for L and I then depend on L, and I respectively, and not only on the cell densities T and V. Of course, this particular choice needs to be justified and the model behavior should be analyzed to decide whether it is in agreement with experimental findings. In Figure 2 we show numerical simulations for the modified stochastic model (4) with h 1 and h 2 as in (5) (a = 0.5, b = 0.1). The other parameter values are taken from [33]. No resolution methods are required to ensure that solutions remain nonnegative. The plots are qualitatively similar to the ones for the original model in [33] (p. 455, Fig.1) where reflecting barriers were used. 3.3. Related stochastic models. The stochastic model for HIV-1 population dynamics [33] has been further developed and analyzed. We comment here only on few of these extensions and related works.
In [15], based on numerical simulations A. Kamina et al. compared the model [33] with another stochastic model for early HIV-1 population dynamics which does not involve SDEs. The authors were aware that solutions of the SDE model (3) attain negative values and use absorbing barriers for all model variables in their numerical scheme, i.e., " If any variable resulted in having a negative value after each iteration, then zero replaces this negative value " (see [15], p.189). The simulations in [15] resemble the behavior of the sample paths in [33]. However, the same questions arise regarding their reliability.
In [37] Y. Yuan and L.J. Allen proposed stochastic models for viral infection and immune responses in the early stages of infection. One of the considered classes is formulated as SDE system and closely related to the stochastic model (3). The underlying deterministic model is the same, except that immune response is included and actively infected and latently infected cells are not distinguished.
The model variables are the density of target cells T , the actively infected cells I and the virions V that satisfy the following system where W 1 , . . . , W 5 are independent standard Wiener processes and dW i denote the corresponding Itô differentials. The constants a, b, c, k, α, γ, ρ, λ, δ V , δ I , δ V and µ are positive. This model describes virion release through bursting, assuming that 2154 there is a coupling of the release of free virions and the burst of infected cells. An alternative release strategy is budding, where the death of infected cells and the release of virions are assumed to be independent processes. The stochastic model for budding in [37] is based on the same deterministic model, but the stochastic perturbations are slightly different, where all constants are positive and W 1 , . . . , W 4 are independent standard Wiener processes. Both stochastic models exhibit the same problems as system (3).

Proposition 2.
While the underlying deterministic model preserves non-negativity, for Itô's as well as Stratonovich's interpretation, the solutions of the stochastic models (6) and (7) attain negative values with positive probability.
Proof. The statements are a direct consequence of Corollary 1.
Based on the positivity criterion the SDE models can be modified to accomplish viability. One possibility is to omit the stochastic perturbations violating the conditions in Corollary 1 and to consider the following SDE model for bursting Similarly, viable stochastic models for budding can be obtained. However, as we pointed out previously, it has to be analyzed and justified whether such models resemble the behavior of the biological system and are in agreement with experimental data. We close this section with a remark about the SDE model for HIV internal dynamics studied in [7] which is closely related to the models we discussed. The underlying deterministic system is the same as in [37], but without the effect of immune response. The stochastic perturbations, however, are essentially different, where all constants are positive and W 1 and W 2 are independent standard Wiener processes. As we immediately observe the model preserves non-negativity. In [7] the qualitative behavior of solutions was studied, in particular, the non-negativity of solutions was shown and their asymptotic behavior analyzed.

Stochastic models for virus transmission and infectious diseases.
In this section we consider the stochastic model for virus transmission by A.J. Arenas et al. [2] and the model for the spreading of infectious diseases by Z. Huang et al. [13].

Models for virus transmission.
Stochastic models for the transmission of respiratory syncytial virus in the region of Valencia, Spain, were proposed and numerically studied by A.J. Arenas et al. in [2]. The underlying deterministic system is a classical SIRS model, where the variable S denotes the susceptibles, I the infected and R the recovered individuals, where the constants µ, ν and γ are positive and β is a continuous, positive function.
In order to take random fluctuations into account A.J. Arenas et al. suggested SDE extensions. Perturbations of the birth rate µ led to the system while perturbations of the transmission β rate gave rise to the system where the constants α, δ and δ 0 are positive and W denotes a standard scalar Wiener process. The first stochastic model does not preserve the invariance property of the underlying deterministic model, whose solutions are non-negative and bounded by a maximum value. The solutions I and R of the stochastic extension (9) are non-negative, but S attains negative values with positive probability. Furthermore, the model does not preserve upper bounds for the solutions.
All variables of the SDE model (10) are non-negative and the sum S + I + R is bounded by 1.

Proof.
We observe that f 1 , f 2 and f 3 satisfy for all S, I, R, t ≥ 0, which implies that the ODE system (8) preserves non-negativity. Moreover, the total population T := S + I + R satisfies the ODĖ Since f (1) = 0, Theorem 2.2 implies that T is bounded by one, which concludes the proof for the deterministic system. In model (9) the stochastic perturbations in the governing equations for R and I satisfy the conditions in Corollary 1, but the stochastic perturbation in the first equation does not. Consequently, R and I remain non-negative, but S can attain negative values. The conditions necessary for upper bounds in Theorem 2.2, however, are never satisfied, not even for S since there is no upper bound for R.
For system (10) all variables remain non-negative. Moreover, the total population T satisfies the same ODE as for the deterministic system, which implies that T is bounded by 1. Consequently, S, I and R are bounded by 1.
Using the stochastic invariance criterion we can modify system (9) to obtain a viable SDE model. One possibility is the following system which is close to the original model, but preserves non-negativity and upper bounds for the solutions, where the constants β, α, µ C , µ F , η, ρ and γ are positive (e.g., see [9]). Moreover, (V * , C * , F * ) denotes the strictly positive equilibrium point of the system. Z. Huang et al. take random fluctuations of the environment into account by perturbing the parameters β, µ C , µ F and ρ and derive the following stochastic extension of Marchuk's model where the constant σ is positive, and W denotes a standard Wiener process. They study the local and global stochastic stability of the system using tools from the theory of random dynamical system. The problem of non-negativity of solutions, however, is not addressed and is a fundamental property for the viability of the model. Proof. This is a direct consequence of Corollary 2.2.
A viable stochastic model that is close to the original one is the followinġ It takes the necessary conditions in Corollary 1 into account, and its solutions remain non-negative. 5. Stochastic Hodgkin-Huxley type models. Our next example is the stochastic model for cerebellar granule cell excitability proposed by A. Saarinen et al. in [28]. This model was subject of our previous work [6], where we showed that solutions attain undesired values and proposed a new class of admissible models. We give a short summary of the results and refer to [6] for further details and simulations.

5.1.
A stochastic model for cerebellar granule cell excitability. A stochastic model for cerebellar granule cell excitability was proposed and numerically studied in [27,28]. As many biophysical models of neurons it is based on the well-known deterministic Hodgkin-Huxley formalism [11], which qualitatively describes the conduction and excitation in nerves. Such models are commonly formulated as systems of deterministic ODEs. The behavior of neurons and voltage-dependent ion channels, however, is known to be stochastic in nature, which motivates the modeling approach in [28]. The mathematical model is given as system of SDEs for the dependent model variables x i , which represent the gating variables for the specific ion channels, the transmembrane potential V and the intracellular calcium concentration C, dx 9 = f 9 (V, C, x 9 )dt + σ 9 dW 9 , dV = F (t, V, x 1 , . . . , x 9 )dt, The reaction functions in the equations for the gating variables are and the rate functions for activation α i and inactivation β i are continuous and positive. The system of SDEs is interpreted in the sense of Itô, W i (t), t ≥ 0, denote standard scalar independent Wiener processes, dW i the corresponding Itô differentials, and the parameters σ i are positive and constant, i = 1, . . . , 9. For the concrete form of the interaction functions F and G and the complete description of the model we refer to [27,28]. This model extends a previous deterministic model for cerebellar granule cell excitability by adding the stochastic terms σ i dW i (t) in the governing equations for the gating variables x i , i = 1, . . . , 9. Ion channel stochasticity has been detected experimentally and is due to the thermal interaction of molecules constituting an ion channel. It can be observed as random opening and closing of an ion channel at an experimentally fixed membrane potential (see [28,10]).
The gating variables x 1 , . . . , x 9 describe the opening and closing rates of the specific ion channels and necessarily take values within the unit interval [0, 1]. While the unperturbed deterministic model certainly ensures this property, independent of the choice of parameters σ i serious difficulties arise for the stochastic model (11).
If σ i = 0 for some i ∈ {1, . . . , 9}, then the set K is not invariant for the stochastic model (11) and solutions attain values outside of the unit interval with positive probability. This is valid for Itô's as well as Stratonovich's interpretation.
Proof. We observe that It was observed in numerical simulations and pointed out in [28] that the solutions attain undesired values. Indeed, the invariance property cannot hold if the parameters σ i are taken to be constant. To solve this modeling problem was mentioned as a challenge for future work.
Similar difficulties occurred for the stochastic Hodgkin-Huxley model developed in [10], where the stochastic perturbations have the same structure (see [10], p. 2071). Loops were included in the numerical code to project undesired values of the gating variables back into the admissible range. However, as we already mentioned such resolution techniques modify the original model and it is questionable if the resulting simulations truly reflect the model behavior.

5.2.
A modified, admissible model. We obtain viable stochastic models if we replace the constants σ i by appropriate functions g i , that ensure the desired invariance of the unit interval. To be more precise, we propose models of the form dx 9 = f 9 (V, C, x 9 )dt + g 9 (t, V, C, x) dW 9 , where x = (x 1 , . . . , x 9 ), and the stochastic perturbations g i : [0, ∞[×R 11 → R satisfy g i (t, y) = 0 for y ∈ K such that y i ∈ {0, 1}, for all t ≥ 0 and i = 1, . . . , 9. For instance, an admissible choice is Proposition 6. The modified stochastic model (12) ensures that the gating variables x i take values within the interval [0, 1], for all i = 1, . . . , 9. This is valid for Itô's and for Stratonovich's interpretation of SDEs.
Proof. This is a direct consequence of Theorem 2.2.

6.
A stochastic JAK-STAT intracellular signaling pathway. In this section we discuss an example of different nature. The modeling problem and lack of invariance is not caused by adding noise terms as in the previous cases, but by an intermediate step.
The stochastic system was proposed by C. Surulescu and N. Surulescu in [29] and is based on the deterministic model for cellular signaling pathways developed by I. Swameye et al. in [31].
6.1. Underlying deterministic models. Complex intracellular signaling networks control processes like cell proliferation, differentiation and survival. A mathematical model for the core module of the signaling pathway of the Janus family of kinases (JAK)-signal transducer and activator of transcription (STAT) was proposed in [31]. It is formulated as the system of ODEs for the four different STAT5 populations,ẋ where EpoR is the amount of activated erythropoietin receptors, the cytoplasmic STAT5 is composed of the unphosphorylated STAT5 x 1 , the phosphorylated monomeric STAT5 x 2 and the phosphorylated dimeric STAT5 x 3 , while x 4 denotes the nuclear STAT5. The model describes how upon receptor recruitment monomeric STAT5 is tyrosine phosphorylized and converted into x 2 , how it dimerizes into x 3 and migrates to the nucleus, i.e., transforms into x 4 . The parameters k 1 , k 2 and k 3 are positive.
We could discuss the invariance of rectangular subsets for model (13), but for simplicity focus on the non-negativity of solutions. Proof. It can be easily verified that the interaction functions satisfy the conditions in Theorem 2.2.
The model (13) was further developed by I. Swameye et al. in [31] to allow for nucleocytoplasmic cycling, i.e., including the possibility that STAT5 is exported from the nucleus and reactivated in the cytoplasm. The extended model is of the formẋ where the constant k 4 is positive, x τ 3 (t) = x 3 (t−τ ) and the sojourn time τ of STAT5 in the nucleus is modeled by a fixed time delay τ > 0.
We do not focus on this model for the moment, but will get back to it later. In [29] C. Surulescu and N. Surulescu used a delay chain approach to replace the system of delay differential equations by a system of ODEs with an additional component. Furthermore, to include random effects in the process of receptor recruitment a noise term in the equation for the state variable x 1 was added. We analyze the resulting stochastic model in the following section.
Proposition 8. For the stochastic JAK-STAT model (15) the set is not invariant. More precisely, x 4 attains unrealistic negative values while all other variables remain non-negative.
Proof. The stochastic perturbation in the equation for x 1 is admissible and fulfills the conditions in Theorem 2.2. The conditions on the deterministic interaction functions hold for all variables, except for x 4 . Indeed, it would require that which is certainly not satisfied.
The simulations in Figure 3 illustrate one critical situation where the variable x 4 attains unrealistic negative values. We remark that the modeling problem is not caused by adding the stochastic perturbation, which actually is of admissible type if positivity is concerned. The problems arise from the system of delay differential equations (14) that was used to derive the underlying deterministic model in (15). 3. An admissible stochastic model. Using the invariance criterion the stochastic JAK-STAT model (15) can be modified such that it preserves the nonnegativity of solutions. One possibility is to consider the model where θ(t) is a positive continuous deterministic function and W a standard Wiener process. The simulations in Figure 4 illustrate the behavior of solutions, where we plot ten sample paths for the stochastic system (16).
7. Stochastic models for population growth. Finally, we discuss one of the simplest and most studied SDE models, namely, a scalar SDE describing population growth in random environments. We first address general stochastic growth models and then models for tumor growth. 7.1. Population growth in randomly fluctuating environments. Let the variable N denote the size of a population, e.g., a cell or biomass density. If we assume the population lives in an environment that is subject to random fluctuations, its growth can be modeled by the scalar SDE where f (N ) is the average per capita growth rate, and the effect of random fluctuations in the environment is described by the function σ(N ). Most commonly used growth functions include functions of logistic type where α > 0 denotes the growth rate and β > 0 the carrying capacity, or Gompertz growth functions of the form In both cases, for the unperturbed deterministic system the population size is non-negative and bounded by its maximum value β. Deterministic and stochastic growth models like (17) have been widely studied and, depending on the parameters involved, predictions about the species' persistence and extinction have been deduced. Most stochastic models consider constant noise terms σ(N ) = σ 0 . Even in this simple case, Itô's and Stratonovich's interpretation lead to different solutions and to different biological implications, see [3,4,20,34] and the references therein. A solution to this controversy was suggested in [3].
More general noise terms were considered by C.A. Braumann in [4], however, besides other hypotheses he assumes that the function σ satisfies He mentions that models that "may have some 'privileged' population sizes for which the noise intensity is zero (i.e., for which the growth rate is 'miraculously' unaffected by environmental fluctuations)" are not so realistic (see [4], p. 632). However, based on the invariance criterion, it might be more appropriate to require that the stochastic model does not only preserve positivity, but also the maximum bound for the population size. This is guaranteed if the function σ satisfies σ(0) = 0 and σ(β) = 0.
For instance, we could choose where γ is a positive constant. Such models have been mentioned by M. Turelli in [34] (p. 162 and p. 166), but, up to our knowledge, have not been further investigated.
7.2. Models for tumor growth and chemotherapeutic treatment. We now discuss a stochastic model for tumor growth, that falls into the general class of growth models (17). Different deterministic ODE models for the growth of cancer tumors and the effects of chemotherapy have been proposed and studied. Commonly used are Gompertz type growth equationsṗ for the number of tumor cells p, where µ and λ are positive constants; µ denotes the maximal number of tumor cells and λ the growth rate (see [32]). During chemotherapy the tumor growth kinetics is perturbed by the cytotoxic drug, which can be described by a loss term G : R → R of the from where the time-dependent function v is non-negative. It represents the dose of the drug v, and k 1 and k 2 are positive constants. Taking into account the effect of chemotherapy leads to the model In [19], H. Lisei and D. Julitz introduce a stochastic version of (18) and consider the SDE model where σ is a positive constant and W denotes a standard Wiener process. The existence, uniqueness and non-negativity of solutions for problem (19) can easily be established, since a change of variables transforms the SDE with multiplicative noise σpdW into a family of random ODEs. Arguments for the specific form of the stochastic perturbation have not been given and it may have been chosen, since it allows for explicit computations. Using the stochastic invariance criterion, we can construct more general SDE models that take random behavior in the drug therapy into account and preserve not only the non-negativity, but also upper bounds for the solutions.
Proof. By assumption the stochastic perturbations satisfy the conditions for invariance in Theorem 2.2. We observe that the deterministic part fulfills f (0) = 0 and f (µ) = −λG(v)µ ≤ 0, which implies the invariance of the interval [0, µ] for the variable p. Moreover, it can directly be seen that the conditions for the deterministic part in the second equation are fulfilled.
In Figure 5 we plot ten sample paths for the extended stochastic model (20), 7.3. Viable stochastic models for tumor growth. We shortly mention one viable stochastic model for tumor growth that respects positivity and the upper bound for the population density. This model was suggested and analyzed by G. Rosenkranz in [26] and has the form with deterministic growth function and stochastic perturbation where the positive constants α and γ are the relative rates of neoplastic cell production and destruction, and β > 0 denotes the upper limit for the tumor population X. The constants ω 2 and τ 2 measure the relative strength of the external and internal fluctuations, and W denotes a standard Wiener process.
The simulations in Figure 6 illustrate the behavior of solutions, where we plot ten sample paths of the process X for θ = γ = β = ω = 1, α = 10 −17 , τ = √ 1.5 × α and X(0) = 0.5. Related stochastic models for tumor growth were studied by R. Lefever and W. Horsthemke in [17]. The underlying deterministic model is the same as in [26], but the stochastic perturbations are of a different form. The authors analyzed the stability of two stochastic models and compared the behavior of solutions depending on Itô's and Stratonovich's interpretation. In particular, the models are of the form where the deterministic growth function is the same as in (21), and the constant σ is positive. We observe that solutions of the first SDE are non-negative, but do not ensure the maximum bound for the tumor population. However, independent of Itô's or Stratonovich's interpretation, the interval [0, β] is invariant for the second model.

7.4.
A stochastic immunogenic tumor model. Finally, we mention one stochastic model that extends the scalar tumor growth model and takes immune response into account. It was proposed and studied in numerical simulations by R. Horhat et al. in [12], and its behavior numerically analyzed in [14]. The underlying deterministic system is a slight modification of the model [16], which describes the growth of tumor cells subject to immune attack by effector cells. The stochastic model is formulated for the tumor cells X and effector cells Y and takes the form dX = (a 1 X(1 − a 2 X) − XY )dt + (c 11 (X − X * ) + c 12 (Y − Y * ))dW 1 , where a 1 , a 2 , b 1 , b 2 , b 3 and c ij , i, j = 1, 2, are positive constants, (X * , Y * ) denotes an equilibrium point of the unperturbed deterministic system, and W 1 and W 2 are standard scalar Wiener processes. The stochastic extension does not preserve the invariance property of the underlying deterministic model. Proof. This is a direct consequence of Theorem 2.2.

Remark 2.
For certain parameter ranges the deterministic model also preserves an upper bound for the effector cells. However, since the parameter values chosen in [16] lie outside of this range we formulated the invariance property in the previous proposition with respect to the set [0, 1 a2 ] × [0, ∞[. As in previous sections we could use the stochastic invariance criterion to suggest modified SDE models that preserve non-negativity and, if desired, also upper bounds for the tumor and effector cells. 8. Conclusion and perspectives. We formulated a general invariance criterion for systems of SDEs which is valid for Itô's as well as Stratonovich's interpretation. It provides explicit necessary and sufficient conditions that can easily be verified in applications, and hence, constitutes a powerful tool for model validation. When deriving stochastic versions of established deterministic models difficulties can arise, and we gave a variety of examples for stochastic models that produce unphysical values. The stochastic invariance criterion can directly be applied to derive viable stochastic models and yields resolutions for these modeling problems. We illustrated its use in many biological and biomedical applications, pointed out problems of previous stochastic models and proposed modified admissible models. Numerical simulations illustrate our results and observations. The problem of invariance for systems of SDEs can be brought into a more general framework, namely, the study of structure preserving SDEs: Assuming that a deterministic model possesses certain properties, e.g., invariant manifolds, first integrals, symmetries, flow invariance, etc., can we characterize the set of stochastic perturbations, seen as Itô or Stratonovich equations, that preserve these properties?
Such questions can be handled in many different ways. For variational structures, like Lagrangian or Hamiltonian structures, symmetries and first integrals, we refer to [5]. This article is based on the formalism of stochastic embeddings which is explicitly designed to analyze these properties.
As we already mentioned, most of the models based on SDEs have been studied by numerical experiments. With respect to invariance the validation of models by numerical methods is a difficult task, not only in the stochastic, but also in the deterministic case. Even if the underlying continuous model preserves invariance several numerical methods produce values lying outside of the admissible range. Consequently, a natural demand is to analyze whether numerical schemes preserve the invariance property and to develop and characterize such admissible numerical integrators. This will be subject of future work.