Dynamics of a diffusive vaccination model with therapeutic impact and non-linear incidence in epidemiology

In this paper, we study a more general diffusive spatially dependent vaccination model for infectious disease. In our diffusive vaccination model, we consider both therapeutic impact and nonlinear incidence rate. Also, in this model, the number of compartments of susceptible, vaccinated and infectious individuals are considered to be functions of both time and location, where the set of locations (equivalently, spatial habitats) is a subset of with a smooth boundary. Both local and global stability of the model are studied. Our study shows that if the threshold level the disease-free equilibrium is globally asymptotically stable. On the other hand, if then there exists a unique stable disease equilibrium . The existence of solutions of the model and uniform persistence results are studied. Finally, using finite difference scheme, we present a number of numerical examples to verify our analytical results. Our results indicate that the global dynamics of the model are completely determined by the threshold value .


Introduction
Mathematical and epidemiological models are important tools for analysing various real world phenomena in health science and epidemiology. For infectious diseases, many mathematical and epidemiological models have been studied by researchers to understand the effect of vaccination for controlling the spread of infectious diseases. Diffusive vaccination models are useful models for analysing the impact of vaccination for infectious diseases. Moreover, diffusive vaccination models are useful for getting information about how to control the reasoning individuals.
It is known that vaccine works with the immune system. Evidently as the disease can not provide immunity, so not the vaccination. As a result, most of the diseases have a recovered/immune stage for which vaccination is successful. Some other bacteria can remain in [10] studied the following deterministic epidemic model with non-linear incidence H(I) = ( In the above model, the authors introduced the parameter c, the therapeutic treatment coverage of infectious individuals I(t) removed to S(t) compartment. Note that the above model is an SIS model and it was shown that the effectively treated infectious individuals return to the susceptible compartments and behaves similarly. The authors also observed realistically that q 2 ≤ q 1 from the fact that vaccination can reduce or eliminate the incidence of infection. Also, Gumel and Moghadas [10] analysed the corresponding characteristic equation and studied the local stability of its disease-free and disease equilibria and the optimal vaccine coverage threshold needed for disease control and eradication analytically. In 2014, Buonomo et al. [5] constructed suitable Lyapunov functions and established global stability of disease-free and disease equilibrium of the above system (1) by using LaSalle's invariance principle [16]. The authors also presented optimal vaccination and treatment strategies to minimize both the disease burden and intervention.
Recently, many researchers have considered spatial structure as a central factor because it affects the spatial spreading of disease [1,2,6,14,15,24,38,39]. In this paper, we propose a spatially dependent vaccination model which is a diffusive version of the above model (1), where we consider the individual movements of all three compartment cells. We strongly believe that our proposed model is a more general and realistic biological and epidemiological model. Throughout the paper, we use the following notation for simplicity: A = × (0, ∞) and ∂A = ∂ × (0, ∞). In the following, we present our proposed spatially dependent vaccination model with nonlinear incidence S(x, t) − (m + n) S(x, t) + cI (x, t) in A, with the following initial values and the zero-flux Neumann boundary conditions where ∂ ∂ω denotes the outward normal on ∂ . The Neumann boundary conditions imply that the populations do not move across the boundary ∂ or the population going out and coming in are equal on the boundary. It is also noted that S(x, t), V(x, t), I(x, t) are the number of compartments of susceptible individuals, vaccinated individuals and infectious individuals at time t > 0 and in location x ∈ , respectively. The notion is a spatial habitat in R n with a smooth boundary ∂ , is the usual Laplacian Operator, and δ 1 , δ 2 and δ 3 are the diffusion rates of susceptible, vaccinated and infectious compartments respectively. Since the model monitors dynamics of population, it follows that all its dependent variables and parameters, for examples, a, b, c, m, n, q 1 and q 2 must be non-negative as in the non-spatial model (1). We also set the upper bound of c as c < c * , which can be found in the proof of Lemma A.1.
A schematic representation of the model (2) is shown in the following Figure 1.
One of the fundamental issues in the study of infectious diseases via mathematical and epidemiological models is to find the stability of the two constant equilibria, that is, diseasefree equilibrium and disease equilibrium. In this paper, we study both local and global stability of our model. Our study shows that if the threshold level R 0 ≤ 1, the disease-free equilibrium E 0 is globally asymptotically stable. On the other hand, if R 0 > 1 then there exists a unique stable disease equilibrium E * . The existence of solutions of the model and the uniform persistence results for the model are studied. Finally, using finite difference scheme, we present a number of numerical examples to verify our analytical results. Our results indicate that the global dynamics of the model are completely determined by the threshold value R 0 .
The paper is organized in the following manner. In Section 2, we present disease-free and disease equilibrium respectively. Moreover, we present basic reproduction number in Section 2. We present our main results in Section 3. In Section 4, we present a number of numerical examples to verify our analytical results using finite difference scheme. Bifurcation results are also supported with parameter varying graphs. In section 5, we present existence and uniqueness of the solution of the system (2), local and global steady states along with responsible constraints are presented. Uniform persistence theorems for the model (2) are also highlighted as an interplay of our study. Finally, Section 6 discloses the summary of the results.

Preliminaries
For a deep look in the dynamics of the system (2), in this section, we will keep an eye on the basic reproduction number, the expected number of secondary cases reproduced by one infected individual in its entire infectious period.

Disease-free equilibrium
To define the disease-free equilibrium (S 0 , V 0 , I 0 ) of the system (2), we write the diffusion rates δ i = 0, since disease-free equilibrium state is not spatially dependent; then It is noted that for the disease-free equilibrium, we consider the count of compartments of infectious individuals I 0 = 0. Then we find, This gives the disease-free equilibrium: Let us now find the disease equilibrium of the governing system (2).

Disease equilibrium
In the case of equilibrium state, we have the disease equilibrium (S * , V * , I * ), where the diffusion rates δ i = 0. Then we write (2) as Here, the number of compartments of infectious individuals I * = 0. Then, we find the count of susceptible individuals in the form and the vaccinated individuals V * = n(1 + I * ) 2 (a + cI * ) (bq 1 I * + (m + n)(1 + I * ))(bq 2 I * + m(1 + I * )) . ( 8 ) Then, for the count of infectious individuals, we get the following polynomial of degree two where The real positive roots of (9) define the count of infectious individuals I * ; where the constant term of the quadratic Equation (9) Thereby, when R 0 > 1, we get the unique disease equilibrium E * (S * , V * , I * ) of the model (2). Now, from (7) and (8) we claim that and similarly for I * The proof of these claims are given in Lemma A.1 in Appendix.

Basic reproduction number
The Jacobian matrix of the linearized model (2) at E 0 is: with eigenvalues λ 1 = −(m + n), λ 2 = −m and λ 3 = Since all the model parameters are positive, it can be easily observed that λ 1 , λ 2 < 0. Thus, the equilibrium E 0 is locally asymptotically stable provides λ 3 < 0. Hence, by the definition of basic reproduction number [3], R 0 of (2) is For the sake of comprehension and clarity, we state our key results in the following section.
is non-positive or is dominated by negative values in the responsible Lyapunov function.

Remark 3.1:
See the last part of the proof of Theorem 3.4 in Section 5. The result in [5] that corresponds to Theorem 3.4, and on whose proof the proof of Theorem 3.4 is based, simply requires c = 0.

Theorem 3.5:
The proofs of the Theorems 3.1-3.5 are formulated through a series of steps in the Section 5. At this stage, first we want to justify all the key results by considering several numerical examples.

Examples and applications
For numerical verification for our analytic work, we choose finite-difference method based on Crank-Nicolson implicit time difference [8,17]. We can nicely observe the simulation part of the model (2) by using some graphical presentations. We take the initial conditions as: and the boundary condition is: Let us assume the diffusion rates δ 1 = δ 2 = δ 3 = 1.
Then the formula (10) gives us the basic reproduction number as R 0 = 22.9545 > 1 which ensures by Theorem 3.4 that, these values of parameters leads us to the disease equilibrium results as shown in Figure 3.

Parameter bifurcation observations
Now we are interested to know how the system (2) responses for different values of the system parameters.   (Figure 4), we clearly see that the disease is being extincted faster as c is increasing. But when c is more than 24.0 then c has no valuable effect for the disease for this parametric setup and more interestingly we get a cusp at t = 0.

From these Figures
Though, in our system (2) we assumed c to be non-negative anyhow; but if the disease causing environment still predominates, then we may consider c to be negative, for example, c ∈ [−m, 0). And in that scenario, we get the following results, Figure 5 shows that, if c is negative i.e. in disease causing environment when basic reproduction number R 0 < 0.001, then it is a disease-free equilibrium while R 0 > 0.001 reveals disease equilibrium. We also see the infection is increasing in a constant rate very roughly when R 0 is undefined in the case of c = −m.
Here, in Figure 6, we clearly observe the impacts of vaccination coverage parameter n over susceptible (S) and vaccinated (V) individuals. Susceptible (S) count converges to a minimum level and vaccinated (V) count increases to a maximum level as n growing large. But infectious (I) count remains approximately same for each cases.

Existence and uniqueness of solution
In this portion, we prove the existence and uniqueness of the solution of the system (2) by learning the algorithm partially from a similar study of Xu et al. [36].
Let us denote the subset of R 3 with vectors x ≥ 0 as R 3 + and X := C( , R) be a Banach space with the supremum norm · X . Also we define X + := C( , R 3 + ) then (X, X + ) is a strongly ordered space. Suppose that is the C 0 semigroups associated with δ 1 − (m + n), δ 2 − m and δ 3 − (m + c) subject to the Neumann boundary conditions, respectively. Then it follows that for any ρ ∈ C( , R) and t ≥ 0 where, i , i = 1, 2, 3 are the Green functions associated with δ i , i = 1, 2, 3, subject to the Neumann boundary conditions, respectively. It then follows from [29] that the function is compact and strongly positive. Particularly, is a strongly continuous semigroup.
Using these operators, we can write (2)-(4) as the following integral equation where, It can also be rewritten as the following abstract differential equation where, u = (S, V, I) and ρ = (S 0 , V 0 , I 0 ). Since F(ρ) is local Lipschitz continuous on X + , it then follows that for any ρ ∈ X + , (11) admits a unique noncontinuous mild solution u(·, t, ρ) such that u(·, t, ρ) ∈ X for all t in its maximum interval of existence. Moreover, it follows from ( [35], Corollary 2.2.5) that u(·, t, ρ) is a class solution of (2) with Neumann boundary conditions (4) for all t > 0. Further, by the scalar parabolic maximum principle, we see from the equation in (2) that S(x, t), V(x, t) and I(x, t) are all non-negative. Therefore, we obtain the following basic result on solution of the governing system (2)-(4).

Lemma 5.1: For any initial value function
Next, we show that the solution of the system (2)-(4) with the initial value function ρ ∈ X + actually exists globally, that is, σ = ∞. To this end, we need the following result ( [21], Lemma 5.1).
Consider the following reaction-diffusion equation where D > 0, A > 0 and d > 0 are positive constants.

Lemma 5.2:
The system (12) admits a unique positive steady state v * = A d which is globally attractive in C( , R). Now we are ready to produce the proof of the Theorem 3.1.
By (14), for any ρ ∈ X + , we see that there exist some t 1 = t 1 (ρ) > 0 such that Now, according to (13), as the first equation of (2) is local Lipschitz continuous on X + , it can easily be said that, for any ρ ∈ X + , there exist some t 1 = t 1 (ρ) > 0 such that Then by the similar argument as above, we also show that there are M i > 0, independent of the choice of ρ ∈ X + , and Therefore, the non-negative solution of (2)-(4) is ultimately bounded with respect to the maximum norm. This means that the solution semiflow (t) : X + → X + defined by ( (t)ρ)(x) = u(x, t, ρ), x ∈ , is point dissipative. In view of [35], (t) is compact for any t > 0. Thus, [11] implies that (t) : X + → X + , t ≥ 0, has a global compact attractor in X + . This completes the proof.

Analysis of local steady states
In this section, we want to explain the local stability of the equilibria for the system (2). Thus we consider the proof of our second result, Theorem 3.2. (2) at E 0 , we get

Proof of Theorem 3.2.: By linearizing the system
Then, we can obtain the following characteristic polynomial where, λ is the eigenvalue which determines temporal growth, I is the 3 × 3 identity matrix and L is the wave-number [24]. Then, we have Now, it is clear that It follows from R 0 < 1 that E 0 is locally asymptotically stable.
In the following, we prove the second part of the theorem. Linearizing the system (2) at E * , we obtain where, Then we obtain the following characteristic equation where, Now, let us take then we can get These lead us to the following conclusion By the Routh-Hurwitz criterion, we know that all eigenvalues of (16) have negative real parts. It means that the disease equilibrium E * of system (16) is locally asymptotically stable when R 0 > 1.

Global stability analysis
In this section, we investigate the global stability of the two constant equilibria E 0 and E * in the case of a bounded domain in which (S(x, t), V(x, t), I(x, t)) is an arbitrary positive solution of the system (2). First, let us consider the following shortcuts for convenience In case of global analysis, we consider the Lyapunov functional and the results varies with basic reproduction number. We stated two important results in the earlier Section 2.
At this phase, we are in stable setting to establish the Theorem 3.3 as long as the basic reporduction number R 0 ≤ 1.

Proof of Theorem 3.3.: Let define a Lyapunov function as
where, Calculating the time derivative of W 1 (x, t) along the solution of (2) gives Then from (2), we can write But, as a = (m + n)S 0 and mV 0 = nS 0 , we can write By Green's formula and Neumann boundary conditions (4), we get Similarly, Again, by Green's formula and the Neumann boundary conditions (4), we have the Green's first identity as By the same arguments, we also can write and Then using the above arguments, we have Recall the Equation (A6) which is described in the proof of Theorem A.1 (Appendix) Since c > 0, then using (23) the last integral of (22) satisfies Hence, dV 1 dt < 0 whenever R 0 ≤ 1. And, when S = S 0 , V = V 0 , I = 0; we calculate, dV 1 dt = 0 and vice-versa. Consequently, the singleton E 0 is the greatest compact invariant set in {(S, V, I) ∈ C( , R 3 + ) : dV 1 dt = 0}. Then, LaSalle's invariance principle [12] refers to lim t→∞ (S(x, t), V(x, t), I(x, t)) → E 0 ; which means, whenever R 0 ≤ 1, the disease-free equilibrium E 0 = (S 0 , V 0 , 0) is globally asymptotically stable. This establishes Theorem 3.3.
In a similar manner, it is stated that the disease equilibrium of (2) is globally asymptotically stable and the proof is prescribed as follows: Proof of Theorem 3.4.: Let us define a Lyapunov function as where, Calculating the time derivative of W 2 (x, t) along the solution of (2) gives Then from (2), it can written as Note that from (6), we have a = bq 1 I * 1 + I * S * + (m + n)S * − cI * , and by substituting these in (24) yields For writing convenience, let assume, f (I) = I 1+I such that Applying the Green's formula and zero Neumann boundary conditions, we obtain We know the arithmetic mean is greater than or equal to the geometric mean. Consequently, for all S > 0, V > 0 and I > 0, we find When Case (a) is true for all t ∈ (0, ∞), or for at-least large t > t 1 or t → ∞, the situation is clearly in favour and the result is well established.
But for Case (b) to be true, the integral function Z coincides with our expected result if the rest part of (25) equates or dominates on Z for all t ∈ (0, ∞), or for at-least large t or t → ∞.

Uniform persistence
By linearizing the third equation of system (2) at E 0 , the disease-free equilibrium, we get the followings: Then referring the arguments as in the proof of ( Now substitutingρ(x) ≡ 1 and the values of S 0 , V 0 into (5.4) we obtain the principal eigenvalue of (26) corresponding to which there is the unique positive eigen-functionρ(x) ≡ 1.
Thus, observing this equation we can claim the following lemma:
Proof: From the system (2), it is clear that S(x, t, ρ) > 0 and V(x, t, ρ) > 0 in A for any ρ ∈ X + . Then, Now applying ( [21], Lemma 1) and the comparison principle, we get Then there exists a t 1 > 0 such that Consequently, the second equation of system (2) follows that lim t→∞ inf V(x, t) ≥ an 2(bq 1 + m + n)(bq 2 + m) , uniformly for x ∈ , Finally, from the third equation of the system (2), we can write By the strong maximum principle and the Hopf boundary Lemma [25], this validates the second part.
After the completion of the above arguments, we obtain the results for disease persistence as described in Theorem 3.5 in Section 2. Now, it is time to produce the last result, Theorem 3.5 when the disease are persisting.
Since Theorem A.1 from appendix proves existence of global solution for the system (2) with distinct diffusion rates, the persistence theorem is also true for the system (2) where the diffusion rates (δ 1 , δ 2 , δ 3 ) are not identical and we describe the following statement as a remark.

Conclusion
In this manuscript, a spatially dependent vaccination model is proposed for infectious diseases. We have studied analytic inter-locution of disease-free equilibrium, disease equilibrium, basic reproduction number, existence and uniqueness of the solution of the corresponding system, local stability, global stability and uniform persistence theorem for the system. We present a number of numerical examples to verify our analytical results. It is shown that the numerical solution of the system corresponds to the analytical results. Our study may help to predict the upcoming probable results of treatments via vaccination and therapy against malignant diseases. Now, we take G = I − R 0 such that Therefore, I ≤ R 0 . ( A 8 ) Which proves that is invariant [9,23,31].
This completes the proof. Average number of contact partners q 1 Transmission probability of susceptible individuals