Dynamics of an Intra-host Diffusive Pathogen Infection Model

In this paper, we first propose a diffusive pathogen infection model with general incidence rate which incorporates cell-to-cell transmission. By applying the theory of monotone dynamical systems, we prove that the model admits the global threshold dynamics in terms of the basic reproduction number ($\mathcal{R}_{0}$), which is defined by the spectral radius of the next generation operator. Then, we derive a discrete counterpart of the continuous model by nonstandard finite difference scheme. The results show that the discrete model preserves the positivity and boundedness of solutions in order to ensure the well-posedness of the problem. Moreover, this method preserves all equilibria of the original continuous model. By constructing appropriate Lyapunov functionals for both models, we show that the global threshold dynamics is completely determined by the basic reproduction number. Further, with the help of sensitivity analysis we also have identified the most sensitive parameters which effectively contribute to change the disease dynamics. Finally, we conclude the paper by an example and numerical simulations to improve and generalize some known results.


Introduction
Over the past few decades, there has been a great effort in the mathematical modelling of within-host pathogen infection models. These models have been used to describe the dynamics inside the host of various infectious diseases such as HIV, HCV, HBV and HTLV, etc. The classical model for within-host virus dynamics is a system of three ordinary differential equations [1,2], where the key assumption is that cells and viruses are well mixed, and ignores the mobility of viruses. To study the influences of spatial structures of virus dynamics, Wang and Wang [3] proposed the following diffusive system by assuming the motion of virus follows the Fickian diffusion [4] (1) ∂I(x, t) ∂t =β 1 SV − d I I, where S, I and V are the concentrations of susceptible or uninfected cells, infected cells and free pathogens at the position x at time t, respectively. The susceptible cells are produced at a constant rate Λ and are infected by free virions at a rate β 1 SV . Parameters d S , d I and d V represent the death rates of uninfected cells, infected cells and free virus, respectively. The free virions are produced from the infected cells at a rate αI. D 3 is the diffusion coefficient and ∆ is the Laplacian operator. Notice that the above system (1) only focus on virus-to-cell spread in the bloodstream even though some literatures reveal that cell-to-cell (infected source cell and a susceptible target cell) transmission is vital to spread of virus in vivo [5][6][7][8]. An understanding of viral cell-to-cell spreading will enhance our ability to intervene in the efficient spreading of viral infections. For more information on dealing with target cell dynamics and cellto-cell transmission one can refer [9][10][11][12][13][14][15][16][17] and references therein. On the other hand, the bilinear incidence rate is a simple description of the infection in system (1). As mentioned in [14,18], a general incidence rate may help us to gain the unification theory by the omission of unessential details. Hence, inspired by the aforementioned work, in this article we propose the following pathogen infection model on domain Q = R + × Ω, ∂I(x, t) ∂t = D 2 ∆I + Sf (V ) + Sg(I) − (γ + d I )I, x ∈ Ω, t > 0, x ∈ Ω, t > 0.
Here D 1 and D 2 are the diffusion coefficients of susceptible and infected cells, respectively [19] with γ is the lysis rate of infected cells [20]. The incidences are assumed to be the nonlinear responses to the concentrations of virus particles and infected cells, taking the forms Sf (V ) and Sg(I), where f (V ) and g(I) denote the force of infection by virus particles and infected cells and satisfy the following properties [21]: (iii) the per capita infection rates by virus particles and infected cells will slow down due to certain inhibition effect since (3) implies that f (V ) V ≤ 0 and g(I) I ≤ 0. Obviously, the incidence rate with above conditions contains the bilinear and the saturation incidences such as f (V ) = β 1 V or β 1 V 1 + V and g(I) = β 2 I or β 1 I 1 + I , where incidence rates β 1 , β 2 > 0. In this paper, we consider the system (2) with initial conditions as follows where S 0 (x), I 0 (x), V 0 (x) ∈ C 2 (Ω) ∩ C 0 (Ω), and homogoneous Neumann boundary conditions where Ω is an open bounded subset of R n with piecewise smooth boundary ∂Ω and ν being the unit outer normal to ∂Ω. Generally, the exact solution for a system like (2) is difficult or even impossible to be determined. It is a natural requirement of an adequate numerical method that it possess the discrete equivalents of the qualitative properties the continuous system satisfies. However, how to select the proper discrete scheme so that the global dynamics of solutions of the corresponding continuous models can be efficiently preserved is still an open problem [22]. Actually, Mickens has made an attempt in this regard, by proposing a robust non-standard finite difference (NSFD) scheme [23], which has been widely employed in the study of different models [24][25][26][27][28][29]. Motivated by the work of [23], we apply the NSFD scheme to discretize system (2) and obtain (6)  Here, we set the spatial domain x ∈ Ω = [a, b] where a, b ∈ R and ∆x = (b − a)/M be the space step-size generating M equal sub-interval over the domain and the length of uniform time intervals be ∆t. We denote approximations of S(x n , t k ), I(x n , t k ) and V (x n , t k ) by S k n , I k n and V k n , respectively at each mesh point {(x n , t k ), n ∈ {0, 1, · · · , M }, k ∈ N} with x n = a + n∆x and t k = k∆t. The discrete initial and boundary conditions respectively. Our aim is to show that the discretized system (6) which derived by using NSFD scheme can efficiently preserves the global asymptotic stability of the equilibria to the original system (2). The rest of this paper is organized as follows. In Section 2, we study the dynamical behavior of the continuous system (2), such as the existence of positive solutions and its uniqueness, the existence of equilibria, basic reproduction number, local stability and global stability. In Section 3, we investigate the global dynamics of discrete system (6). Numerical simulations are carried out to validate the theoretical results in section 4 and a brief conclusion finishes the paper.
2. Dynamical behavior of the system (2) 2.1. Existence, Uniqueness and Positivity. To discuss the dynamical behavior of system (2), first we give the definition of upper and lower solution.
) are a pair of upper and lower solution to the problem (2), ifŠ ≤Ŝ,Ǐ ≤Î,V ≤V inΩ × [0, ∞) and the following differential inequalities hold: for (x, t) ∈ Ω × (0, ∞) and It is easy to see that 0 = (0, 0, 0) and K = (K 1 , K 2 , K 3 ) are a pair of coupled lower-upper solutions to problem (2), where and d = min{d S , d I }. Using the following lemma provided by Redinger [30], we get the existence and uniqueness of the solution.

