STOCHASTIC VIRUS DYNAMICS WITH BEDDINGTON-DEANGELIS FUNCTIONAL RESPONSE

Stochastic virus dynamics modeled by a system of stochastic differential equations with Beddington-DeAngelis functional response and driven by white noise is investigated. The global existence of positive solutions and the existence of stationary distribution are proved. Upper and lower bounds of the pathwise and asymptotic moments for the positive solutions are sharply estimated. The absorbing property in time average is shown and the moment Lyapunov exponents are proved to be nonpositive.


Introduction
In this paper we shall consider the stochastic viral dynamics modeled by the following stochastic differential equations dx = λ − δx − βxz 1 + ax + bz dt + σ 1 x dB 1 (t), dy = βxz 1 + ax + bz − qy dt + σ 2 y dB 2 (t), dz = (ky − γz)dt + σ 3 z dB 3 (t), (1.1) where the susceptible cells whose density is denoted by x(t) are generated at a constant rate λ, die at a density-proportional rate δx, and become infected with a rate βxz/(1+ax+by) in terms of the Beddington-DeAngelis functional response, the infected cells of the density y(t) are produced at the same rate from the susceptible cells and die at a density-dependent rate qy, while the virus-free cells z(t) are released from the infected cells at a rate ky and die at a rate γz, cf. [4,8,16]. All these involved parameters and the intensity coefficients σ i 's of the stochastic driving terms are positive constants.
In the classical infectious disease model, a bilinear incidence as βxz was often used. However, the actual incidence rate is most likely not linear over the entire range of the densities x and z. The nonlinear incidence rate with this Holling type II response βxz/(1 + ax + bz) was introduced by Beddington [1] and DeAngelis et al. [3] in the study of HIV-1 virus and uninfected CD4 + T cells. This incidence rate is more reasonable than a bilinear rate because it incorporates the behavioral change and crowding effect of the infective cells or particles.
As indicated in recent researches on epidemic dynamics, infection dynamics, as well as population dynamics, cf. [2,4,5,14,16,17] and references therein, it is found that the deterministic models are subject to some limitations in modeling the virus transmission of infectious or epidemic diseases. Virus or disease transmission dynamics are influenced by the random effect of environmental noise, immunodeficiency, antibiotic resistance, even weather conditions, or uncertain fluctuations. These researches also showed that stochastic differential equation models with the white noise generated by Brownian motion can provide some additional degree of realism in dealing with continuous fluctuations of randomness affecting the birth and death rates, transmission coefficients, and other parameters in the system . The stochastic models with this kind of Wiener noise likely serve better in predicting the future dynamics than the deterministic models. It is commented in [5, page 879] that adding the stochastic Wiener noise to the parameters of relevant rates is a well-established way of introducing stochastic environmental noise into biologically realistic population dynamic models. Thus it is reasonable to include the additive or multiplicative white noises to the death rates as −δx + σ 1 x dB 1 (t), −qy + σ 2 y dB 2 (t) and −γz + σ 3 z dB 3 (t) in the model equations (1.1) for tackling the effect of environmental fluctuations.
The global stability of virus dynamics for the corresponding deterministic system was analyzed in [8] and the following results are proved: When the comprehensive reproductive ratio of the virus the disease-free equilibrium E 0 = (λ/δ, 0, 0) is globally asymptotically stable; When R 0 > 1, then the endemic equilibrium E 1 = (x 0 , y 0 , z 0 ) is globally asymptotically stable, where Dynamics of a stochastic predator-prey system with Beddiongton-DeAngelis functional response and multiplicative or additive noise was studied in [9,18]. It was shown that if the noise is large, both species in the system may go to extinction, which does not occur for the deterministic system. Besides, in [21] and [27], stationary distribution of stochastic population systems and stochastic Hopf bifurcation of the predator-prey system with Beddington-DeAngelis response have been investigated. Stochastic dynamics of predator-prey system with Leslie-Gower response and lévy jumps was studied in [19].
The main objective of this article is to study the pathwise, time-averaging, and asymptotic dynamics for the almost surely positive solutions in terms of construction of Lyapunov functions, proving the existence of a stationary distribution, conducting sharp estimates of stochastic upper and lower bounds of trajectories, application of the ergodic property, and evaluation of the moment Lyapunov exponents in regard to stochastic stability.
The basic concepts and some underlying results for stochastic ordinary differential equations can be found in [10] and [22]. In this paper, let (Ω, F, {F t } t≥0 , P ) be a complete probability space with a filtration {F t } t≥0 , which is right continuous and F 0 contains all P -null sets. Assume that B i (t), i = 1, 2, 3, are independent standard Brownian motion defined on this probability space. We use the notation For an n-dimensional stochastic differential equation driven by the white noise, which will be called nonautonomous if functions f and g are explicitly dependent of time t, and continuously differentiable in x to the second order and in t to the first order. Define the differential operator L associated with (1.3) by Then the action of the operator L on where V t = ∂V ∂t , V x = ( ∂V ∂x1 , · · ·, ∂V ∂xn ) and V xx = ( ∂ 2 V ∂xixj ) n×n . By the Itô formula, we have the following differential formula, which corresponds to the the total derivative in deterministic case, Definition 1.1. [10,21] Let P X0,t (·) be the probability measure induced by a stochastic process {X(t)} t≥0 in R n + (or R n ) over a probability space (Ω, F, P) with initial state X(0) = X 0 , namely, where B(·) stands for the σ-algebra of all Borel sets. If there is a probability measure µ(·) on the measurable space (R n + , B(R n + )) such that in distribution for any X 0 ∈ R n + , then X(t) has a stationary distribution µ(·).
In Section 2, we shall prove the global existence of positive solutions to the system of the model equations (1.1) for the stochastic virus dynamics. In Section 3, it will be shown that a stationary distribution exists. In Section 4, we present the pathwise and asymptotic moment estimation of the upper and lower bounds of the positive solutions. In Section 5, we show the absorbing property of in time average. In Section 6, we show that the moment Lyapunov exponents are nonpositive for this dynamical system.

