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.

1. 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–11]).

Li and Muldowney [3] proved global stability of a SEIR-like model. However, there are some special cases of SEIR-like 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 R0 = 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-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 COVID-19) 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 three-dimensional 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.

2. 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 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 μ. 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

3. 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 Ω = {(s, e, i, r, N) ∈ ℝ+5 : 0 ≤ s, 0 ≤ e, 0 ≤ i, 0 ≤ r, s + e + i + r ≤ 1, N ≤ Λ/μ}.

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 = μ + μi into the first, second, and third equations, we have the following reduced system:

It is easy to show that the region Ω = {(s, e, i) ∈ ℝ+3 : 0 ≤ s, 0 ≤ e, 0 ≤ is + e + i ≤ 1}, where the model equation (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 P0 = (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 R0 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 ℱ and the transition term 𝒱 are, respectively, given by

The linearized matrices F and V computed at the disease-free equilibrium P0 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 R0 does not contain parameter p which accounts for exogenous reinfection. This indicates that R0 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 d2 is always positive, d0 > 0⇔R0 < 1, and d0 < 0⇔R0 > 1. Now, by Descartes Rule of Signs, it follows that there is a unique endemic equilibrium whenever d0 < 0: two positive endemic equilibria if d0 > 0, d1 < 0, d12 − 4d0d2 > 0 and no positive endemic equilibria otherwise. Moreover, a bifurcation point occurs when d0 > 0, d1 < 0, and d12 − 4d0d2 = 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 R0) collide, R0 = R0. The quantity R0 given in the appendix equation (A.2) can be easily obtained by setting d12 − 4d0d2 = 0 (see the appendix for detailed derivation). Now taking into account the above discussion on equation (12), we deduce Lemma 1.

Lemma 1 .

For d1 > 0, then, the model equation (4) admits no positive real equilibria whenever R0 < R0 < 1

For d1 < 0, then, the model equation (4) admits two positive endemic equilibria P1 and P2 whenever R0 < R0 < 1

For d1 < 0, then, the model equation (4) admits a unique positive endemic equilibrium P∗ whenever R0 ≥ 1

The occurrence of two positive endemic equilibria when R0 < 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).

4. Proof of Existence of Backward Bifurcation

Define

Then, Theorem 2 follows:

Theorem 1 .

The model equation (4) exhibits backward bifurcation at R0 = 1 whenever p > p

The model equation (4) exhibits forward bifurcation at R0 = 1 whenever p < p

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: y1 = s, y2 = e, y3 = i, and the total population n = ∑3y. Define Y = (y1, y2, y3) (T denotes transpose), such that the model equation (4) can be rewritten as dY/dt = F(y), where F = (f1, f2, f3). Hence, it follows

Now, let and choose as the bifurcation parameter. Note that at R0 = 1,

Then, the Jacobian matrix of equation (14) evaluated at DFE is given as

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(, respectively, denoted by v = (v1, v2, v3) and w = (w1, w2, w3).

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

so that

The stability of the model system (4) switches at the transcritical point R0 = 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 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 R0 < 1, (i.e., backward bifurcation). If p < p, then the model will have forward bifurcation at R0 = 1.

The type of bifurcation that occurs at R0 = 1 is primarily determined by exogenous reinfection parameter p whenever it exceeds a certain critical value denoted by p. 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, R0 < R0 < 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 R0 < 1.

Figure 1

Illustration of type of bifurcation. (a) Represents forward bifurcation. Parameters used include μ = 0.016, μ = 0.1, k = 0.001, γ = 2, q = 0.05, c = 45, β ∈ {0.42,0.45}, and p = 0.1245 < p = 0.1252. (b) Represents backward bifurcation. Parameters used are the same as in (a) except p = 0.15 > p = 0.1252. β∗ = 0.444 corresponds to R0 = 1. In both figures, the red solid line represents unstable equilibria while the blue solid line represents stable equilibria.

It is imperative to mention that in the presence of exogenous reinfection, global stability of the endemic equilibrium P∗ when R0 > 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 R0 > 1.

5. Global Stability Using the Poincare-Bendixson Property

Since our objective is to establish global stability of the unique endemic equilibrium P∗ when R0 > 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].