Equilibria and Basic reproduction number.
It is easy to verify that system (2) always has a disease-free equilibrium E 0 (S 0 , 0, 0) with S 0 = Λ d S , and if exists the endemic equilibrium E * (S * , I * , V * ) satisfies In order to find the basic reproduction number (R 0 ) for the system (2), we obtain the following linear system at E 0 for the infected classes: Substituting I(x, t) = e λt ψ 2 (x) and V (x, t) = e λt ψ 3 (x) into (10), we obtain the following cooperative eigenvalue problem: By [31](Theorem 7.6.1), we conclude that (11) has a principal eigenvalue λ(S 0 , f (0), g (0)) with a positive eigenfunction. Now we are in a position to apply the ideas and the theory in [32] to define R 0 for the model (2). Let T : C(Ω, R 2 ) → C(Ω, R 2 ) be the solution semigroup of the following reaction-diffusion system: x ∈ Ω, t > 0, Thus, with initial infection Ψ(x) = (ψ 2 (x), ψ 3 (x)), the distribution of those infection members becomes T (t)Ψ(x) as time evolves. As in [32], the matrices F and V defined as Therefore, the distribution of total new infections is Then, we define It is clear that L is a positive and continuous operator which maps the initial infection distribution Ψ to the distribution of the total infective members produced during the infection period. Applying the idea of next generation operators [32], we define the spectral radius of L as the basic reproduction number By some calculations, we obtain that where R 01 and R 02 are partial basic reproduction numbers induced by virus-to-cell transmission and cell-to-cell transmission, respectively. The following theorem now prove regarding the meaningful steady states.
It is easy to proof for the case R 0 < 1. Consider R 0 > 1, it follows from (9) that Hence, equation G(I) = 0 has at least one positive root I * ∈ 0, Λ γ + d I . That implies the existence of positive equilibrium of the system (2). In order to show that the positive equilibrium is unique, we use Then According to equation (3), we have which implies that G (I * ) < 0. If there exists the second positive equilibrium E (S , I , V ), then one has G (I ) < 0. But which contradict the conditions (13). This completes the proof.
Proof. The linearization of system (2) at E 0 can be expressed by Therefore, the characteristic equation at E 0 is It is obvious that (14) has an eigenvalue λ 1 = −(d S + µ i D 1 ). The other two eigenvalues λ 2 and λ 2 are roots of It is easy to see that and Since R 02 < R 0 < 1, we have λ 2 + λ 3 < 1 and λ 2 λ 3 > 0. This gives that Re(λ 2 ) < 0 and Re(λ 3 ) < 0. Thus, all eigenvalues of (14) have a negative real parts when R 0 < 1. Hence, E 0 is locally asymptotically stable. This completes the proof.
Now we turn our attention to the endemic equilibrium E * .
Proof. Linearizing system (2) at E * gives Thus, the characteristic equation at E * is where From (9) and (13), we get , . Hence, , Hence, Q 2 Q 1 −Q 0 > 0. Then, by using Routh-Hurwitz criterion we claim that all eigenvalues of (15) have negative real parts. Thus, the endemic equilibrium E * of system (2) is locally asymptotically stable when R 0 > 1. This completes the proof.
2.4. Global Stability. Now, we discuss the global stability of the equilibria for the system (2) by considering Lyapunov functional based on the Volterra function Φ(x) = x − 1 − ln x. Clearly, Φ(x) ≥ 0 for all x > 0 and the equality holds if and only if x = 1. In presence of diffusion, the aim is to show that every solution of the system (2) with a positive initial value that is different from the equilibrium point will converge to the equilibrium.
Theorem 2.6. If R 0 ≤ 1, then the disease-free equilibrium E 0 of system (2) is globally asymptotically stable.
Proof. Define a Lyapunov function and B is a positive constant to be determined later. Then, along the solutions of the system (2), we have By equation (3) and choosing B = (γ + d I − S 0 g (0))/α, we obtain Using Green's formula and the Neumann boundary conditions in (5), we obtain and Ω ∆Idx = Ω ∆V dx = 0.
Using above conditions, we obtain Therefore, dL(t) dt ≤ 0 whenever R 0 ≤ 1. It follows that the largest invariant subset of dL(t) dt = 0 is the singleton E 0 . By LaSalle's Invariance Principle [35], the infection-free equilibrium of the system (2) is globally asymptotically stable when R 0 ≤ 1.
Next, we turn our attention to show the global stability of the endemic equilibrium E * .
Theorem 2.7. Consider a Lyapunov function Then, H(t) is non-negative and is strictly minimized at the unique equilibrium (S * , I * , V * ), i.e. it is a valid Lyapunov function. Hence, E * = (S * , I * , V * ) is globally asymptotically stable.
Proof. According to (2), we have Using Green's formula and the Neumann boundary conditions in (5), we obtain Using above conditions, we conclude that Furthermore we have dH(t) dt = 0 only at steady state E * = (S * , I * , V * ). Therefore, by Lyapunov's direct method, the steady state solution E * is globally asymptotically stable.