Global Existence of Positive Solutions
First we prove the existence and uniqueness of a global positive solution for the initial value problem of the SDE (1.1).
Theorem 2.1. Under the condition q ≥ k, for any initial data (x 0 , y 0 , z 0 ) ∈ R 3 + , there exists a unique positive solution (x(t), y(t), z(t)) to the system (1.1) such that the solution will remain in R 3 + for all t ≥ 0 with probability one. Proof. Since the coefficients of the equations in (1.1) are locally Lipschitz continuous, for any ω ∈ Ω and any given initial data (x 0 , y 0 , z 0 ) ∈ R 3 + , there is a unique local solution (x(t), y(t), z(t)) on t ∈ [0, τ e ), where [0, τ e ) is the maximal existence interval pathwise depending on the initial data.
We now show that the solution exists globally, namely, τ e = ∞ a.s. Let an integer m 0 > 0 be sufficiently large such that (x 0 , y 0 , z 0 ) ∈ [1/m 0 , m 0 ] 3 . For each integer m > m 0 we define the stopping time It suffices to prove that τ ∞ = ∞ for any initial data a.s. If this statement is false, then there would be constants ε ∈ (0, 1) and T > 0 such that Consider the function V (x, y, z) defined by Using the Itô formula, we have where the condition q ≥ k implies that (2.5) Therefore, we obtain By integrating both sides of (2.6) from 0 to τ m ∧ T and then taking the expectation on both sides, it yields Let Ω m = {τ m ≤ T }, then P (Ω m ) ≥ ε. Note that for every ω ∈ Ω m at least one of x(τ m , ω), y(τ m , ω), z(τ m , ω) equals to m or 1/m, then V (x(τ m , ω), y(τ m , ω), z(τ m , ω)) is no less than Thus we get It follows that where I Ωm (ω) is the indicator function of Ω m (ω). Let m → ∞. Then it leads to a contradiction ∞ > V (x(0), y(0), z(0)) + KT = ∞.
Therefore we have τ e = ∞ for any given initial data almost surely, which completes the proof. According to [22,Theorem 2.9.1], this positive solution X(t) = (x(t), y(t), z(t)) T is a Markov process with time-homogeneous transition function.

