Lyapunov Functions for Tuberculosis Models with Fast and Slow Progression

The spread of tuberculosis is studied through two models which include fast and slow progression to the infected class. For each model, Lyapunov functions are used to show that when the basic reproduction number is less than or equal to one, the disease-free equilibrium is globally asymptotically stable, and when it is greater than one there is an endemic equilibrium which is globally asymptotically stable.


1.
Introduction. According the World Health Organization, one third of the world's population is infected with tuberculosis (TB), leading to between two and three million deaths each year. For most individuals infected with TB, the immune system is able to control the causative agent, Mycobacterium tuberculosis, but not eliminate it. These individuals are not infectious and suffer no symptoms, although they usually test positive on a skin test. However, it is possible that after a latent period of years or decades, these individuals may become symptomatic and infectious. There is also a smaller fraction of individuals for whom the progression to active TB is much faster. This fast progression is particularly common for individuals with a compromised immune system.
In [1], two models that include fast and slow progression of TB are studied, but the global stability of these models has remained an open problem. Because the models are three-and five-dimensional, a rigorous demonstration of global stability has not been expected. More recently a survey of TB modelling [2] presented a list of "challenging questions" which included resolving the global stability of the three-dimensional model in [1].
In this paper, the global dynamics of both the three-dimensional model and the five-dimensional model of fast and slow progression in [1] are resolved through the use of Lyapunov functions.
The Lyapunov functions used in this paper to demonstrate the stability of the endemic equilibria are of the same form as those used recently in [4] and [5] to determine the global dynamics of SEIR, SEIS, and SIR models. It should be noted, though, that those works do not include both fast and slow progression. This distinction is important since the progression term is the non-linearity which is the source of any non-trivial dynamics. Before [4] and [5], a Lyapunov function of the same form was used in [3] and [7], for example, to study Lotka-Volterra systems for predator-prey interactions.

2.
A Simple Model with Fast and Slow Progression.

Model Formulation.
A population is divided into three classes based on epidemiological status; individuals are classified as either susceptible, exposed, or infectious. The sizes of these groups are represented by S, E and I, respectively. All recruitment is into the susceptible class, and occurs at a constant rate Λ. Death for reasons that are not related to the disease is proportional to the population size, with rate constant µ. The additional death rate due to disease affects only class I and has rate constant d. The time before exposed individuals become infectious is assumed to satisfy an exponential distribution, with mean waiting time 1 k . Thus individuals leave class E for class I at rate kE. In each time unit, a susceptible individual has on average βI contacts that would suffice to transmit the disease. Thus the rate at which susceptibles are infected is βSI. A fraction p of the newly infected individuals are assumed to undergo fast progression directly to the infectious class, while the remainder are latently infected and enter the exposed class. The transfer diagram is given below.
This yields the following set of differential equations: The behaviour of this system is studied in [1], but the global stability is not treated there; we resolve the global dynamics here. Using the next generation matrix method [8], the basic reproductive number for system (2.1) is Remark 1. This is the same expression as given in [1] but differs from that which appears in [2], which is apparently missing a factor β µ+d from the final term in their expression (6).
Under the flow described by (2.1), the region D = {(S, E, I) ∈ R 3 ≥0 : S + E + I ≤ Λ/µ} is positively invariant. Furthermore, each solution in R 3 ≥0 approaches D and so we may restrict our analysis to this region. For R 0 ≤ 1 the only equilibrium is the disease-free equilibrium P 0 = (Λ/µ, 0, 0) ∈ D. For R 0 > 1, P 0 is present as is an additional equilibrium P * = (S * , E * , I * ) ∈ D with S * , E * , I * > 0. Any solution which has an initial condition in R 3 ≥0 for which E +I is positive, immediately moves to the interior of the positive orthant. On the other hand, if the initial condition in R 3 ≥0 satisfies E = I = 0, then the solution limits to P 0 . 2.2. Global Stability of P 0 for R 0 ≤ 1. Consider the Lyapunov function Differentiating with respect to time yields On D, we have S ≤ Λ/µ, and so with equality only at P 0 . For R 0 ≤ 1, this gives with equality only if I = 0. By LaSalle's extension to Lyapunov's method [6], the limit set of each solution is contained in the largest invariant set for which I = 0, which is the singleton {P 0 }. and Note that equations (2.3) and (2.5) are only valid for p = 1. Nonetheless, the Lyapunov function which follows (with A and B given by equation (2.9)) does work for p = 1, but it is necessary to do the calculation separately for that case. As that case is much simpler than the case where p = 1, we omit it. We now follow the approach used in [4] and [5], and consider the Lyapunov function where A and B are constants that will be given below. Differentiating U gives