Dynamical behavior of the Discretized model (6)
In preceding section, we have shown that the global asymptotic stability of the equilibria for the continuous system (2) is completely determined by the basic reproduction number R 0 by constructing appropriate Lyapunov functionals. This arises a natural question that whether the global asymptotic stability of the equilibria of the discrete system (6) can be preserved. In this section, we will discuss this problem. Clearly, the discretized system (6) has the same two steady states as the continuous system (2). In the following theorem, we show the system (6) is non-negative and bounded. Proof. The positivity of the solutions of the discretized system (6) can be proved using the M-matrix theory [26]. From the first equation of system (6), we get It is clear that A k is a strictly diagonally dominant matrix. Thus, the first equation of system (6) is equivalent to From the second equation of the system (6), we have Since B is a M-matrix, we get Similarly, from the third equation of system (6), we have Since C is a M-matrix, we obtain Since all parameters in the system (6) are positive, it is easy to see that the solution remains non-negative for all k ∈ N. Next, we prove the boundedness of the solution. Define a sequence {G k } as follows: It follows from the first two equations of system (6) that where d = min{d S , d I }. Hence, we have By mathematical induction, we obtain lim sup k→+∞ Q k ≤ N (M + 1).
By the third equation of the system (6), we get Since {G k } is bounded, there exists a positive constant ξ such that M n=0 I k n ≤ ξ. Thus, we have This completes the proof.
3.1. Global Stability. In this section, we establish the global asymptotic stability of the steady states E 0 and E * of the discrete system (6), by constructing discrete Lyapunov functions.
which gives that {L k } k∈N is a monotone decreasing sequence. That means there exists a constantL ≥ 0 such that lim k→+∞ L k =L and we have lim k→+∞ (L k+1 − L k ) = 0. Arguing on the lines of [26] to system (6), we obtain lim k→+∞ S k n = S 0 , lim k→+∞ I k n = 0 and lim k→+∞ V k n = 0, for all n ∈ {0, 1, · · · , M } and when R 0 ≤ 1. This completes the proof.
Next, we discuss the global stability of the endemic equilibrium E * when R 0 > 1.
Proof. We define the following discretized Lyapunov function Obviously, H k ≥ 0 for all k ∈ N with the equality holds if and only if S k n = S * , I k n = I * and V k n = V * for all n ∈ {0, 1, · · · , M } and k ∈ N. Using ln x ≤ x − 1, the difference of H k Using the similar arithmetic-geometric inequality (16) for S, I and V , we get By assumption (A2), we have the following inequality Hence, using above conditions, we conclude that This implies that H k is a monotone decreasing sequence, then there exists a constantH such that lim k→∞ H k =H, and we have lim k→∞ (H k+1 − H k ) = 0. Which conclude that By the first equation of system (6), we obtain Taking k → +∞ in the above equality, we have Similarly, we get This completes the proof.
Thus, using the Theorems 3.2 and 3.3, we conclude that the discretized system (6) exhibits dynamic consistency with the continuous system (2) as well as the global asymptotic stability of both the uninfected and the infected steady states.