Stationary Distribution
The system of stochastic autonomous differential equations (1.1) can be written in the vector-matrix form Its diffusion matrix is given by Here we shall prove the existence of a stationary distribution for a positive solution of the system (1.1). We mention the following assumption in cf. [10, p.107].
Assumption U. There is a bounded open domain U ⊂ R 3 + with regular boundary and the following properties: I. In a neighborhood of U , the smallest eigenvalue of the diffusion matrix A(x, y, z) is uniformly bounded away from zero.
II. For any X = (x, y, z) ∈ R 3 + \U , the mean time τ at which a path from X reaches the set U is finite and sup K E X [τ ] < ∞ for every compact set K ⊂ R 3 + . Lemma 3.1. Suppose the Assumption U is satisfied. Then the Markov process X(t, X 0 ) given by the pathwise solution of (3.1) has a stationary distribution µ(·) such that µ(R 3 \ R 3 + ) = 0 and for any nonnegative continuous function Q(X) it holds that Moreover, for any integrable function F (X) with respect to µ, it holds that The proof of this lemma is referred to Theorems 4.4.1 and 4.4.2 in [10].
Theorem 3.1. If there exists a positive equilibrium point (x 0 , y 0 , z 0 ) ∈ R 3 + for the deterministic system (1.2), then there exists a stationary distribution with respect to the positive solutions of the system (1.1).
Thus the first condition in the Assumption U is also satisfied. Therefore, by Lemma 3.1, the system of SDE (1.1) has a unique stationary distribution µ(·) on R 3 + . The ergodic property in terms of this stationary distribution will be shown in Section 5.

Asymptotic Moment Estimation
First we prove an auxiliary lemma useful for dynamic estimation of the moments of the positive solutions to the stochastic viral system (1.1).
Proof. Divided the equation (4.1) by v 1− 1 p (t), we get , then we come up with the linear ordinary differential equation and the solution (4.2) is obtained. We now derive an asymptotic upper bound for the p-th moment of each component solution of the system of SDE (1.1).
Theorem 4.1. Suppose the following conditions are satisfied, Let (x(t), y(t), z(t)) be a solution of the system of stochastic viral equations (1.1) with any initial data (x 0 , y 0 , z 0 ) ∈ R 3 + . Then for any p > 1 it holds that .
Proof. By the Itô formula, for the x-component solution of (1.1) we have Consequently, (4.5) Taking the expectation on both sides of (4.5), we have Therefore, the p-moment of the x-component solution satisfies the differential equation Define u(t) = E[x p ] and we have According to the condition (4.3), we have δ + Hence, (4.10) By Lemma 4.1 and its proof, the solution of the following Bernoulli equation .
Thus the solution of (4.11) is given by According to the condition (4.3), we have q − σ 2 2 2 (p − 1) > 0. Due to the exponential decay and by the comparison theorem, we get (4.12) Finally, for the z-component solution we have and (4.14) Similarly, since γ + (4.15) The proof is completed.
Next we derive specific pathwise upper bound and lower bound for each of the three component solutions of the system (1.1) for the study of stochastic virus dynamics and stability. Theorem 4.2. Every positive solution (x(t), y(t), z(t)) of the system (1.1) with any initial data (x 0 , y 0 , z 0 ) ∈ R 3 + satisfies the estimates (4.16) where the upper and lower bounds Φ u (t), Φ l (t), Ψ u (t), Ψ l (t), Γ u (t) and Γ l (t) are the stochastic processes shown by (4.17) through (4.22) below.
Proof. Since the solutions are positive in probability one, from (1.1) we have We can assert and verify that the stochastic process given by is the unique solution of stochastic differential equation with the initial condition Φ(0) = x 0 . By the comparison theorem for the pathwise solutions of stochastic differential equations, it holds that For the y-component solution, we have dy = βxz 1 + ax + bz − qy dt + σ 2 y dB 2 (t) ≥ (−β − qy) dt + σ 2 y dB 2 (t).
Here we can claim and check that the stochastic process given by is the unique solution of the initial value problem with Ψ(0) = y 0 . By the comparison argument, we obtain On the other hand, Similarly, by the comparison argument, we can get the stochastic process given by as the upper bound for the y-component solution, For the component solution z(t), we have Accordingly, the stochastic process given by is the unique solution of the initial value problem with Γ(0) = z 0 and we have On the other hand, Thus we have (4.21) Finally, by the first stochastic differential equation in (1.1) and z(t) ≤ Γ u (t), we notice that Hence we can deduce by the similar argument that (4.22) The proof is completed.

