Condition for Global Stability for a SEIR Model Incorporating Exogenous Reinfection and Primary Infection Mechanisms

A mathematical model incorporating exogenous reinfection and primary progression infection processes is proposed. Global stability is examined using the geometric approach which involves the generalization of Poincare-Bendixson criterion for systems of n-ordinary differential equations. Analytical results show that for a Susceptible-Exposed-Infective-Recovered (SEIR) model incorporating exogenous reinfection and primary progression infection mechanisms, an additional condition is required to fulfill the Bendixson criterion for global stability. That is, the model is globally asymptotically stable whenever a parameter accounting for exogenous reinfection is less than the ratio of background mortality to effective contact rate. Numerical simulations are also presented to support theoretical findings.


Introduction
Mathematical models, in particular, models tracking dynamics of infectious diseases, are of utmost importance due to their application in the assessment of public health policies by national and international agencies. In such models, one of the intriguing aspects that often occurs is when modellers need to know whether the disease will disappear or will permanently remain in the population. This question is answered mainly through investigating the asymptotic stability of the disease-free equilibrium (DFE) as well the endemic equilibrium. It is already known that if the DFE is globally asymptotically stable, then the disease eradication is assured regardless of the initial number of infected individuals in the population [1]. An influx of infected cases may trigger an isolated epidemic outbreak, but they may not make the disease endemic in the population [2]. In contradiction to DFE, if the endemic equilibrium is globally asymptotically stable (GAS) and a few infected individuals are initially introduced, then, the disease will be permanently present in the population.
In the sequel, there are two major methods used in the analysis of global stability of endemic equilibrium: Lyapunov direct method and geometric method. Although the Lyapunov direct method is often used in proving global stability of infectious diseases models, it is sometimes difficult to use because it requires an auxiliary function which is hard to construct. This is because there are no existing general methods for constructing such Lyapunov functions. Moreover, Lyapunov functions for models with parameter(s) that induce bistability phenomena may not even exist. The second method, sometimes referred to as geometric approach to global stability, is a generalization of the Poincare-Bendixson criterion for systems of n ordinary differential equations. This method was developed by Li and Muldowney [3,4] in midnineties to address problems encountered with the Lyapunov direct method. The technique is extensively being used to analyze global properties of mathematical models emanating in mathematical epidemiology as well as in other biomathematical contexts. For instance, its applications can be seen in toxicant-population interaction models, Lotka-Volterra models incorporating delay [5,6], and mathematical epidemic models modelling dynamics of HIV in a human host [7]. Due to mathematical technicalities involved with the geometric approach, it seems to work quite well with low-dimensional mathematical models such as SIR (Susceptible-Infective-Recovered) and SIRS (Susceptible-Infective-Recovered-Susceptible) which can be reduced to bidimensional models (i.e., n = 2). Moreover, this method has also been found to be more appropriate in SEIR-like (Susceptible-Exposed-Infective-Recovered) which are represented by a system of four ordinary differential equations (ODEs). This is because their dynamics can be reduced to a three-dimensional system (see [3,[7][8][9][10][11]).
Li and Muldowney [3] proved global stability of a SEIRlike model. However, there are some special cases of SEIRlike models that may not be entirely explained by the global thresholds obtained in [3]. Thus, building on the work of Li and Muldowney, we consider new dynamics that may alter the global stability conditions. For, instance if individuals in incubation period which in this case in compartment E experience reinfection, then, new dynamics are likely to arise which may alter the global conditions. Several diseases such as SARS (Severe acute respiratory syndrome), HIV (Human immunodeficiency virus), Ebola, malaria, influenza, and tuberculosis (TB) have an incubation period, and during this stage of infection, individuals can be reinfected by the same or a new strain of the disease (exogenous reinfection).
Castillo-Chavez and Song [12] reviewed the earliest mathematical models of TB dynamics that appeared in the 1960s. Their review investigated several epidemiological factors that are pertinent to TB transmission dynamics. These included the role played by close and casual contacts on TB dynamics, the role of demographic factors (births and deaths), and the role played by intervention strategies implemented to mitigate TB. Further, they extended their review by considering other models such as cell-based models for TB transmission at the immune system level as well as Markov chain models primarily focusing on TB projection. Although these authors considered a SEIT (Susceptible-Exposed-Infective-Treated) model to analyze the local stability of the disease-free equilibrium in the neighbourhood of R 0 = 1, using the famous Castillo-Chavez method derived in [12], they never considered the global properties of the persistent equilibrium. Again, the SEIT proposed model did not incorporate the primary progression pathway in the analysis of the local stability using the Castillo-Chavez method.
It is imperative to note that the Castillo-Chavez et al. [12] method for local stability is used to verify whether a dynamical system exhibits the phenomenon of backward bifurcation. The SEIT model reviewed in [12] obtained a condition for backward bifurcation which is different from the one obtained in the model studied here simply because they omitted the primary progression pathway. The exclusion of primary progression pathway and lack of investigating the global properties of the persistent equilibrium motivated us to attempt to understand the global properties of a SEIR model. Moreover, the model studied here goes beyond understanding the phenomenon of backward bifurcation which is well documented in literature (see [13][14][15][16][17][18] and the references therein); instead, we attempt to investigate the global properties of a SEIR model incorporating both exogenous reinfection and primary progression infection pathways which was not studied in [12]. The global stability analysis of the persistent equilibrium performed on the proposed SEIR model using the geometric approach is the central feature that distinguishes our model from Castillo-Chavez and Song's [12] paper.
Many infectious diseases in particular those caused by bacterial and viral infections do not render permanent immunity after recovering from the first episode. As a result, they are characterized by partial or complete loss of immunity and subsequent exogenous reinfection. Exogenous reinfection can occur in diseases such as tuberculosis (TB) and flu. Often upon initial infection, individuals pass through different stages before they reach the symptomatic stage where they manifest disease symptoms. Some individuals may pass through a latent stage or dormant state (as in TB or  where the disease stays inactive for a specified period of time while some may progress directly to the infectious stage (fast progression or primary progression). An individual infected with TB can either pass through a latent stage or progress directly to the infectious stage (primary progression).
The novelty of this study is the incorporation of exogenous reinfection and primary progression infection mechanisms which were excluded in previous threedimensional SEIR-like models used to study global stability (see [3,7,11,19]). Thus, our main task is to investigate global properties of a SEIR model including exogenous reinfection and primary progression infection processes using geometric approach. According to our knowledge, the global stability results deduced in this model have not been previously obtained.

Mathematical Model
The host population is partitioned into four compartments, the susceptible, exposed (latent), infective, and recovered, with subpopulations denoted by S, E, I, and R, respectively. The total population at time t is represented by NðtÞ = SðtÞ + EðtÞ + IðtÞ + RðtÞ. We assume a frequency-dependent incidence rate (standard incidence rate) as the force of infection. Here, c represents the host-host contact, β is the probability that a contact results in transmission, and IðtÞ/NðtÞ is the prevalence of infection, sometimes referred to as the "frequency" of infection [20]. After initial infection, a proportion q of susceptible individuals can move directly to the infectious stage (primary progression) while the rest move to latent compartment (slow progression), while in latent compartment, individuals are assumed to be exogenously reinfected, hence moving to the infectious stage at a rate pλE. Precisely, the subpopulation of susceptibles is generated by recruitment through births and immigration at a rate Λ. The population is decreased due to contact with infectious individuals at a rate λS. The exposed 2 Computational and Mathematical Methods in Medicine subpopulation is generated due to slow progression of the disease at a rate ð1 − qÞλS. Moreover, individuals leave the exposed subpopulation due to exogenous reinfection and endogenous reactivation at rates pλE and kE, respectively. The infected subpopulation is generated by both fast progression and exogenous reinfection at rates qλE and pλE: It is diminished when individuals recover from disease due to therapy at rate γI and disease-related death rate μ d . The recovered subpopulation is generated by recovery of infected individuals at a rate γ. All individuals in each subpopulation experience natural death at the same rate through the background mortality rate μ: The equation representing transfer of individuals among the compartments is now given by

Existence of the Equilibria
For mathematical tractability, we transform model equation (2) in terms of proportions of the individuals in each compartment. Thus, we define s = S/N, e = E/N, i = I/N, and r = R/N as proportions for the compartments S, E, I, and R, respectively. Now differentiating the corresponding proportions with respect to time, it is not difficult to deduce that s, e, i, r, and N satisfy the following differential equations: where Ω = fðs, e, i, r, NÞ ∈ ℝ 5 It is easy to observe that the variable r in equation (3) does not affect other equations; thus, we can drop the fourth equation. Again, all equations in (3) depend on N. So replacing Λ/N = μ + μ d i into the first, second, and third equations, we have the following reduced system: It is easy to show that the region Ω = fðs, e, iÞ ∈ ℝ 3 (4) is biologically feasible is positively invariant, where ℝ 3 + represents the nonnegative cone of ℝ 3 . In the absence of disease, the model equation (4) has an intrinsic equilibrium point P 0 = ð1, 0, 0Þ. This equilibrium is the disease-free equilibrium point (DFE).
The basic reproduction number is defined as the number of secondary infections generated by one newly infected individual when introduced into an entirely susceptible population at the disease-free equilibrium point during its mean infective period [1]. The basic reproduction number denoted by R 0 is calculated following the next-generation operator method as developed in Van den Driessche and Watmough [1]. Note that the nonlinear term with new infections F and the transition term V are, respectively, given by The linearized matrices F and V computed at the diseasefree equilibrium P 0 yield Hence, the next generation matrix is One of the eigenvalue of the above matrix is zero while the other (dominant eigenvalue) gives the basic reproduction number. Thus, the basic reproduction number is given as It is easy to see that R 0 does not contain parameter p which accounts for exogenous reinfection. This indicates that 3 Computational and Mathematical Methods in Medicine R 0 alone may not sufficiently explain dynamics of model system (4); rather, additional restrictions are needed. Now, we proceed to obtain the nonzero steady states when disease is present in the population. Setting the equations of system (4) to zero and simplifying lead to the following expressions written in terms of i * Substituting expression (9) into the third equation of system (4) leads to the following: where The root i * = 0 corresponds to a scenario where there is no disease in the population (DFE). The other nonzero roots of model system (4) can be obtained by solving for i * from and substituting them in (9). From equation (12), it is not difficult to see that d 2 is always positive, d 0 > 0 ⇔ R 0 < 1, and d 0 < 0 ⇔ R 0 > 1: Now, by Descartes Rule of Signs, it follows that there is a unique endemic equilibrium whenever d 0 < 0: two positive endemic equilibria if d 0 > 0, d 1 < 0, d 2 1 − 4d 0 d 2 > 0 and no positive endemic equilibria otherwise. Moreover, a bifurcation point occurs when d 0 > 0, d 1 < 0, and d 2 1 − 4d 0 d 2 = 0 (i.e., the point where the two positive endemic equilibria collide, leaving the DFE as the only equilibrium point). At the point where the two positive equilibria (which we shall denote by R c 0 ) collide, The quantity R c 0 given in the appendix equation (A.2) can be easily obtained by setting d 2 1 − 4d 0 d 2 = 0 (see the appendix for detailed derivation). Now taking into account the above discussion on equation (12), we deduce Lemma 1.

Lemma 1.
(i) For d 1 > 0, then, the model equation (4) admits no positive real equilibria whenever R 0 < R c 0 < 1 (ii) For d 1 < 0, then, the model equation (4) admits two positive endemic equilibria P 1 and P 2 whenever R c 0 < R 0 < 1 (iii) For d 1 < 0, then, the model equation (4) admits a unique positive endemic equilibrium P * whenever R 0 ≥ 1 The occurrence of two positive endemic equilibria when R 0 < 1 hints the possibility of backward bifurcation phenomenon. Hence, the next section of the work is dedicated on detemining the type of bifurcation exhibited by model equation (4).

Proof of Existence of Backward Bifurcation
Define Then, Theorem 2 follows:

Theorem 2.
(i) The model equation (4) exhibits backward bifurcation at R 0 = 1 whenever p > p c (ii) The model equation (4) exhibits forward bifurcation at R 0 = 1 whenever p < p c Proof. To show that model system (4) exhibits backward bifurcation phenomena, we apply the Center Manifold approach as outlined by Castillo-Chavez and Song in [12]. For clarity and understanding of the center manifold theory, the model equation (4) variables are transformed as follows: y 1 = s, y 2 = e, y 3 = i, and the total population n = ∑ 3 j=1 y j . Define Y = ðy 1 , y 2 , y 3 Þ T (T denotes transpose), such that the model equation (4) can be rewritten as dY/dt = FðyÞ, where F = ðf 1 , f 2 , f 3 Þ. Hence, it follows Now, let cβ =β and chooseβ as the bifurcation parameter. Note that at R 0 = 1, Then, the Jacobian matrix of equation (14) evaluated at DFE is given as Computational and Mathematical Methods in Medicine Withβ = β * , the transformed system (14) has a hyperbolic equilibrium point (i.e., it has a simple eigenvalue with zero real part, and all other eigenvalues are negative). Hence, we can apply the center manifold theory in [12] to analyze dynamics of the transformed system (14) nearβ = β * . It is easy to obtain the right and left eigenvectors of the Jacobian matrix J ðP 0 Þ , respectively, denoted by v = ðv 1 , v 2 , v 3 Þ T and w = ðw 1 , w 2 , w 3 Þ.

Right eigenvectors
Now, we proceed to obtain the associated bifurcation parameters, a and b, as described in Theorem 4.1 of [12], where To obtain the associated bifurcation coefficient a, we first obtain the nonvanishing partial derivatives of model system (14) evaluated at DFE. Hence, it follows that so that where Moreover, the nonvanishing partial derivatives associated with b are The stability of the model system (4) switches at the transcritical point R 0 = 1: According to Theorem 4.1 of [12], it is stated that if both bifurcation coefficients a and b are positive, then the system exhibits backward bifurcation. Notice that a > 0 if and only if p > p c implying that if this condition hold then the model equation (4) will enter into a bistability regime where there is coexistence of two positive endemic equilibria when R 0 < 1, (i.e., backward bifurcation). If p < p c , then the model will have forward bifurcation at R 0 = 1.
The type of bifurcation that occurs at R 0 = 1 is primarily determined by exogenous reinfection parameter p whenever it exceeds a certain critical value denoted by p c . Thus, the disease can persist even though the reproduction number is less than one. In such case, the disease can only be eliminated if the basic reproduction number is decreased below a certain threshold. That is, R 0 < R c 0 < 1. Figures 1(a) and 1(b) show the associated forward bifurcation and backward bifurcation, respectively. It is clear that incorporation of exogenous reinfection does induce new dynamics in model system (4) when R 0 < 1.
It is imperative to mention that in the presence of exogenous reinfection, global stability of the endemic equilibrium P * when R 0 > 1 may not be automatically guaranteed. The next section of the work is dedicated on finding under what condition the endemic equilibrium point P * is globally asymptotically stable whenever R 0 > 1:

Global Stability Using the Poincare-Bendixson Property
Since our objective is to establish global stability of the unique endemic equilibrium P * when R 0 > 1 in the presence of exogenous reinfection parameter p, we first briefly outline the general mathematical framework of the procedure as developed in M. Li and J. Muldowney [3,19].

Computational and Mathematical Methods in Medicine
Suppose the map y ↦ f ðyÞ is a C 1 function for y in an open subset D ⊂ ℝ n , and consider the following autonomous dynamical system: Let yðt, y 0 Þ be the solution to equation (24) satisfying yð0, y 0 Þ = y 0 . Now, we make the following basic assumptions: (H1) D is simply connected.
(H2) There exists a compact absorbing set K ⊂ D.
(H3) Equation (24) has a unique equilibrium y * in D. Now under the stated assumptions (H1)-(H3), y * is said to be globally stable in D if it is locally stable and all trajectories in D converge to the same equilibrium y * . That is, system (24) has no nonconstant periodic solutions. It is important to mention that the major role for global stability is determined by the Bendixson criteria. For n ≥ 2, a Bendixson criterion refers to a condition satisfied by field f which precludes the existence of nonconstance periodic solutions of equation (24). When n = 2 (i.e., the planar case), the classical results (Poincare-Bendixson theorem and Dulac criteria; see [21]) adequately provide such global conditions. For n ≥ 3, a remarkable approach for proving global stability can be traced in the work due to Li and Muldowney [3,4,19]. In their paper, they showed that if conditions (H1)-(H3) hold and differential equation (24) fulfills a Bendixson criterion that is robust under C 1 local ε − perturbations of f at all nonequilibrium nonwandering points for system (24), then y * is globally stable in D provided it is stable. Note that a function g ∈ C 1 ðD → ℝ n Þ is called a C 1 local ε − perturbation of f at y 0 ∈ D if there exists an open neighbourhood U of y 0 in D such that the support sup ð f − gÞ ⊂ U and jf − gj C 1 < ε, where j f − gj C 1 = sup f|f ðyÞ − gðyÞ|+|f y ðyÞ − g y ðyÞ| : y ∈ Dg. Also, a point y 0 ∈ D is said to be nonwandering for system (24) if for any neighbourhood U of y 0 in D and there exists arbitrary large t such that U ∩ yðt, UÞ ≠ ∅. An example is any equilibrium, alpha limit point, or omega limit point which is nonwandering. We now state the new Bendixson criterion robust under C 1 local ε − perturbations and based on the use of the Lozinski measure as developed in [19]. Now, consider the differential equation (24) under the stated assumptions (H1)-(H3). Let PðyÞ be a n 2 × n 2 matrix-valued function which is C 1 for y ∈ D, and consider where P f is the directional derivative of P in the direction of the vector field f in system (24) and is defined as and J ½2 represents the second additive compound matrix J (that is, Df ðyÞ = JðyÞ). In [22] where relation of compound matrices to differential equations is established, it is shown that for an arbitrary n × n matrix J = J i,j , J ½2 is a n 2 × n 2 matrix. For a special case n = 3, the second additive compound matrix J ½2 can be obtained as   Computational and Mathematical Methods in Medicine Now consider the following quantity q 2 given as where ρðAÞ is the Lozinski measure of A with respect to vector , and is defined as (see [23,24]). In the paper [19], they proved that if conditions (H1) and (H2) are satisfied, then q 2 < 0, indicating that there are no orbits giving rise to simple closed rectifiable curve in D that is invariant for system (24), (that is, periodic orbits, homoclinic orbits, and heteroclinic cycles). Furthermore, it has been remarked in [19] that under the stated assumptions (H1)-(H3), quantity q 2 < 0 implies the local stability of equilibrium point y * . As a result, Theorem 2 is true.
Theorem 3 (see [19]). Assuming that conditions (H1)-(H3) hold, then, the equilibrium point y * is globally asymptotically stable in D if a function PðyÞ and a Lozinski measure ρ exist such that quantity q 2 < 0.
It is straightforward to see that whenever R 0 > 1, there exist unique and positive endemic equilibria P * (see Lemma 1) for model system (4). The method outlined above requires that (i) the endemic equilibrium P * is unique in the interior of Ω (i.e., condition H3 holds) and (ii) in the interior of Ω there exists an absorbing compact set (condition H2 holds). The model equation (4) studied here with the assumption that R 0 > 1 fulfills conditions H1-H3. It is easy to prove that when R 0 > 1, the disease-free equilibrium P 0 is unstable (see [1]). The instability of the disease-free equilibrium P 0 combined with P 0 ∈ δΩ signals uniform persistence [25]. That is, there exists a positive constant k 0 > 0 such that for every solution ðsðtÞ, eðtÞ, iðtÞ, rðtÞÞ of system (4) with ðsð0Þ, eð0Þ, ið0Þ, rð0ÞÞ in the interior of biologically feasible region, Ω satisfies Because of boundedness of the region Ω, uniform persistence is equivalent to the existence of a compact set in the interior of Ω which is absorbing for (4) (see [26]). Hence, condition H1 is satisfied. Also, it is shown that whenever R 0 > 1 the model system (4) has only one equilibrium P * in the interior of Ω, so that condition H3 is verified. Theorem 4. Supposing the conditions R 0 > 1 and 0 < p < μ/cβ are satisfied, then the unique endemic equilibrium P * corresponding to the differential equation (4) is globally asymptotically stable with respect to solutions of (4) originating in the interior of Ω.
Proof. Now, for our model system (4), the task involves showing that quantity q 2 is less than zero or establishing condi-tion(s) that may lead to q 2 < 0. In the interior of the biologically feasible region Ω, suppose conditions (H1)-(H3) hold and let y = ðs, e, iÞ. Let f ðyÞ be the vector field of system (4). The Jacobian matrix J = ∂f /∂y associated with a general solution yðtÞ of our system (4) can be obtained as The second additive compound matrix J ½2 of J is given as where Now based on the model system (4), we choose a suitable vector norm |· | in ℝ 3 and a 3 × 3 matrix-valued function PðyÞ. We set P as : ð34Þ It follows that

Computational and Mathematical Methods in Medicine
The matrix A = P f P −1 + PJ ½2 P −1 is thus obtained as Now, matrix A can be rewritten in block form as where Following [19], we let ðu, v, wÞ represent the vectors in and let ρ represent the Lozinskii measure with respect to this norm. Applying the method of approximating the ρðAÞ as given in [24], we have where g 1 = ρ 1 A 11 ð Þ+ A 12 j j, Note that in equation (41), jA 12 j and jA 21 j are operator norms of A 12 and A 21 with respect to the l 1 vector norm when they are regarded as mapping from ℝ 2 to ℝ and ℝ 2 to ℝ, respectively. ρ 1 ðA 22 Þ represents the Lozinskii measure of the 2 × 2 matrix A 22 with respect to the l 1 norm in ℝ 2 . To obtain ρ 1 ðA 22 Þ, we sum the absolute value of the off-diagonal elements to the diagonal one in each column of A 22 and then take the maximum of two sums. Now, it follows that , Now we have Note that in system (4), second and third equations can, respectively, be written as Now substituting equality (44) into equation (43) and equality (44) into equation (43) leads to Hence, we now have Integrating both sides of equation (47) at the same time for t > t, we get and the Bendixson criterion given by equation (28) is verified. Observe that q 2 < 0 under condition that 0 < p < μ/cβ: Thus, the proof of Thoerem (3) is complete.

Numerical Analysis
Numerically, it is possible to verify the validity of Theorem 4. Recalling that for model equation (4) to be globally stable, the condition 0 < p < μ/cβ must be fulfilled. Now, we numerically investigate two scenarios numerically: that is, where condition p < μ/cβ and where p > μ/cβ.

Conclusion
A mathematical model is proposed, and its global stability property is analyzed using the geometric approach. The model incorporated exogenous reinfection and primary progression infection processes which were excluded in previous SEIR-like models used to investigate global stability through the Bendixson criterion. The analysis of endemic equilibrium points reveals that the model exhibits backward bifurcation phenomena where two positive endemic equilibria coexist when basic reproduction number is below one.
To investigate global stability, we applied geometric approach and found that the Bendixson criterion cannot be satisfied unless the exogenous reinfection parameter p is less than the ratio of background mortality rate (μ) to effective contact rate ðcβÞ. Interestingly, the primary progression parameter q does not impact the condition for global stability for a disease model that follows Susceptible-Exposed-Infective-Recovered stages. Numerically, we verified the validity of global asymptotic stability condition in Theorem 4 and found that it does act as the stability condition for the proposed SEIR model with exogenous reinfection and primary progression infection processes.