STOCHASTIC DYNAMICS AND SURVIVAL ANALYSIS OF A CELL POPULATION MODEL WITH RANDOM PERTURBATIONS

. We consider a model based on the logistic equation and linear ki-netics to study the eﬀect of toxicants with various initial concentrations on a cell population. To account for parameter uncertainties, in our model the coeﬃcients of the linear and the quadratic terms of the logistic equation are aﬀected by noise. We show that the stochastic model has a unique positive solution and we ﬁnd conditions for extinction and persistence of the cell popu- lation. In case of persistence we ﬁnd the stationary distribution. The analytical results are conﬁrmed by Monte Carlo simulations.


1.
Introduction. Cell-based in vitro assays [27] are efficient methods to study the effect of industrial chemicals on environment or human health. Our work is based on the cytotoxicity profiling project carried by Alberta Centre for Toxicology in which initially 63 chemicals were investigated using the xCELLigence Real-Time Cell Analysis High Troughput (RTCA HT) Assay [26]. We consider a mathematical model represented by stochastic differential equations to study cytotoxicity, i.e. the effect of toxicants on human cells, such as the killing of cells or cellular pathological changes.
The cells were seeded into wells of micro-electronic plates (E-Plates), and the test substances with 11 concentrations (1:3 serial dilution from the stock solution) were dissolved in the cell culture medium [20]. The microelectrode electronic impedance value was converted by a software to Cell Index (n), which closely reflects not only cell growth and cell death, but also cell morphology. The time-dependent concentration response curves (TCRCs) for each test substance in each cell line were generated [26] and based on these curves the toxicants in the present study were divided in 10 groups [30]. In Fig. 1 we display the TCRCs for the toxicant monastrol.
The success of clustering and classification methods depends on providing TCRCs that illustrates the cell population evolution from persistence to extinction. In [1] we consider a model represented by a system of ordinary differential equations to determine an appropriate range for the initial concentration of the toxicant. The model's parameters were estimated based on the data included in the TCRCs [1].  Let n(t) be the cell index, which closely reflects the cell population, C o (t) be the concentration of internal toxicants per cell, and C e (t) be the concentration of toxicants outside the cells at time t. We suppose that the toxicants do not exist in the cells before experiments, so C o (0) = 0, and that C e (0) is equal to the concentration of toxicant used in the experiments. We assume that the death rate of cells is linearly dependent on the concentration C o of internal toxicants and we consider linear kinetic, so we get the following deterministic model [1]: dC e (t) dt = λ 2 2 C o (t)n(t) − η 2 2 C e (t)n(t) Here β > 0 denotes the cell growth rate, γ = β K , where K > 0 is the capacity volume, α > 0 is the cell death rate, λ 2 1 represents the uptake rate of the toxicant from environment, η 2 1 is the toxicant input rate to the environment, λ 2 2 is the toxicant uptake rate from cells, and η 2 2 represents the losses rate of toxicants absorbed by cells.
The deterministic model (1)-(3) is a special case of the class of models proposed in [5], and it is related to the models considered in [7,11,15]. However, since we consider an acute dose of toxicant instead of a chronic one, the analysis of the survival/death of the cell population is different from the one done in the previously mentioned papers.
We have noticed that, for the toxicants considered here, the estimated values of the parameters η 1 , η 2 , λ 1 , and λ 2 verify η 2 1 η 2 2 − λ 2 1 λ 2 2 > 0 [1]. In this case we have 0 < C e (t) ≤ C e (0), 0 ≤ C o (t) ≤ , and n(t) > 0, for all t ≥ 0. (see Lemma 3.1 in [1]). Moreover from Theorem 3.2 in [1] we know that lim t→∞ C e (t) exists and its value determines the asymptotic behavior of the system: then |n| 1 = ∞ 0 n(t)dt < ∞ and the population goes to local extinction: In practice we usually estimate a parameter by an average value plus an error term. To keep the stochastic model as simple as possible, we ignore the relationship between the parameters β and γ, and we replace them by the random variables By the central limit theorem, the error terms may be approximated by a normal distribution with zero mean. Thus we replace equation (1) by a stochastic differential equation and, together with equations (2) and (3), we get the stochastic model Here σ i ≥ 0, i = 1, 2 are the noise intensities. (Ω, F, {F t } t≥0 , P) is a complete probability space with an increasing, right continuous filtration {F t } t≥0 such that F 0 contains all P-null sets, and B i , i = 1, 2 are independent standard Brownian motions defined on the above probability space. Several versions of a stochastic logistic equation similar with (5) were considered in [18], [19], [8], [9], [10] and [21]. The system of stochastic differential equations (5)-(7) is closely related with the stochastic models in a polluted environment considered in [15], [16], and [24]. However, for the models considered in these papers, instead of the equations (6) and (7), C o (t) and C e (t) obey two linear equations without any terms involving n(t). Moreover, instead of a combination of linear and quadratic terms as in (5), in [15] only a linear stochastic term is considered, and in [16] two stochastic competitive models are considered including exclusively either linear stochastic terms or quadratic stochastic terms.
In this paper we extend the methods applied in [15] and [16] to find conditions for extinction, weakly persistence, and weakly stochastically permanence for the model (5)- (7). In addition to this we focus on the ergodic properties when the cell population is strongly persistent. The main contribution of this paper is the proof that n(t) converges weakly to the unique stationary distribution. If only one of the noise variances σ 2 1 , σ 2 2 is non-zero, we also determine the density of the stationary distribution. For the study of the ergodic properties we apply techniques used for stochastic epidemic models in [4], [28], [29] and [23], and for a stochastic population model with partial pollution tolerance in a polluted environment in [25].
In the next section we prove that there is a unique non-negative solution of system (5)-(7) for any non-negative initial value. In section 3 we investigate the asymptotic behavior, and in section 4 we study the weak convergence of n(t) to the unique stationary distribution using Lyapunov functions. Numerical simulations that illustrate our results are presented in section 5. The last section of the paper contains a short summary and conclusions.
2. Existence and uniqueness of a positive solution. We have to show that system (5)-(7) has a unique global positive solution in order for the stochastic model to be appropriate. Let R + = {x ∈ R : x ≥ 0}, and R * + = {x ∈ R : x > 0}.
Since equations (6) and (7) are linear in C o and C e we have Let's define the differential operator L associated with the system (5)-(7) by where For any given initial value x(0) ∈ D the system (5)-(7) has a unique global positive solution almost sure (a.s.), i.e. P{x(t) ∈ D, t ≥ 0} = 1.
Proof. The proof is similar with the proof of theorem 3.1 in [29]. Since the coefficients are locally Lipschitz continuous functions, there exists a unique solution on [0, τ e ), where τ e is the explosion time ( [3]). To prove that the solution is in D and τ e = ∞ we define the stopping time where m > m 0 and m 0 > 0 is a positive integer sufficiently large such that n(0) ∈ [1/m 0 , m 0 ], 0 ≤ C o (0) ≤ m 0 , and C e (0) ∈ [1/m 0 , m 0 ]. Here we set inf ∅ = ∞.
Obviously {τ m } is increasing and let τ ∞ = lim n→∞ τ m , where 0 ≤ τ ∞ ≤ τ e a.s.. From formula (8) it is easy to see that C o (t) ≥ 0 for any t < τ ∞ . We show that τ ∞ = ∞ a.s., so τ e = ∞ a.s. and the solution is in D for any t ≥ 0 a.s. Assume that there exists T > 0, and > 0 such that P (τ ∞ ≤ T ) > . Thus there exists an integer We define the C 3 − function V : D → R * + as follows We get Omitting some of the negative terms, for any x ∈ D we have Since f is continuous on (0, ∞) and lim n→∞ f (n) = −∞ it can easily be shown that Using Itô's formula (10) forṼ and taking expectation we have for any m ≥ m 1 : Notice that for any ω ∈ Θ m , m ≥ m 1 we have V (x(τ m , ω)) ≥ b m = min{V (y)|y = (y 1 , y 2 , y 3 ) has the components y 1 or y 3 equal with m −1 or m, or y 2 = m}. Hence Thus we have proved by contradiction that τ ∞ = ∞.
Here we focus on the case when n(0) > 0, we have only an acute dose of toxicant C e (0) > 0, C o (0) = 0, and the external concentration of toxicant C e (t) is never larger than C e (0). For this we have to impose some conditions on the parameters. Similarly with the deterministic case we obtain the following results (for completion the proofs are included in Appendix A and Appendix B).
Definition 3.1. The population n(t) is said to go to extinction a.s. if lim t→∞ n(t) = 0 a.s..
< 0 a.s. then the population n(t) goes exponentially to extinction a.s.
Proof. The proof is similar with the proof of Theorem 6 in [16]. We start with some preliminary results. By Itô's formula in (5) This means that we have Notice that the quadratic variation [17] of Now we do the proof for part a. Using the exponential martingale inequality (Theorem 7.4 [17]) and Borel-Cantelli lemma ( [22], pp. 102), and proceeding as in the proof of Theorem 6 in [16] we can show that for almost all ω there exists a random integer n 0 = n 0 (ω) such that for all n ≥ n 0 we have Hence, for all n ≥ n 0 and all 0 ≤ t ≤ n we have Substituting the above inequality in (12) we get for all n ≥ n 0 , and any 0 < n − 1 ≤ t ≤ n. Since lim t→∞ Next we prove part b. Suppose that P (Ω) > 0 where Ω = {lim sup t→∞ n(t) ≤ 0}. From Theorem 2.1 we know that n(t) > 0 , t ≥ 0 a.s., so P (Ω 1 ) > 0 where Ω 1 = {lim t→∞ n(t) = 0}, and Ω 1 ⊆ Ω. Thus, for any ω ∈ Ω 1 we have Moreover, from the law of large numbers for local martingales (Theorem 3.4 in [17]) there exists a set Ω 2 ⊆ Ω 1 with P (Ω 2 ) > 0 such that for any ω ∈ Ω 2 we have From (12) we get: > 0 a.s., we have a contradiction with (13), so lim sup t→∞ n(t) > 0 a.s.
We have the following result regarding the expectation of n(t).
Proof. Using Itô's formula in (5) we get: Let for any m > m 0 , where m 0 > 0 was defined in the proof of Theorem 2.1. Obviously η m ≥ τ m , m > m 0 , where τ m is given in (11). In Theorem 2.1 we have proved that lim m→∞ τ m = ∞ a.s., so we also have lim m→∞ η m = ∞ a.s.. Taking expectation in (14) we get: Thus, there exists a constant K 1 > 0 such that sup t≥0 E[n(t)] ≤ K 1 .
Proof. For any > 0, set c 1 ( ) = K 1 / , where the constant K 1 > 0 is given in the previous lemma. From the Markov's inequality [22] we obtain Hence, from Lemma 3.6 we get By Itô's formula in (5) we get for any real constant c: , and we get: Taking expectation in (16) and using Lemma 3.6 we get: where η m was defined in (15). Letting m → ∞ we get For any > 0 set c 2 ( ) = /M 2 . From Markov's inequality we have Thus lim inf t→∞ P(n(t) ≥ c 2 ( )) ≥ 1 − , and this inequality and Corollary 1 implies that n(t) is stochastically permanent.
4. Stationary distributions. The deterministic system (1)-(3) has a maximum capacity equilibrium point (K, 0, 0) , where K is the capacity volume ( [1]). For the stochastic system (5)- (7), (K, 0, 0) is not a fixed point, and, when the cell population is persistent, we no longer have lim t→∞ n(t) = K. In this section we study the asymptotic behavior of n(t) when lim t→∞ C o (t) = 0 a.s.. For stochastic differential equations, invariant and stationary distributions play the same role as fixed points for deterministic differential equations. In general, let X(t) be the temporally homogeneous Markov process in E ⊆ R l representing the solution of the stochastic differential equation where B r (t), r = 1, . . . , d are standard Brownian motions. We define the operator L associated with equation (17): Let P (t, x, ·) denote the probability measure induced by X(t) with initial value X(0) = x ∈ E: P (t, x, A) = P (X(t) ∈ A|X(0) = x), A ∈ B(E), where B(E) is the σ−algebra of all the Borel sets A ⊆ E. Definition 4.2. The Markov process X(t) is stable in distribution if the transition distribution P (t, x, ·) converges weakly to some probability measure µ(·) for any x ∈ E.
It is clear that the stability in distribution implies the existence of a unique stationary measure, but the converse is not always true [2]. We have the following result (see lemma 2.2 in [29] and the references therein). Lemma 4.3. Suppose that there exists a bounded domain U ⊆ E with regular boundary, and a non-negative C 2 −function V such that A(x) = (A i,j (x)) 1≤i,j≤l is uniformly elliptical in U and for any x ∈ E \ U we have LV (x) ≤ −C, for some C > 0. Then the Markov process X(t) has a unique stationary distribution µ(·) with density in E such that for any Borel set B ⊆ E lim t→∞ P (t, x, B) = µ(B) for all x ∈ E and f being a function integrable with respect to the probability measure µ.
c. There exists a constant C 1 > 0 such that sup t≥0 E[X(t)] ≤ C 1 and, for any Proof. The proofs for a. and b. can be done similarly with the proof of Theorem 2.1, using the C 2 -function V : The proof of c. is analogous with the proof of Lemma 3.6.
Let P X (t, x, ·) denote the probability measure induced by X(t) with initial value X(0) = x ∈ R * + , t ≥ 0. In the following theorem, using Lemma 4.3, we show that the Markov process X(t) is stable in distribution.
Theorem 4.5. If σ 2 1 < 2β then the Markov process X(t) has a unique stationary distribution µ 1 (·) with density in R * + such that for any Borel set B ⊆ R * for all x ∈ R * + and f being a function integrable with respect to the probability measure µ 1 .
Proof. We consider the C 2 -function V : Since LV (·) is a continuous function on R * + and LV (0) = Then U is a bounded domain, and LV (x) < −C for any x ∈ R * + \ U , where C > 0 is the minimum between C 1 and C 2 . Notice that A(x) = σ 2 1 x 2 + σ 2 2 x 4 is uniformly elliptical on U , so the assumptions of Lemma 4.3 are met. Therefore, the Markov process X(t) has a unique stationary distribution µ 1 (·) and it is ergodic.
Proof. We follow the same idea as in the proof of Theorem 2.4 in [28]. From theorem 4.5 we know that X(t) w → t→∞ µ 1 , where µ 1 is a probability measure on R * + . By the Continuous Mapping Theorem [22], Y (t) = 1/X(t) also converges weakly to a probability measure ν 1 on R * + , the reciprocal of µ 1 . We will show that Indeed, if we denote ξ(t) = N (t) − Y (t), then ξ(0) = 0 and from equations (20) and (21) we get The solution of the previous linear equation is given by (see chapter 3, [17]) Obviously ξ(t) ≥ 0, t ≥ 0, a.s., and this means that we have Y (t) ≤ N (t) for any t ≥ 0 a.s. Similarly, using equations (21) and (22), we can show that Y (t) ≤ Y (t) for any t ≥ 0 a.s.. Secondly we show that for any 0 < < From equations (20) and (22) we get
a. If σ 1 = 0 the equation (21) become Let define It can be easily shown that where G 1 is given in (27). So, by Theorem 1.16 in [13], Y (t) is ergodic and with respect to the Lebesgue measure its stationary measure ν 1 has density Thus, by Theorem 4.5, X(t) = 1/Y (t) is ergodic and its stationary measure µ 1 is the reciprocal of the measure ν 1 , so with respect to the Lebesgue measure has density p(x) = p 1 (1/x)/x 2 given in equation (26). Notice that we also have a.s..
On the other hand, if σ 2 1 < β, η 2 1 η 2 2 − λ 2 1 λ 2 2 > 0, and lim inf t→∞ n(t) > 0 a.s., then according to Theorem 4.7 n(t) w → t→∞ µ 1 . Indeed, since lim inf t→∞ n(t) > 0 a.s., then ∞ 0 n(t)dt = ∞ a.s., and from the proof of Theorem 2.3 we know that lim t→∞ C o (t) = lim t→∞ C e (t) = 0, so the assumptions of Theorem 4.7 are satisfied. 5. Numerical simulations. First we illustrate numerically the results obtained in section 3 regarding survival analysis. We consider a cell population exposed to the toxicant monastrol as in the experiments described in [1] . The parameters' values for this toxicant are estimated in [1]: β = 0.074, K = 18.17, η 1 = 0.209, λ 1 = 0.177, λ 2 = 0.204, η 2 = 0.5, and α = 0.016. Notice that for this toxicant we have η 2 1 η 2 2 − λ 2 1 λ 2 2 > 0. We solve numerically the system (5)-(7) using an order 2 strong Taylor numerical scheme [12].  One of the applications of the mathematical model is for finding the threshold value for C e (0) at which the population becomes extinct. This value depends on the initial value n(0) and for the deterministic model (1)-(3) can be found numerically (see also Fig. 3 in [1]). From Theorem 3.5 we can see that large values of the noise variance σ 2 1 result in population extinction, so we expect that the presence of noise will lower the values of the threshold.
We illustrate this for the model with initial values n(0) = 2.  found numerically, and it is approximately C det e (0) = 382.2. We can see in Fig. 2 (a) that for the stochastic model with σ 1 = 0.01, σ 2 = 0 and initial value C e (0) = 380 the population goes to extinction, while in the deterministic case (σ 1 = σ 2 = 0) the population is persistent for these initial values. According to Fig. 2 (a), in the stochastic case the threshold value for this simulation is C stoch e (0) ∈ (375, 380). Similar results are obtained for the stochastic model with σ 1 = 0, σ 2 = 0.001 and are presented in Fig. 3 (a). For this simulation the threshold value in the stochastic case is C stoch e (0) ∈ (379, 380). Notice also that the results displayed in Figs Fig. 5) and the graph of the corresponding density given in Corollary 2 (a). In Fig. 6 we do similar plots for σ 2 = 0, several values of σ 1 , and the density of the corresponding Gamma distribution in Corollary 2 (b). We can notice that the histograms give very accurate approximations for the densities in Corollary 2. Also, in both Fig. 5 and Fig. 6, when the values of σ 1 or σ 2 increase, the histograms become right skewed. Moreover, for large values of σ 1 or σ 2 the population becomes extinct and µ 1 = δ 0 (see also Theorem 3.5).
6. Conclusions. We present a stochastic model to study the effect of toxicants on human cells. To account for parameter uncertainties, the model is expressed as a system of coupled ordinary stochastic differential equations. The variables are the cell index n(t), which closely reflects the cell population, the concentration C o (t) of internal toxicants per cell, and the concentration C e (t) of toxicants outside the cells at time t. There are a few papers that consider similar stochastic models for population dynamics, but they mainly study conditions for extinction and persistence.
Here we focus on the ergodic properties when the population is persistent. We first prove the positivity of the solutions. Then we investigate the influence of noise on the cell population survival. When the noise variances σ 2 1 or σ 2 2 are sufficiently large, the population goes to extinction. Numerical simulations show that, for the stochastic model, the population goes to extinction at threshold values C stoch e (0) below the deterministic threshold value C det e (0). Furthermore, increasing the noise variances σ 2 1 or σ 2 2 results in a lower value C stoch e (0) at which the population becomes extinct.
Moreover, we prove that when the noise variance σ 2 1 is sufficiently small and the population is strongly persistent, then the cell index converges weakly to the unique stationary probability distribution. Increasing the noise intensity causes a right skewness of the stationary distribution.
Here we illustrate our results for the toxicant monastrol. We have also considered other toxicants from the experiments described in [1] classified in various clusters [30]. We have noticed that the cluster type does not change the type of stationary distribution, nor has an effect on the behavior of the distributions in response to increased noise variances.
Appendix A. Proof of Lemma 2.2.