Suppose the map y ↦ f(y) is a C1 function for y in an open subset D ⊂ ℝ, and consider the following autonomous dynamical system:

Let y(t, y0) be the solution to equation (24) satisfying y(0, y0) = y0. 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 C1 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 ∈ C1(D → ℝ) is called a C1 local ε− perturbation of f at y0 ∈ D if there exists an open neighbourhood U of y0 in D such that the support sup(f − g) ⊂ U and |f − g| < ε, where |f − g| = sup{∣f(y) − g(y)∣+∣f(y) − g(y)∣:y ∈ D}. Also, a point y0 ∈ D is said to be nonwandering for system (24) if for any neighbourhood U of y0 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 C1 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 matrix-valued function which is C1 for y ∈ D, and consider where P 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, J[2] is a matrix. For a special case n = 3, the second additive compound matrix J[2] can be obtained as

Now consider the following quantity given as where ρ(A) is the Lozinski measure of A with respect to vector norm ∣·∣ in ℝ, , and is defined as (see [23, 24]). In the paper [19], they proved that if conditions (H1) and (H2) are satisfied, then , 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 implies the local stability of equilibrium point y∗. As a result, Theorem 2 is true.

Theorem 2 (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 .

It is straightforward to see that whenever R0 > 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 R0 > 1 fulfills conditions H1–H3. It is easy to prove that when R0 > 1, the disease-free equilibrium P0 is unstable (see [1]). The instability of the disease-free equilibrium P0 combined with P0 ∈ δΩ signals uniform persistence [25]. That is, there exists a positive constant k0 > 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 R0 > 1 the model system (4) has only one equilibrium P∗ in the interior of Ω, so that condition H3 is verified.

Theorem 3 .

Supposing the conditions R0 > 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 Ω.

Now, for our model system (4), the task involves showing that quantity is less than zero or establishing condition(s) that may lead to . 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

It follows that

The matrix A = PP−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 . Now, for the norm |·| in ℝ3, select 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

Note that in equation (41), |A12| and |A21| are operator norms of A12 and A21 with respect to the l1 vector norm when they are regarded as mapping from ℝ2 to ℝ and ℝ2 to ℝ, respectively. ρ1(A22) represents the Lozinskii measure of the 2 × 2 matrix A22 with respect to the l1 norm in ℝ2. To obtain ρ1(A22), we sum the absolute value of the off-diagonal elements to the diagonal one in each column of A22 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

Thus, we have

Hence, we now have

Integrating both sides of equation (47) at the same time for , we get and the Bendixson criterion given by equation (28) is verified. Observe that under condition that 0 < p < μ/cβ. Thus, the proof of Thoerem (3) is complete.

6. 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β.

By choosing the following set of hypothetical parameters, μ = 0.017, μ = 0.1, k = 0.001, γ = 2, q = 0.05, p = 0.0005, c = 60, β = 0.519 and initial conditions s = 1, e = 0, i = 0.00001, with these parameter values, condition p < μ/cβ is satisfied. That is, p = 0.0005 < μ/cβ = 0.00055. Figure 2 shows that model system (4) is globally asymptotically stable whenever condition p < μ/cβ is satisfied. That is, there are no periodic solutions

Figure 2

Illustration of nonexistence of periodic solutions when condition p < μ/cβ holds true. The corresponding R0 with the given parameter values is R0 = 1.5118 > 1.

Now, we choose another set of parameter values that satisfy the condition p > μ/cβ, namely, μ = 0.017, μ = 0.1, k = 0.001, γ = 2, q = 0.05, p = 0.135, c = 60, β = 0.519, and k0 = 0.019. The initial conditions used are the same as in Case (i) above. Figure 3 depicts that model system (4) has sustained periodic solutions if condition p > μ/cβ

Figure 3

Illustration of existence of periodic solutions when p > μ/cβ. The corresponding R0 with the given parameter values is R0 = 1.5118 > 1.

7. 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.