Absorbing Property in Time Average
Here we prove the absorbing property of the positive solutions in time average for the stochastic viral system (1.1) based on the asymptotic moment estimation and the stationary distribution shown in previous sections.
where (x(t), y(t), z(t)) is any solution of the system (1.1) with initial data in R 3 + and L i (p), 1 ≤ i ≤ 3, are given in (4.4). Consequently, the p-th moments of the pathwise positive solutions exist for p > 1 with respect to the stationary distribution µ.
Proof. By the Birkhoff ergodic property, for any N > 0, it holds that Using the dominated convergence theorem and by Theorem 4.1, we have E lim where L i (p), 1 ≤ i ≤ 3, are given in (4.4). Therefore, Let N → ∞, then we obtain where we used (3.2) and (3.3) in Lemma 3.1. Thus (5.1) is proved.

Moment Lyapunov Exponent
Lyapunov exponent and moment Lyapunov exponent are important characteristics for study of the almost sure stability and the moment stability of a stochastic dynamical system generated by stochastic differential equations. However, the actual evaluation of the moment Lyapunov exponents will be very difficult except by few approximation approaches of the asymptotic expansions, cf. [10,Appendix B].
Here for this stochastic virus model, we can make a calculation in the following result based on the upper and lower bound estimates. Theorem 6.1. Under the condition δ > σ 2 1 /2, q > σ 2 2 /2 and γ > σ 2 3 /2, the meansquare moment Lyapunov exponent of the positive solution trajectories of the system (1.1) is nonpositive, where X(t, X 0 ) = (x(t, x 0 ), y(t, y 0 ), z(t, z 0 )) is the solution with any initial data X 0 ∈ R 3 + . Proof. Step 1. (4.17) shows that Note that is a geometric Brownian motion, which is a solution of the linear SDE Accordingly we have Then by the Fubini theorem and the Jensen inequality of integral as well as the time-shifting property of Brownian motion, we can derive In view of the simple property log(1 + x) ≤ x for all x > −1 and δ > σ 2 1 /2, it follows that lim sup Step 2. For the y-component of the solution, from (4.19) we get Since B 1 (t), B 2 (t), B 3 (t) are independent Brownian motions, from (6.3) and (6.4) we can deduce Substitute (6.7) into (6.6). In view of the condition δ > σ 2 1 /2 and q > σ 2 2 /2, we obtain Therefore, similar to (6.5), we have lim sup where .

Conclusion
The aim of this paper is to investigate the virus dynamics modeled by the susceptibleinfected-removed (SIR) equations with the Beddington-DeAngelis functional response and the stochastic multiplicative and independent white noises. Here we studied the pathwise, time-averaging, and asymptotic dynamics for the almost surely positive solutions of this system. Starting with the proof of global existence of the positive solutions by construction of a linear-log Lyapunov function and using the stochastic Itô formula, the existence of a stationary distribution with respect to the positive solutions is established by construction of another sophisticated Lyapunov function. Through conducting sharp estimates, we obtained the asymptotic upper and lower bounds of the moments of stochastic trajectories.
Based on the asymptotic moment estimation and the stationary distribution as well as the ergodicity, we further proved the absorbing property in time average for the solution trajectories of this stochastic virus model. Finally we made the nonpositive evaluation of the mean-square and any order moment Lyapunov exponents. These results will be useful for a potential study in regard to prediction of virus persistence and extinction asymptotically in a long run.
Recently there have been stochastic epidemic models driven with noise of Lévy processes (sometimes called telephone noise or telegraph noise [2,11,23,25]) proposed and studied. Many good results on epidemic dynamics and prediction of persistence and extinction of diseases or virus for SIR and SIRS type equations with Markov switching have been reported, cf. [6,7,12,13,25,26] and references therein. Besides there are researches of longtime dynamics on the budworm growth and predator-prey models with harvesting and distributed delays also involving stochastic noises of regime-switching, cf. [2,15,20,28]. It is expected that the Beddington-DeAngelis dynamics with Lévy type noise will also be investigated.