Abstract
We consider a family of periodic SEIRS epidemic models with a fairly general incidence rate of the form , and it is shown that the basic reproduction number determines the global dynamics of the models and it is a threshold parameter for persistence of disease. Numerical simulations are performed using a nonlinear incidence rate to estimate the basic reproduction number and illustrate our analytical findings.
1. Introduction
Epidemiological models in mathematics have been recognized as valuable tools in analyzing the dynamics of an infectious disease nowadays. They are used to describe the spread of disease and also to make control measures known to avoid its persistence, for example, via vaccination terms or treatment terms. These models consider the total population divided into compartments, given by the biological assumptions on the model and represented by functions depending on time . The most common categories used are susceptible (), infected (), recovered (), exposed (), quarantined (), and vaccinated (), and the dynamics of model is given by transmission rates from a compartment to another. We have then indicated that the models could be of type , , , , , , and so forth.
To ensure that the model can give a justified qualitative description of the disease, the choice of the incidence rate plays an important role. An incidence rate is defined as the number of new health related events or cases of a disease in a population exposed to the risk in a given time period. Some examples are the bilinear incidence rate, the saturated incidence rate, or a general incidence rate. The bilinear incidence rate has been repeatedly used by several authors. It is given by , where is the transmission rate and the product represents the contact between infected and susceptible individuals (based on the law of mass action). It was introduced by Kermack and McKendrick [1] in 1927, and even when it is mathematically simple to use, it faces multiple problems and challenges when it is used to describe disease propagation among gregarious animals or persons [2], because it goes to infinity when becomes larger. In order to improve the modelling process to study the dynamics of infection among a large population, Capasso and Serio [3] in 1978 introduced a saturated incidence rate by studying the Cholera epidemic spread in Bari, given by , where is the transmission rate and the saturation constant. Unlike the bilinear incidence, saturated incidence does not grow up without a limit, but it goes to a saturation limit as goes to infinity. Multiple types of saturated incidence have been used in the literature; see, for example, [2] for a list of them. To avoid the use of a single incidence function, the use of a general incidence rate that includes a family of particular functions with similar properties has become a topic of interest by several authors (see, e.g., [4–8]).
The basic reproduction (represented by ) is defined as “the average number of secondary cases produced by a single infected case when it is introduced in a susceptible population” and it has an important role in the study of disease transmission. In biological terms, usually when this number is less than one, the disease is eradicated from population, but when it is greater than one, the infection persists. Mathematically, it is of interest to compute a threshold parameter with the properties of the basic reproduction number. A method to compute this number for certain compartmental disease models is via the next-generation matrix method developed in [9]; however, it is not useful when the model presents time periodic seasonal terms. Authors like [10, 11] have defined its basic reproduction number for periodic models as an average, to give some results about extinction or persistence of infection. However Bacaër and Guernaouni in [12] introduced the definition of basic reproduction number for periodic environments, and, later, Wang and Zhao [13] made a formal definition of it, via the monodromy matrix.
In the present work, we focus on a family of SEIRS epidemic models with a time periodic seasonal term, improving the model of Moneim and Greenhalgh in [14], by introducing an incidence rate with a general function taken from [4] and the references therein.
We propose the following SEIRS model:where is the total population size, with denoting the fractions of population that are susceptible, exposed, infected, and recovered, respectively. is the transmission rate and it is a continuous, positive -periodic function. () is the vaccination rate of all newborn children. is the vaccination rate of all susceptibles in the population and it is a continuous, positive periodic function with period , where is an integer. is the common per capita birth and death rate. , and are the per capita rates of leaving the latent stage, infected stage, and recovered stage, respectively. It is assumed that all parameters are positive constants.
Bai and Zhou in [5] answered some open problems stated in [14], they also showed that their condition is a threshold between persistence and extinction of the disease via the framework established in [13]. They assumed that the incidence was bilinear. In our study, the nonlinear assumptions on function are listed below (see [4]).(A1) is continuously differentiable.(A2), and for all .(A3).
Under these assumptions, function includes various types of incidence rate; in particular, when , we are on the bilinear case considered in [14].
In addition, we assume the following extra conditions (see [15]).(A4).(A5)There exists such that when ,
This set of assumptions on the function allows for more general incidence functions than the bilinear one, like saturated incidence functions and functions of the form ; in particular, in the case when , they represent psychological or media effects depending on the infected population. In this last case the incidence function is nonmonotone on . (A3) regulates the value of comparing it with the value at of a line containing the origin of slope (note that this line varies as increases), (A4) requires a concave at the origin, and (A5) imposes the geometrical condition that in a small neighborhood of the origin must lie between the tangent line of at and a concave parabola tangent to at .
We consider a family of epidemic models with periodic coefficients and general incidence rate in epidemiology. Then we show that the global dynamics of solutions is determined by the basic reproduction number , generalizing the results in [5]. The layout of this paper is as follows: In Section 2, we prove the existence of a disease-free periodic solution and we introduce the basic reproduction number via the theory developed in [12, 13]. In Section 3, we adapt the arguments given in [5] to prove that the disease-free periodic solution of system (1) is globally asymptotically stable if and it is persistent when . Finally, in Section 4, we give some numerical simulations of our results, making a comparison between our basic reproduction number and the average reproduction number used by several authors (see, e.g., [10, 11]).
2. The Basic Reproduction Number
First of all, we prove nonnegativity of the solutions under nonnegative initial conditions.
Theorem 1. Let . The solution of (1) with is nonnegative in the sense that , , and satisfies , with constant.
Proof. Let ; then, adding all equations of system (1), we can see that , so the value of is constant. Now, set as the solution of system (1) under initial conditions . By the continuity of solutions, for all of and that have a positive initial value at , we have the existence of an interval such that for . We will prove that .
If for a and other components of remain nonnegative at , then implying that whenever the solution touches the -axis, the derivative of is nondecreasing and the function does not cross to negative values. Similarly, when for a and other components remain nonnegative, we have When for a and other components remain nonnegative, Finally, when for a and other components remain nonnegative, Therefore, whenever touches any of the axes , , , or , it never crosses them.
In order to make the analysis of the model in a simpler way from now on, we make a reduction of dimension in system (1) making , obtaining the following:
The dynamics of system (1) is equivalent to that of (7); moreover, due to positivity of solutions, we have , so we study the dynamics of system (7) in the region
A disease-free periodic solution can be found for (7). To find it, set ; then, from the first equation of (7) we can obtain the following initial value problem:
From [5, 14], the equation above admits a unique positive -periodic solution given bywhere
Therefore, is a disease-free periodic solution of (7); moreover, from [5] we have that ; therefore, lives in .
Using the notation of [9], we sort the compartments so that the first two compartments correspond to infected individuals. Let and define(i): the rate of new infection in compartment ,(ii): the rate of individuals into compartment by other means,(iii): the rate of individuals transfer out of compartment .
System can be written aswhere ,
Linearizing system (12) around the disease-free solution, we obtain the matrix of partial derivatives , where
Using Lemma of [9], we part and and set
For a compartmental epidemiological model based on an autonomous system, the basic reproduction number is determined by the spectral radius of the next-generation matrix (which is independent of time) [9]. The definition of basic reproduction number for nonautonomous systems has been studied for multiple authors; see, for example, [12, 13]. Particularly, Wang and Zhao in [13] extended the work of [9] to include epidemiological models in periodic environments. They introduced the next infection operator given bywhere is the ordered Banach space of all periodic functions from to , which is equipped with the maximum norm. is the initial distribution of infectious individuals in this periodic environment, and is the evolution operator of the linear periodic system:meaning that, for each , the matrix satisfies
is the distribution of accumulative new infections at time produced by all those infected individuals introduced before , with kernel . The coefficient in row and column represents the expected number of individuals in compartment that one individual in compartment generates at the beginning of an epidemic per unit time at time if it has been in compartment for units of time, with [16].
Let is an eigenvalue of if there is a nonnegative eigenfunction such that
Therefore, the basic reproduction number is defined asthe spectral radius of . The basic reproduction number can be evaluated by several numerical methods and approximations [15–17]; in Section 4 we discuss this topic.
3. The Threshold Dynamics of
3.1. Disease Extinction
Theorem 2. Let be defined as (20); then the disease-free periodic solution is asymptotically stable if and unstable if .
Proof. We use Theorem of [13] and check conditions (A1)–(A7). Conditions (A1)–(A5) are clearly satisfied from the definitions of and given in Section 2. We prove only conditions (A6) and (A7). Define and let be the monodromy matrix of system(A6). Let be a fundamental matrix for system , with defined as before and periodic; the monodromy matrix is given by . The general solution of (22) is so and . Note that , so and Due to the fact that is a constant, its eigenvalue is itself and for .(A7). Solving the system , we arrive at the general solution so Computing , we have Clearly, for .
Note 1. Due to the fact that is a fundamental solution of a periodic system, we can always choose it such that , so the monodromy matrix satisfies . This property is used in further analysis.
In order to prove the global stability of the disease-free periodic solution, we enunciate some useful definitions and some lemmas.
Let be continuous, cooperative, irreducible, and -periodic matrix function, and the fundamental matrix of system . Denote by the spectral radius of .
Lemma 3. Let . Then there exists a positive, -periodic function such that is a solution of (see proof in Lemma 2.1 of [18]).
Lemma 4. Function of model (1) satisfies , .
Proof. Using assumptions on function , we haveso function decreases and then .
Lemma 5. Let be a solution of system (7) with initial conditions , and the disease-free periodic solution of (7); then
Proof. Proof is similar to Lemma of [14]. satisfies the first equation of system (7); then Let ; then Using Gronwall’s inequality , Taking limits in both sides, we obtain that .
Now, we are able to enunciate our theorem for global stability of disease-free periodic solution.
Theorem 6. The disease-free periodic solution of system (7) is globally asymptotically stable if .
Proof. From Theorem 2 we have that is unstable for and asymptotically stable for , so it is sufficient to prove that any solution with nonnegative initial conditions approaches as tends to infinity.
Let ; from Lemma 5 we have so there exists a such that for all which implies that . Then, from the definition of supremum, we have that for all Then, we have proved that for all we can find a such that for all .
Now, using Lemma 4, for we can find a such that for We consider the following perturbed subsystem:which can be rewritten as with defined in (15) andMatrix is -periodic, cooperative, irreducible, and continuous. Using Lemma 3, if , then there exists a positive and -periodic function such that is solution of system (38). Note that for all , function is also a solution of system (38) with initial condition at .
Choose a and such that ; then from (37) we have that and using a comparison principle (see, e.g., [19] Theorem ), we have for all .
From Theorem of [13], iff . By the continuity of the spectrum for matrices (see [20], Section II.5.8), we can choose small enough so that and then (see Note 1). Thus, using positivity of solutions and comparison, And similarly for , we obtain thatWe need only to prove that approaches . At disease-free periodic solution , where satisfies Thus, satisfiesLet be arbitrary and . Due to (43) we can find a such that for ; moreover, we can find a such that for . Then, let ; we have for Multiplying in both sides by and integrating from to , we obtainSo, , where is arbitrarily small. Then , and using similar arguments for and , we can find a with for . Also, from (43), we can find with for , so, for , we haveOr, equivalently, , with being arbitrarily small, and this implies that . We conclude by comparison and using Lemma 5 that , completing the proof.
Theorem 6 shows that disease will completely disappear as long as . Thus, reducing and keeping below the unity would be sufficient to eradicate infection, even in a periodic environment and a general incidence rate.
3.2. Disease Persistence
Uniform persistence is an important concept in population dynamics, since it characterizes the long-term survival of some or all interacting species in an ecosystem [21].
In this section we consider the dynamics of the periodic model when . We will show that actually is a threshold parameter for the extinction and the uniform persistence of the disease. Our results are inspired by [5, 15, 18, 22].
Let be the Poincaré map associated with system (7); that is, where is defined in (8) and is the unique solution of system (7) with . We define the following sets:
Note that is not the boundary of , but it is a standard notation of persistence theory.
Lemma 7. Set is positively invariant under system (7).
Proof. Let , that is, , and letbe the solution of (7) with . Due to nonnegativity of solutions and assumptions on function and , we have Using a comparison theorem (see, e.g., [19] Appendix B.1), we have for all Similarly, so, Therefore, remains on for all .
To use persistence theory developed in [21], we show thatwhere
Let and be the solution that passes through that initial condition. We have that , with being a solution of (9) and being a solution that satisfies the initial condition. By uniqueness of solutions we have , so lives on .
Now, if , we want . We prove an equivalent sentence: if , then it does not belong to . Consider an initial point ; then , or . Suppose that and ; then holds
By continuity of and sign of derivative of , we have that, for small , , so, for , . Using invariance of (Lemma 7) we have for all . Finally, for a such that , we have and this implies (56). By the existence of a disease-free periodic solution (proved in Section 2), it is clear that there is one fixed point of in given by ([23]).
Now, we are in a position to introduce the following result of uniform persistence of the disease.
Theorem 8. Let ; then there exists an such that any solution of (7) with initial values satisfies
Proof. We first prove that is uniformly persistent (see Definition from [21]) with respect to , because this implies that the solution of (7) is uniformly persistent with respect to (see [21], Theorem ). Clearly, is relatively open in , so is relatively closed.
Define we show that
By Theorem of [13], if and only if . Choose an small enough with the property (see Appendix A). For , let us consider the following perturbed equation:System above admits a unique positive -periodic solution of the formwhit , which is globally attractive for all solutions of (61) (see Appendix B for proof), and withSince is continuous in , for all there is a such that for we have . Moreover, by continuity of solutions with respect to initial values we can find for all an such that if , then Therefore, for established before, we can find small enough such that , .
Again, by continuity of solutions with respect to initial values, for this small , there exists a such that for all with we have , .
We now claim thatBy contradiction, suppose thatWithout loss of generality, we can assume that for all (see Appendix C). From the discussion above, , and .
For any , let , where and is the greatest integer less than or equal to . Then, we get If we set , then we have , , , and from the first equation of (7) and Lemma 4 we arrive atwhich is exactly the equation in (61). Since the unique periodic solution of (61) is globally attractive, we have for solution of (61) that . So for given before, there exists such that for all or equivalently . Moreover, from previous analysis, ; therefore, using comparison principle on (68) we arrive atfor .
We have , and is fixed small, so we can take and use assumption (A5) in Introduction (see Appendix D) to obtainwhere are defined in (15), is defined in (40), and By Theorem of [13], we have iff . By continuity of spectrum (see [20] Section II), we can find such that Consider the auxiliary system then, using Lemma 3 there exists a solution of (71) with the form , with . Choose a and a small number such that . Using comparison principle we get , which implies and . This leads to a contradiction.
The claim above shows that is weakly uniformly persistent with respect to Note that has a global attractor (see Lemma 5). It follows that is an isolated invariant set in , . Every orbit in converges to and is acyclic. By the acyclicity theorem on uniform persistence for maps ([21] Theorem and Remark ), it follows that is uniformly persistent with respect to ; that is, there exists such that any solution of (7) satisfies ,
4. Numerical Simulations
In this section we provide some numerical simulations to illustrate the results obtained in our theorems and compare them with previous results.
To improve previous models used in references, we use a particular functionwhich includes the case used in [5]. One can check that function (75) satisfies conditions (A1)–(A5). Using this function, system (7) is rewritten as
Set an initial population and take time in years. Suppose per year, corresponding to an average human life time of 50 years. Following [5] take the parameters as follows: per year, per year, , and . Choose the periodic transmission as , with being the transmission parameter, and the periodic vaccination rate . Both functions have period .
There exists multiple methods for computing the basic reproduction number, via numerical approximations, or finding a positive solution of the equation (see Theorem of [13]). In order to compare our work with previous works, we approximate the basic reproduction number with its average value , used by several authors as the reproduction number (for example [10, 11]), so definewhere is given by (15) and with being the average of functions defined as . Computing each average, we obtain so for
Following Theorem of [13], to compute , let , be the evolution operator of the systemthat is, for each , , , and . With this operator, is the unique solution of
Example 1. To illustrate our results, fix . Computing , we have , which is a first approximation of . To solve system (80) numerically, we substitute the terms of expression of in (10): The previous integral cannot be computed analytically, so we approach using Taylor expansion around 0 (remember that we want so solve , where ), so even when we cannot find an explicit expression for , the Taylor expansion is a good way to estimate it in . It could be of interest to also use an approach of around and compare the results with those obtained in the present work (see Section 5 for a discussion about this topic).
Setting an initial value and letting , we solve system (80) numerically for each (using initial conditions and , to satisfy ), and compute until . With previous process we arrive at for and for ; therefore . Using a finer step size to have more accuracy, we arrive at .
Set initial values as , , , and .
There exist multiple numerical methods to compute and plot the solutions of nonautonomous differential equations; see, for example, the Adomian method, the homotopy analysis method, or the modified homotopy methods (see, e.g., [24, 25]). For this work we use Matlab algorithms (ODE 45) to graph the solution of system (76) with these initial conditions. Figures 2 and 3 shows the results. We can see that goes to zero, while tend to stabilize; also is tending to with values between 54,000 and 56,000 (see Figure 1); this shows the results obtained in Theorem 6.
(a)
(b)
Example 2. Now, choose . As we can see in Figures 4 and 5, the solutions of system (1) remain persistent when tends to infinity; this fact suggests that from Theorem 8. In fact, if we compute the basic reproduction number and its average (using the process described in example 1), and ; therefore it is bigger than one. In fact, this shows the results of persistence obtained in Theorem 8.
5. Conclusion
In this paper we presented a model with seasonal fluctuation with a general incidence function that includes the bilinear case (studied by [5]) and a family of saturated incidence rate of the form . We proved the existence of a disease-free periodic solution and defined the basic reproduction number , proving that it is a threshold parameter for disease, in the sense that when , the disease-free periodic solution is globally asymptotically stable, and when , the disease is persistent. A next step of this work is to consider a family of incidence rates more generally, changing by and trying to obtain results of persistence and stability similar to the ones obtained in this work. Another interesting topic is to ask what the behavior of system at is, in order to complete the analysis that we have made.
Several authors (e.g., [10, 11]) define as an average, which we denoted as to distinguish between it and the basic reproduction number defined by [13], via the monodromy matrix (which is a real threshold parameter for extinction and persistence of disease). We compute , approximate (with the help of Taylor theorem), and compare these values, obtaining that is not equal to ; moreover in both examples (similar comparisons can be observed also in the works made by [13, 17]). This fact suggests that the use of for persistence overestimates the threshold. To emphasize this conclusion, it would be helpful to find an example where but and then compute the solutions to observe the behavior (we affirm that the disease will go extinct due to Theorem 6).
To obtain the estimation of we used a code in Maple, which is based on numerical computing of until , where , is the step size, and the initial estimation is taken as . For this approximation we have used a Taylor expansion of the periodic solution ; another interesting possibility could be varying the approximation used for , for example, changing the Taylor approach of around instead of . The graphs of the solutions were obtained with ODE 45 from Matlab, but other methods can be used to improve them, for example, Adomian methods or homotopy methods [24, 25]. The Maple code used to estimate is available for anyone who wants to use it.
Appendix
A. Assumption on Used in Theorem 8
Note that has a positive minimum value (it is periodic, positive, and continuous, so it is bounded for and then for all ) and we can choose a with the property , sufficiently small such that .
B. Periodic Solution of (61)
For each , (61) used in the proof of Theorem 8 is
Solving the equation above, we arrive at the general solutionwhere . We shall examine the behavior of an arbitrary solution . For each , we can use an initial time with initial point and see that
Since is a periodic function, where . ThenAnd using the change of variable , we have
Equation (B.6) gives a recursive relationship between the solution at and after times. If we set , then for each solution this relationship is described by with being on the right side of (B.6). If we take and , two different values of , then
Then, is a contracting map, and by Banach fixed point theorem has a unique fixed point such that or, equivalently, . This fixed point can be found for any that is a solution of a differential equation with arbitrary initial condition at any time . The fixed point has the form Thus, define the function
is a periodic function with period and is continuously differentiable with respect to . One can check (by computing the derivative) that is a solution of differential equation, so by existence and uniqueness of solutions it can be rewritten as with initial condition
If we suppose the existence of another periodic solution , then using (B.6) we arrive at , by uniqueness of solutions , and the periodic solution is unique. Computing the difference , we have so, . Therefore, every solution converges to .
C. Assumption on Used in Theorem 8
Let . If then we have . For all there exists a such that if , then . In particular, for we have or, equivalently, for . Moreover, for all with , we have . Therefore, , .
We can take as initial condition and, therefore, making our assumption valid.
So, we can assume without loss of generality that for all .
D. Expression (71)
From system (7) , with for , so Using assumption (A5) for and positivity of , we have also Therefore, and , so and ; applying this we arrive at This expression can be written as (71).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This article was supported in part by Mexican SNI under Grants nos. 15284 and 33365.