Numerical Results
In this section, we present numerical example that illustrate and confirm the findings of this study for the linear incidence function such as f (V ) = β 1 V and g(I) = β 2 I. Thus the system (2) becomes The basic reproduction number of the system (17) is .
The system (17) possesses an uninfected steady state E 0 Λ d S , 0, 0 and also has an infected At first sensitivity analysis is used to determine the response of the model to variations in its parameter values. In the present case, focus is given to determining how changes in the model parameters impact the basic reproduction number. This is done through the Latin hypercube sampling (LHS) and the partial correlation coefficients (PRCC) to determine the relative importance of the parameters in R 0 for the disease transmission [36]. In such a scenario, it is more appropriate to treat each parameter as a random variable, distributed according to an appropriate probability distribution. We assume that our model parameters are normally distributed although it is quite possible that some parameters are constant towards a particular value such as recruitment rate (Λ) and death rate (d T ) of susceptible cells. PRCC reduces the non-linearity effects by rearranging the data in ascending order, replacing the values with their ranks and then providing the measure of monotonicity after the removal of the linear effects of each model parameter keeping all other parameters constant [37]. The corresponding Tornado plots based on a random sample of 1000 points for the six parameters in R 0 are shown in Figure 1. The horizontal lines represent the significant range of correlation, i.e., |PRCC| > 0.5. The sensitivity analysis suggests that the most significant parameters are β 1 and β 2 , an increase in these values will have an increase in the spread of the disease. Hence, these parameters should be estimated with precision to accurately capture the dynamics of the infection. Now, we numerically illustrate the results for global stability for both the steady states. Accordingly, we use two sets of system parameters, one corresponding to R 0 < 1 (when E 0 is globally asymptotically stable for both the continuous and discretized models) and the other for R 0 > 1 (when E * is globally asymptotically stable for both the continuous and discretized models). The numerical simulation is carried out using the NSFD scheme described by the system (6) with initial condition taken as S(x, 0) = 10 7 , I(x, 0) = 100e x , V (x, 0) = 100e x .
For the purpose of illustration of both the scenarios, we choose the equal diffusion coefficients as D 1 = D 2 = 1 mm 2 d −1 and D 3 = 1 mm 2 d −1 [38]. The one-dimensional spatial domain is taken as Ω = [0, 50] and the simulation carried out for a time window of 100 days. The grid sizes used in the spatial and temporal directions are ∆x = 0.5 and ∆t = 1, respectively. The parameter set Λ = 10 7 cells d −1 , β 1 [26], results in R 0 = 0.21 < 1. Thus, Figure 1. Partial rank correlation coefficient (PRCC) results for significance of parameters involved in R 0 in this case, the uninfected steady state E 0 is globally asymptotically stable. It can be observed from Figure 2 that this is indeed the case and the system eventually approaches the uninfected steady state E 0 = (10 8 , 0, 0). For the other scenario, all the parameter values are identical with the exception of β 1 = β 2 = 3 × 10 −10 virion −1 d −1 [26] which renders R 0 = 12.59 > 1. In this case, the infected steady state is stable as can be observed numerically in Figure 3, where the state variables approach the infected steady state E * = (8 × 10 6 , 1.8 × 10 8 , 3.7 × 10 9 ). For both the sets of simulations it can be easily seen that the steady states do not depend on the initial spatial points.
We have also shown that the global asymptotically stable results are dependent only on the parameters of the non-diffusive system because of R 0 and independent of the choices of the diffusion coefficients. This is also illustrated by way of numerical simulations. For illustrative purpose we only show the case for R 0 > 1 using the corresponding parameter values used above. We extend our diffusion coefficient to D 1 = D 2 = D 3 = 100 mm 2 d −1 and see from Figure 4, that the steady states of the model dynamics are very similar to each other in the long run. Similar results can be observed for different combinations of . When R 0 = 12.59 > 1,the disease-free equilibrium E * of system (17) is globally asymptotically stable.
(D 1 , D 2 , D 3 ). Another crucial advantage of using the NSFD scheme over standard finite difference (SFD) scheme is that the positivity of solutions for long time simulation which already been demonstrated in many works [26,39].