Using (2.2), (2.3), and (2.5) to rewrite this, we get
(2.8) We now show that it is possible to choose A = A(p) and B = B(p) such that f is non-positive for all x, y, z ∈ R >0 . In order to make the coefficients of z, y, and xz equal zero we choose Substituting (2.9) into (2.8) and rearranging gives Differentiating equation (2.10) with respect to p yields If x, y, and z are fixed, then ∂f ∂p has constant sign for p ∈ [0, 1]. Thus, f is maximized at p = 0 or at p = 1.
Suppose p = 1. Then filling into (2.10) gives which is less than or equal to zero by the arithmetic mean-geometric mean inequality, with equality if and only if x = 1. Similarly, if p = 0, then (2.10) becomes which is also less than or equal to zero by the arithmetic mean-geometric mean inequality, with equality if and only if x = 1 and y = z. Thus (2.7) implies that U is less than or equal to zero with equality only if S = S * . LaSalle's Extension implies solutions which intersect the interior of R 3 >0 : S = S * } is the set consisting of the endemic equilibrium P * . Thus, all solutions to equation (2.1) which intersect the interior of R 3 >0 limit to P * . In fact, P * is globally asymptotically stable for all non-negative initial conditions for which E + I > 0.

A Detailed Model with Fast and Slow Progression.
3.1. Model Formulation. We now consider a more detailed model of fast and slow progression for tuberculosis, which also appears in [1]. In addition to the structure that is assumed for system (2.1), we include two additional classes: recovered and actively infected but not infectious. The sizes of these two groups are R and T , respectively. Of the newly infected individuals that undergo fast progression to TB, a fraction f are infectious and the remaining fraction 1 − f are non-infectious. A fraction q of the individuals that undergo slow progression after a period of latency are infectious, and the remainder are non-infectious.
Individuals that are actively infected (infectious or non-infectious) can spontaneously recover from the disease with rate constant c, entering the recovered class. Recovered individuals still have the bacterium in their body and can undergo a reactivation of the disease with rate constant 2ω, with half of the reactivated cases being infectious and half being non-infectious.
The behaviour of this model is studied in [1], but the global stability is not treated there; we resolve the global dynamics here.
The basic reproduction number for system (3.2), is given in [1]. R 0 can also be found by using the next generation method of [8]. approaches B and so we may restrict our analysis to this region. The disease-free equilibrium Q 0 = (Λ/µ, 0, 0, 0, 0) is present for all parameter values. For R 0 > 1, there is an additional equilibrium Q * = (S * , E * , I * , T * , R * ) ∈ B, which is called the endemic equilibrium and which satisfies S * , E * , I * , T * , R * > 0. Any solution which has an initial condition in R 5 ≥0 for which E + I + T + R is positive, immediately moves to the interior of the positive orthant. On the other hand, if the initial condition in R 5 ≥0 satisfies E = I = T = R = 0, then the solution limits to Q 0 . 3.2. Global Stability of Q 0 for R 0 ≤ 1. Consider the Lyapunov function is positive. Differentiating with respect to time yields On B, we have S ≤ Λ/µ and so with equality only at Q 0 . For R 0 ≤ 1, this gives where A, B, C, D > 0 will be given below. The calculation to show that V ≤ 0 on R 5 >0 is long and unpleasant, but runs parallel to the corresponding calculation from the previous section. We provide an outline of the calculation here, omitting details.
In the expansion of V , the right-hand side of each of the five equations in (3.2) appears. In the right-hand side of the third and fourth, we replace (c + µ + d) with pf βS * I * + qkE * + ωR * /I * and p(1 − f )βS * I * + (1 − q)kE * + ωR * /T * , respectively, and in the right-hand side of the fifth, c is replaced with (µ+2ω)R * /(I * +T * ). These substitutions come from the fact that I = T = R = 0 at Q * . Also, since S = 0 at Q * , we can eliminate Λ from the expansion of V by using the substitution Λ = βS * I * + µS * . Next, using E = 0 at Q * we make the substitution β = (k + µ)E * / (1 − p)S * I * .
We can now write the time derivative of V as and α = I * I * +T * . From equation (3.2) at Q * , it can be shown that We now work towards showing that V is less than or equal to zero, with equality only at S = S * . In order to show that this is true, we demonstrate that the same statement holds for F . We introduce new variables x = S S * , y = E E * , z = I I * , w = T T * and u = R R * to eliminate S, E, I, T and R. The resulting expression for F still contains the fraction R * E * . By considering (3.2) at Q * , it can be shown that R * Substituting into (3.5) gives (3.8) Next, we choose A, B, C, and D so that none of the variable terms (such as x and xz w ) are positive. To do this, it is useful to first group together many of the terms in F that involve the same variable term, as well as grouping all of the constant terms together. Rearranging (3.8) yields (3.9) The only variable terms that appear in (3.9) with positive coefficients are z, xz, y, u, and w. If the total of the coefficients of any of these variable terms is positive, then it is conceivable that F could be positive. Thus, we choose A, B, C, and D so as to make the coefficients of z, xz, y, u, and w non-positive. This yields a system of five inequalities in four parameters, which must be simultaneously satisfied. It By considering all quantities except p, q, and f to be fixed, it is clear that F is maximized at either f = 0 or f = 1, since ∂ F ∂f = δ 1 p has constant sign. Also, ∂ F ∂q = δ 2 p + δ 4 = δ 2 (p − 1) has constant sign; thus F is maximized at either q = 0 or q = 1. Now, if f and q are also considered to be fixed, it is clear that F is monotone in p and so is maximized at either p = 0 or p = 1. Thus, F (p, q, f ) is maximized at one of the corners of the unit cube [0, 1] 3 . We will now demonstrate that F is non-positive at each of the corners of the unit cube and is therefore non-positive on the entire unit cube, implying that the same is true of F and of V . Suppose p = q = f = 0. Then F (p, q, f ) = δ 5 . Considering the expression for δ 5 given in (3.14) and applying the arithmetic mean-geometric mean inequality separately to the two terms in parentheses, it is clear that Since S must remain constant at S * , S is zero. This implies that I = I * , making I I * equal to one. Thus, T = T * and R = R * . Now, we must have I equal to zero, implying that E is equal to E * . Thus, the only invariant set contained in C is the singleton {Q * }. This shows that each solution which intersects R 5 >0 limits to the endemic equilibrium Q * , giving the following result. If R 0 ≤ 1, then the equilibrium Q 0 is globally asymptotically stable on R 5 ≥0 . If R 0 > 1, then there is a unique endemic equilibrium which is globally asymptotically stable on R 5 ≥0 \ A, with solutions in A limiting to Q 0 .
4. Discussion. The global stability for systems (2.1) and (3.2) has now been resolved. If R 0 ≤ 1, then each solution limits to the disease-free equilibrium; the disease dies out of the population. If R 0 > 1, then there is a unique endemic equilibrium which is globally asymptotically stable among all states for which the disease is present; if disease is present in the population, then it will persist. For each system, there is a step in the proof of the global stability of the endemic equilibrium where there is a linear system of n inequalities in n−1 unknowns (where n is 3 or 5). Assuming that the inequalities are linearly independent (as is the case here), such a situation may have either no solutions, a unique solution, or an infinite number of solutions. The fact that for each of the systems studied here there is a unique solution implies there is some important underlying structure. A more detailed analysis of this phenomenon may yield much more general results.
Furthermore, the coefficients that are used for the Lyapunov function for R 0 ≤ 1 are a scalar multiple of those used for R 0 > 1, but with the coefficient for S set to zero. Again, this points to deeper structure.