Conclusion
In order to investigate the mechanism of virus infection and viral replication, we have carried out the mathematical analysis for a diffusive intra-host virus dynamics model with cell-to-cell transmission, while allowing for a general nonlinear incidence functions. The well posedness and linear stability of equilibria of this model is investigated. The basic reproduction number R 0 is a threshold index which predicts the extinction and persistence of the disease. It is shown that the global stability of the equilibria is completely determined by R 0 : if R 0 ≤ 1, then the disease-free equilibrium E 0 is globally asymptotically stable stable, which means that the virus is eventually cleared and the infection dies out; if R 0 > 1, then the endemic equilibrium E * is globally asymptotically stable. From the expression of the basic reproduction number of virus, we can know that the basic reproduction number could be under-evaluated without considering either the virus-to-cell transmission or cell-to-cell transmission; it is not enough to eliminate the disease by decreasing the basic reproduction number of virus-to-cell transmission due to the existence of cell-to-cell transmission. Our results also imply that diffusion coefficients have no effect on the global behaviors of such virus dynamics model with homogeneous Neumann boundary conditions. Then, by using the NSFD scheme, we derived the discrete counterpart of the continuous model. Our results show that the discretization scheme can preserve the global properties of solutions for original continuous model, including the positivity, ultimate boundedness and global stability of the equilibria with no restriction on the space and time step sizes.
In addition, the model proposed in this paper is an extension of some previous works [40][41][42][43][44][45] and the obtained results improve some known results. It is interesting to improve the model (2) by incorporating logistic growth term for uninfected target cells and a more general infection function. We leave these for future consideration.