Epidemic models with discrete state structures.

The state of an infectious disease can represent the degree of infectivity of infected individuals, or susceptibility of susceptible individuals, or immunity of recovered individuals, or a combination of these measures. When the disease progression is long such as for HIV, individuals often experience switches among different states. We derive an epidemic model in which infected individuals have a discrete set of states of infectivity and can switch among different states. The model also incorporates a general incidence form in which new infections are distributed among different disease states. We discuss the importance of the transmission-transfer network for infectious diseases. Under the assumption that the transmission-transfer network is strongly connected, we establish that the basic reproduction number R 0 is a sharp threshold parameter: if R 0 ≤ 1 , the disease-free equilibrium is globally asymptotically stable and the disease always dies out; if R 0 > 1 , the disease-free equilibrium is unstable, the system is uniformly persistent and initial outbreaks lead to persistent disease infection. For a restricted class of incidence functions, we prove that there is a unique endemic equilibrium and it is globally asymptotically stable when R 0 > 1 . Furthermore, we discuss the impact of different state structures on R 0 , on the distribution of the disease states at the unique endemic equilibrium, and on disease control and preventions. Implications to the COVID-19 pandemic are also discussed.

Introduction

Many infectious diseases have a long course of progression in infected hosts, often through several distinct stages. These include viral infections such as HIV and viral Hepatitis, bacterial diseases such as Tuberculosis and Cholera, parasitic diseases such as Malaria and Schistosomiasis. Infected individual hosts many experience changes in levels of pathogen load, which impacts the individual’s infectivity, and levels of antibodies and immune cells, which impact host immunity, either due to the natural disease progression or through medical interventions. Individual variations of these important factors are key sources of heterogeneity in the host population, and can be essential for understanding the dynamical complexity of the disease transmission and spread.

For an infectious disease, we consider the state of an individual host that represents all or a subset of important epidemiological, immunological, and ecological factors for disease progression and transmission. In the case of HIV-1 infection, an individual host’s state of infection can be simply represented by a vector-valued variable , where in the level of HIV-1 viral load and the level of CD count, common measures used to gauge the HIV-1 progression of a patient. The values of can vary over time as the infection progresses and when the antiretroviral treatments (ART) are given, and they typically exhibit a high degree of individual variations.

The concept of state of infection of an individual host was used in the classical paper of Kermack and McKendric [1], and has since been explored in mathematical models under different terms. In staged progression models of HIV, an infected individual can progress forward through different stages or backward due to treatment [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12]. The stages of disease progression is a good example of the state of an infected individual. The concept of state of infection is, however, more general than that of stages and can unify many different approaches to model disease heterogeneity. For instance, the concept of state or substate can be applied to susceptible hosts to represent different level of host susceptibility to infection or reinfection after recovery. Such a consideration will lead to the investigation of the impact of differential susceptibility on disease transmission (see [4]). The concept of state or substate can also be considered for vaccinated or immune hosts and used to measure individual variations as well as community distributions of level of antibodies and immunity. The concept of state can also be extended to individual characteristics including spatial location, chronological age, development stage, size etc [13]. And the distributions of individual characteristics directly influence the dynamics of population which could be generated as the state space of measures [14], [15], [16]. Interpretation of the traditional epidemic models in a state-space framework was used to track the evolution of Flu in the United States [17] and predict intervention effect for COVID-19 in Japan [18].

The concept of state has also been widely employed to diffusion models in the social and life sciences. For example, the Bass model for diffusion of innovation [19] and its generalizations and variations [20], [21], models for the diffusion of information on social media [22], models for the spread of opinions [23], and models for the spread of innovations and transient fads [24]. Further information on this type of models can be found in more recent references (e.g. [25], [26], [27]).

Recent HIV research on nonprogressors and elite controllers [28], [29], HIV-1 patients who can keep the HIV-1 progression under check without anti-retroviral treatment, suggests that to fully understand the progression of HIV-1 infection, we need to include in the disease state indicators of the immune system such as counts of HIV-1 specific cytotoxic T lymphocytes (CTLs), levels of anti-HIV-1 antibodies, as well as measures of anti-HIV-1 cytokines and other biomarkers in modelling studies. From a transmission viewpoint, the disease state should also include an individual’s sexual preferences and a scoring of risky behaviours. This approach can stimulate a more comprehensive approach to modelling and better inform the prevention and control measures of the HIV epidemic (see e.g. [30], [31]).

The transmission dynamics of the SARS-CoV-2 virus, which is responsible for the current COVID-19 pandemic, provided another good example of the importance of the state structure for epidemics. The significant number of asymptomatic infections is believed to be a key reason why the COVID-19 epidemics have proven to be more difficult to control than the SARS epidemics in 2003 [32], [33]. The substate of asymptomatic infections can be studied in state-structured models. The current vaccines for COVID-19 are known to be highly effective to protect against serious symptoms and diseases, but their efficacy against transmission is still largely unknown [34]. While the vaccines can prevent serious diseases and deaths among seniors and people with underlying health, the vaccination of younger population may produce more asymptomatic disease carriers, whose existence was already known for COVID-19 even before the vaccines were available [33], if the efficacy of the vaccines against transmission is not sufficiently high. The impact of subgroup of asymptomatic carriers has been studied in the literature (e.g. [35]) and more recently for COVID-19 in the class of A-SIR models [36], [37]. Other subgroups that capture the many characteristics of the COVID-19 were also considered in recent COVID-19 modelling, including undetected infections, infected but not infectious (latent), quarantined or isolated individuals, and hospitalized individuals (see e.g. [36], [38], [39] and references therein). Subgroups of infected individuals are special cases of the infection state and can be investigated under the framework of state-structures.

In our paper, we restrict our consideration to the state of infection among infected individuals and allow infected individuals switch their state of infection among finitely many different substates, hence creating a finite discrete state-structure. We formulate and investigate a class of epidemic models with a discrete state structure among infected individuals. Switch among different sub-states is given by a state-transfer matrix, which helps to capture the heterogeneity among the substates. We note that sub-state structures given by different state-transfer matrices may alter the probability distribution of the duration of infectiousness, and we refer the reader to the approach of the linear chain trick and its generalizations (see e.g. [40]) for this linkage. For modelling dynamical processes on networks, we refer the reader to the approximation frameworks of binary-state dynamical processes on networks [41] and its generalization [42].

In addition to carrying out a standard stability analysis of a relatively comprehensive model, the consideration of state has inspired us to investigate a new set of questions and problems that are of both biological and mathematical interests: (1) the impacts of state structures on the basic reproduction number ; (2) the influence of different disease interventions on the population distribution among different disease states; and (3) global stability analysis is of a considerable mathematical challenge for the state-structure models, and the graph-theoretic approach to the construction of Lyapunov functions [43], [44] was further improved and adapted for the proof of global stability.

Results in our paper include as special cases earlier results on multi-staged epidemic models (see e.g. [9] and reference therein). A continuous approach to modelling the state structure was considered in [45], and the resulting model is a system of differential–integral equations with nonlocal effects. In our model (2.1) with a discrete state structure in Section 2, mathematical results on global stability are obtained under more relaxed conditions than those for infinite dimensional systems as in [45].

The paper is organized as follows. We formulate our model in the next section. In Section 3, we derive the basic reproduction number and give its biological interpretations. In Section 4, we carry out mathematical analysis of the global dynamics of the model for the case : the disease free equilibrium is globally asymptotically stable when ; and a unique endemic equilibrium exists and is globally asymptotically stable when . In Section 5, we carry out numerical investigations of the impacts of the state structure on disease dynamics and outcomes.

Model formulation

We assume that the state of an infected individual is a finite set with a scalar index , . For , let denote the number of the infected individuals who are in the th state at time . Let denote the number of susceptible individuals at time , and denote the number of removed individuals at time . Here, an individual is considered removed if the individual no longer participates in the transmission process. This can happen if an individual is in the terminal illness, such as having AIDS in the case of HIV infection, or the individual is permanently immune from reinfection, as in the case of Measles, or individuals who are tested positive for the infection and isolated or hospitalized as in an COVID-19 epidemic. We use to represent the total number of individuals who are active in the infection process.

In the absence of the disease, we assume that the dynamics of population are described by a nonlinear differential equation: , where is a growth function with a carrying capacity. Common forms of are (see [6], [46], [47]) and (see [9], [48]). We also assume that the disease transmission is horizontal and has a general nonlinear incidence function , where denotes the density dependence. Out of the total new infection , we assume that a fraction will enter the th state, and . Let denote the removal rate from , which may include natural death, disease-caused death, and out-migration. A common form of is the exponential removal . To describe the switches among different states, for given , let denote the transfer rate from the th state to the th state, and by convention, . For , the transfer rate from the th state to the removed compartment is given by . The removal rate from the compartment is given by . The preceding assumptions are illustrated in the transfer diagram in Fig. 1.

Fig. 1

The transfer diagram of model (2.1). Here the incidence terms are given by .

Based on these assumptions and using the transfer diagram, we derive an epidemic model with a discrete state structure described by the following system of nonlinear differential equations: Since the removed population is no longer part of the transmission process and the variable does not appear in the equations of and , we consider the following equation for the removed population separately:

The general functions in model (2.1), and are assumed to be sufficiently smooth such that existence and uniqueness of solutions are satisfied. We further make the following biologically motivated assumptions:

There exists such that and , ;

is a nonincreasing positive function;

For , for all ; and if then iff or ;

For , for , , and iff ;

For , for all and ;

For , , and there exists such that for all ;

There exist constants , such that ;

There exist constants such that ;

There exist constants , , .

It can be verified by examining direction of the vector fields on the boundary of that solutions to system (2.1) with nonnegative initial conditions remain nonnegative for and that the model is well defined. From the first equation of system (2.1) and assumptions , , we know that . Assumption implies that , a carrying capacity of the disease-free population. Adding all equations of (2.1) yields that By assumption , we have Let and we assume that . Then where , which implies Therefore, omega limit sets of system (2.1) are contained in the following bounded region: It can be verified that is positively invariant with respect to (2.1).

The function in model (2.1) represents density dependence in the disease incidence. The monotonicity assumption on in is satisfied by the class , used in the literature. It is however a simplifying assumption. When is non-monotone, it was shown in [49] using a simple SIR model that complex dynamics such as backward bifurcations and Hopf bifurcations can occur. Assumption on the incidence function allows many possibilities for disease transmission at different state. For instance, transmission coefficients can be at states where the infected individual is latent or quarantined, and transmission coefficients can be substantially lower at a state where infected individual are either asymptomatic or protected by a vaccine.

Our model (2.1) include many previously studied compartmental models that have infectious, symptomatically infectious, asymptomatically infected, latent, vaccinated, quarantined, and isolated compartments. When and for , our model can be interpreted as a multi-staged model studied in [9], in which all new infections enter the first stage. In staged progression models, the stages are understood to have a linear sequential structure. The transfer terms , when , denote the rate of disease progression (forward); and when , denote the rate of amelioration (backward). For state-structured models, a sequential order among states is no longer essential, and different states can be regarded as parallel structures. A key difference between state-structured models and staged-progression models is that a newly infected individual can manifest in any disease state while often enters into the first disease stage.

Graph theory preliminaries

Graph theory provides the mathematical language and tools for describing the complex structures of state transfers and transmission.

A directed graph consists of a set of vertices and a set of directed edges that connect ordered pairs of vertices. A nonnegative weight can be assigned to the directed edge from vertex to , then we obtain a weighted directed graph associated with weight matrix , . A directed edge exists if and only if . On the other hand, given a nonnegative matrix , we obtain a weighted directed graph with vertices. For simplicity of notation, we will suppress the vertex and edge sets from the notation of a directed graph, only focus on its weight matrix . Directed graph is said to be strongly connected if each pair of vertices are joint by a directed path.

A nonnegative matrix is reducible if there exists permutation matrix such that where and are square matrices. Matrix is irreducible if and only if directed graph is strongly connected [50].

The algebraic Laplacian matrix of is defined as: For each , let denote the co-factor of the th diagonal entry of . Kirchhoff’s Matrix Tree Theorem gives a graph-theoretic description of co-factor [51]. Let be the set of all spanning sub-trees of that are rooted at vertex , and be the weight of , the product of weights on all edges of .

Matrix Tree Theorem [52]

Assume that . Then Furthermore, if is strongly connected, then , for .

Let be a set of functions, and the set of all spanning unicyclic subgraphs of , with a unique directed cycle . Let be the weight of and the set of arcs in . The following tree cycle identity was first derived and proved in [44], [53].

Tree Cycle Identity [44], [53]

Assume that and is given in Proposition 3.1 . Then the following identity holds:

The identity in the next corollary was first given in [44] and follows from the tree cycle identity.

Let be any family of functions, and are given in Proposition 3.1 . Then the following identity holds:

Equilibria and the basic reproduction number

Equilibria

Set By assumptions -, if and only if , for . System (2.1) can be rewritten as the following:

It follows from assumptions - that the disease-free equilibrium always exists. Let denote a possible endemic equilibrium of (4.2). Then, the coordinates are positive solutions to the following system of equilibrium equations:

We note that endemic equilibria of model (4.2) can be of two different forms: either a positive endemic equilibrium for which for all , or a mixed endemic equilibrium for which for some and for the remaining . A positive endemic equilibrium belongs to the interior of and a mixed endemic equilibrium belongs to the boundary of . In the case of a positive equilibrium, the disease is endemic among all states, while in the case of a mixed endemic equilibrium, the disease is only endemic among a subset of states, and the remaining states are disease free. It was known in the literature that complex disease models can have mixed endemic equilibria (see e.g. [54], [55]). It was pointed out in [55] that mixed endemic equilibria exist when the transmission network of the disease is not strongly connected. Let Using assumptions and , we define where and let and .

In weights defined in (4.5), the term comes from the transmission, and the term defines the state-transfer matrix . The directed graph is called the transmission–transfer network of our model (4.2), and the directed graph is the state-transfer network. From (4.5) we see that and have the same vertex set, but not necessarily the same edge set. Furthermore, if is strongly connected, so is , but the inverse may not be true.

We note that the weight matrix of the transmission–transfer network is obtained using the limiting values of weight functions in (4.4) when point approaches . In later sections, we will also consider weight matrices with measured at endemic equilibrium in . By assumptions -, it can be verified that directed graph is strongly connected if and only if the transmission–transfer network is strongly connected.

In Fig. 2, we show a more intuitive visualization of the transmission–transfer network: the vertex is labelled as and an additional vertex is added so that the resulting directed graph with vertices can be easily identified from the transfer diagram of model (2.1) in Fig. 1. We call this new directed graph by adding vertex an augmented directed graph of . In the weight from to in (4.5), represents the populations transfer from state to state . The term is less intuitive as weight from state to . In the augmented directed graph with vertices in Fig. 2, we can see that this term is created indirectly through infection of by as indicated by a dashed edge from to , together with the corresponding incidence in as indicated by a solid edge from to . This is also explained in Fig. 3-(b) and Fig. 3-(d). It can be verified that the weighted directed graph is strongly connected if and only if its augmented directed graph in Fig. 2 is strongly connected.

Fig. 2

The augmented directed graph with a vertex added as a realization of the transmission–transfer network for model (4.2). The vertex is added to demonstrate that the additional weight on the edge from to is created indirectly through infection of by indicated by a dashed edge from to , and the corresponding incidence in indicated by a solid edge from to .

Fig. 3

The transfer diagram (a) of an SEIR model, its augmented directed graph (b), state-transfer network (c) and transmission–transfer network (d).

In Fig. 3, using a simple SEIR model, we illustrate the relations between the model transfer diagram, the augmented directed graph, the state-transfer diagram and transmission–transfer network. We note that the transmission–transfer network in Fig. 3-(d) is strongly connected, while the state-transfer network between the latent state and infectious state in Fig. 3-(c) only contains a single arrow from to , and is not strongly connected.

Assume that the transmission–transfer network of (4.2) is strongly connected. Then the disease-free equilibrium is the only equilibrium of (4.2) on the boundary of .

The proof is similar to that of Theorem 2.3 in [55], in which only a transmission network is considered. It suffices to exclude the existence of mixed endemic equilibria on the boundary of .

Let denote a nonnegative equilibrium of (4.2), where . Using assumption and the first equilibrium equation of (4.3), we have .

If an arc from to exists and , we show that . Assume otherwise that , then assumptions - imply that , and assumption and imply that . Using the equation for we obtain This contradicts the fact that is an equilibrium, and shows .

If an oriented path from to exists and , using statement above repeatedly, we can show that .

Since the transmission–transfer network is strongly connected, for given two states , there exists and oriented path from to and from to . From the preceding arguments, we know that implies . Therefore, for some implies for all , namely, for all , are either all zero or all positive, and thus the only equilibrium on the boundary is the disease-free equilibrium .  □

If is not strongly connected, then there may exist a state that can reach other states but is not reachable from any other state. Then, when , there exist state for which , and mixed endemic equilibria exist on the boundary. We give the following example to demonstrate this.

Consider the epidemic model (4.6) illustrated by the transfer diagram in Fig. 4(a). Its transmission–transfer network is not strongly connected, and it has two boundary equilibria.

Fig. 4

The transfer diagram (a), augmented direct graph (b), and transmission–transfer network (c) of model (4.6) are shown. The node is an unsustainable external source of infection that feeds into the system but does not receive feedback. The term can be interpreted as importations of infections from other regions that have implemented effective control measures. The augmented direct graph in (b) is not strongly connected because other nodes cannot reach , and correspondingly, the transmission–transfer network in (c) is not strongly connected. The model (4.6) has two boundary equilibria.

To see the claims in Example 4.1, we note that the state is not reachable from any other states either by state transfer or transmission, and the transmission–transfer network is not strongly connected (Fig. 4). The node in the transfer diagram (Fig. 4(a)) can be thought of as an external source of infections that feeds into the system , and the external infection source is not sustainable due to effective control measures. This can be the case when other regions have the epidemics under control and the importations into the region from outside will dwindle to zero. From the model equation (4.6) we can see that two boundary equilibria can exist: the disease-free equilibrium , a mixed endemic equilibrium , .

In this paper, we only consider the case when the transmission-transfer network is strongly connected. In the case when is not strongly connected, the dynamics of model (4.2) can be analysed following the approach in [55] using the condensed graph.

The basic reproduction number

The basic reproduction number is the expected number of secondary infections produced by a typical infectious individual during its entire infectious period in an entirely susceptible population. The mathematical definition of the basic reproduction number is by the spectral radius of the next generation matrix, see [56], [57].

Let be defined at the beginning of Section 4.1 and Following [57], the basic reproduction number for model (4.2) is the spectral radius (dominant eigenvalue) of next generation matrix , where denotes the spectral radius of a matrix , and Using (4.7), and , it can be verified that is a nonnegative matrix, and is a nonsingular -matrix with a nonnegative inverse .

Assumptions - are sufficiently weak to allow the nonlinear incidence function and transfer functions , , and to be not differentiable at the disease-free equilibrium. For instance, a common incidence function satisfies these assumptions with , but it is not differentiable at . When these functions are differentiable at , constants are given by the corresponding derivatives at .

In the most common case when , , , and , assumptions - are automatically satisfied with It can be verified that the definition of basic reproduction number in (4.8) includes as special cases the corresponding definitions for simpler models in the literature (see e.g. [5], [6], [7], [8], [9], [46], [47], [58], [59]).

Biological interpretations of

The basic reproduction number in (4.12) can be interpreted as the sum of the expected number of secondary infections contributed by an infected individual at each infection state.

Suppose that assumptions - hold. Let Then , , and . Furthermore, the following results hold:

Nonnegative matrix has a unique positive eigenvalue .

The basic reproduction number satisfies

The matrix can be written as By (4.11), (4.13), we have . Then matrix has rank , and thus its spectral radius is the only positive eigenvalue and is equal to its trace , leading to (4.12).  □

As pointed out in [57], the entry of matrix represents the average time of an infected individual originally at the th state spends in the th state, the entry of matrix is the incidence rate at which an infected individual in the th state cross infect in the th state. Overall, the entry of matrix is the expected number of new infections into the th state made by an infected individual of the th state during its infectious period, accounting for the time this individual stays in all states.

As defined in (4.11), is the sum of the th column of matrix , and represents the expected number of new infections produced by a newly infected individual at the th state during its entire infectious period.

To interpret in (4.12) biologically, consider a newly infected individual, with probability the individual is at the state , , and will produce expected secondary infections during its infectious period. Therefore, the total expected number of secondary infections during its infectious period is

Consider the special case when the proportions of incidence among states is such that , for , namely, all new infections enter the first state. In such a case, disease states can simply represent stages of disease progression. The resulting model agrees with the staged progression model considered in [9]. By (4.12), , which agrees with the basic reproduction number given in [9]. Similarly, if , then

In the next result, we show that the basic reproduction number can also be interpreted using the average time spent by an infected individual at each disease state.

Assume that assumptions - hold. Let Then the following results hold:

and .

is the average infectious period of an infected individual.

The basic reproduction number satisfies

By (4.11), (4.12), Then, (4.15) follows from (4.14).

Since the entry of matrix represents the average time of an infected individual originally at the th state spends in the th state, we know that in (4.14) is the average time an infected individual spends in the th state. For each , denotes the force of infection at the beginning of an epidemic with respect to an infected individual at the th state, and represents the expected number of secondary infections produced by an infectious individual who is at state , and the basic reproduction number is given by (4.15). Moreover, is the average infectious period of an infected individual.  □

Upper and lower bounds on

Let be defined in (4.10). Define Biologically, represents the basic reproduction number of th state when the other states are absent. For instance, when , , , and , then we arrive at the standard expression for the basic reproduction number of a single group SIR model: The following proposition gives a lower and upper bound on the basic reproduction number .

Suppose that assumptions - hold. Then the basic reproduction number defined by (4.12) satisfies the following inequality:

Let for some . By (4.7), (4.9) we obtain that . Since matrix is nonnegative and using (4.11), we have By (4.12), we have . Similarly, we can show , completing the proof.  □

Global dynamics and threshold results

In this section we carry out analysis of the global dynamics of system (4.2), and establish a sharp threshold result: if , then the disease-free equilibrium is globally asymptotically stable in , and the disease dies out; if , and if , then the positive endemic equilibrium is unique and globally asymptotically stable in the interior of . In this case, the disease persists in the population and has a positive stationary distribution among all disease states. The proof relies on construction of global Lyapunov functions. We have adapted the graph-theoretic approach developed in [43], [44] to system (4.2).

Global dynamics when

We first describe the Lyapunov function for the case .

Let is defined by (4.11) . If holds in , and , then is a Lyapunov function for system (4.2) in , namely, along all solutions of (4.2) in .

Differentiating function with respect to , we have Using (4.11), (4.12) we obtain completing the proof.  □

We make the following assumptions:

For holds for all ; if , then ; for , .

For all , ; for all , .

It can be verified that these assumptions are satisfied by common incidence forms and linear transfer functions. Since biologically the incidence function should be nondecreasing in , assumption is satisfied if either or is strictly increasing in . Under assumptions and , we can show that condition (5.1) is satisfied and the function defined in (5.2) is a Lyapunov function as stated in the next result. The proof is technical and is given in the Appendix.

Assume that and hold. Then function is a Lyapunov function for system (4.2) in .

Suppose that assumptions - are satisfied. The following statements hold:

If and , hold, then the disease-free equilibrium is globally asymptotically stable in .

If then the disease-free equilibrium is unstable. Furthermore, if the transmission–transfer network of model (4.2) is strongly connected, then system (4.2) is uniformly persistent with respect to the interior of .

By Corollary 5.1, function is a Lyapunov function for (4.2), and for all . Let be the largest invariant subset of . Then . Let be a solution in . We consider two cases. Case (1): . In this case, if and only if . By assumption , this means that , hence , and . Using assumption , we have Therefore, for all . Case (2): . In this case, implies , and thus by assumption . Moreover, along any solution in , holds. Using the same procedure as in Case (1), we obtain that for all . In both cases, we have shown that . By LaSalle’s Invariance Principle [60], is globally asymptotically stable in when . This establishes part (a).

If , then, by continuity, for sufficiently close to except when . Solutions starting sufficiently close to leave a small neighbourhood of , except for those on the invariant -axis. Since the transmission–transfer network of model (4.2) is strongly connected, is the only equilibrium on the boundary of by Proposition 4.1. By a uniform persistence result, Theorem 4.3 in [61], and using a similar proof of Proposition 3.3 in [62], we can show that instability of is equivalent to the uniform persistence of system (4.2) with respect to the interior of . This completes the proof.  □

Uniform persistence of system (4.2) in and uniform ultimate boundedness of solutions in imply that there exists an equilibrium in the interior of , see [63] (Theorem D.3) or [64] (Theorem 2.8.6).

Suppose that assumptions - hold. Assume that and that the transmission–transfer network of system (4.2) is strongly connected. Then system (4.2) has an endemic equilibrium in the interior of .

The assumption that the transmission–transfer network is strongly connected is necessary for ruling out boundary equilibria other than , and in turn, crucial for the proof of uniform persistence and existence of positive endemic equilibrium when . This condition was neglected in many earlier work on complex epidemic models.

We assume that in this subsection and establish that, when , system (4.2) has a unique endemic equilibrium and it is globally asymptotically stable in the interior of .

System (4.2) is a large-scale system of nonlinear equations. Proof of global stability for systems of this type is typically done using the method of Lyapunov functions. The graph-theoretic approach for constructing global Lyapunov functions developed in [43], [44] has been successfully applied to many complex systems. An application of the approach to multi-staged models was shown in [9]. Comparing to the multi-stage model in [9], system (4.2) has disease incidence terms in all equations and are more interconnected. The graph-theoretic approach needs to be adapted to (4.2) to resolve its global-stability problem.

Let , , be an endemic equilibrium. Using the weight function in (4.5), we define a nonnegative matrix , with The corresponding weighted directed graph has the same vertex and edge sets as the transmission–transfer network in Section 4.1, but the weights are measured at the endemic equilibrium . Using and , we know that is strongly connected if and only if is strongly connected.

Motivated by method in [9], we construct a candidate Lyapunov function of form where is given in Proposition 3.1 and . In the special case when , and for , the Lyapunov function in (5.5) reduces to the Lyapunov function for the staged progression model in [9]. The proof of our next result is technical and provided in the Appendix.

Assume that . Suppose that there exists Lipschitz continuous function such that the following assumptions hold:

For ,

For ,

For ,

Then the function defined in (5.5) is a Lyapunov function for system (4.2) with respect to the interior of .

We remark that assumption is automatically satisfied if transfer functions and are linear. If the incidence function is of separable form , we can choose . Then assumption is satisfied if is strictly increasing in , as in the case for or , , . If the growth function takes the form , then assumption is satisfied if is strictly increasing in .

For our next result, we make the following assumption, which is satisfied by common incidence functions and linear transfer functions.

For each , one of the functions is strictly monotone in .

Suppose that assumptions - , - hold, , and the transmission–transfer network of model (4.2) is strongly connected. If , then the positive endemic equilibrium is unique and globally asymptotically stable in the interior of .

Choose the Lyapunov function defined by (5.5). Then in the interior of by Proposition 5.2. Moreover, implies that and thus by assumption . Since has a global maximum at , we obtain that, for some constant , hold for all , . Therefore, along solutions that stay in the largest invariant subset of , the following relation holds Using the -equation in system (4.2), we have This equation holds only at . Substituting into (5.6), we obtain The monotonicity assumption implies that for all . Thus, the largest invariant set in the set is the singleton . By LaSalle’s Invariance Principle [60], is globally asymptotically stable in the interior of . The global stability also implies that the endemic equilibrium is unique.  □

We note that in the case when all newly infected individuals go into the first state, namely, , Theorem 5.1, Theorem 5.2 give the global-stability results in [9].

We also note that, if is a nonconstant function, the global stability of when remains largely open, even for the case when is monotone. For certain non-monotone functions , it was shown in [49] that multiple endemic equilibria can exist and periodic oscillations can occur, and the global stability results are not expected to be valid.

Numerical investigations and discussions

Results from numerical simulations are provided in this section to demonstrate our theoretical results. Several observations are made regarding the impacts of state structures on the dynamics. As we will show, some of the observations shed interesting lights on the impacts of state structure and raise interesting questions for further theoretical studies.

To focus our study on the state structure, we choose a simple form for functions , , and in model (4.2). We let and consider a special star-shaped directed graph for cross-transmissions among the 5 states, as shown in Fig. 6-(a). The corresponding weight matrix is given by is irreducible. The simplified model is described by the following systems of differential equations: It can be verified that all assumptions - in Section 2, - and - in the previous sections are automatically satisfied by model (6.1). From the results in Section 5 we know that the disease-free equilibrium is globally asymptotically stable if , and that, when , there exists a unique positive endemic equilibrium and it is globally asymptotically stable.

Fig. 6

Directed graphs for transfers among disease states, transmission coefficients of disease states, fractions of new infections into disease states, and distributions of disease prevalence among states at the endemic equilibrium.

Model (6.1) can be regarded as an approximation for the transmission dynamics of HIV infection with treatment among an adult population. We have chosen the following parameter values based on those in [4], [65]: The initial populations are chosen as . We have also selected two sets of transmission coefficients: for which the basic reproduction number ; and for which the basic reproduction number is . Results in Theorems 5.1 and 5.2 are illustrated in Fig. 5.

Fig. 5

Simulations for model (6.1) are shown to demonstrate global threshold results in Theorem 5.1, Theorem 5.2. Parameter values are given in (6.2)–(6.4).

While the parameter values in (6.2)–(6.4) were not fitted from clinical data, they are chosen in comparable scale or order of magnitude to those parameter values in two previous HIV studies [4] and [65]. In [4], detailed discussions were given on the selection of parameter values from the HIV literature for a -stage progression model of HIV, which helped us to choose values for state-transfer parameters that are biologically plausible for HIV progression. In [65], parameter values were fitted for a simple transmission model for the HIV for a population in a remote Chinese village, which helped us to choose values for demographic parameters and transmission coefficients. Results in the discussion are not only of general theoretical interest, but also relevant to HIV progression and control.

The impact of parameters on the distribution of disease prevalence among states

In this section, we discuss the distribution of number of infected individuals among different disease states when the dynamics have stabilized at the endemic equilibrium, as well as important factors that impact this distribution. This is a complicated problem on its own. Using numerical simulations, we provide some interesting observations. Let be the endemic equilibrium, we are interested in factors that influence the distribution of disease prevalence at state given by . A specific question of interest is the following: what factor (or factors) determines that a particular is the highest among all ? or determines the order among ? Important factors we have investigated include transmission coefficients , fraction of disease incidences in state , and transfer rates among disease states.

Observation I: Varying transmission coefficients do not alter the order among the disease prevalences , .

This is intuitively clear since directly influence the disease incidences, which is given in the following form, for the state : We expect that will influence the overall values of all but will not alter the relative order among the disease prevalence . In another word, the highest may not produce the largest .

In Fig. 6, the transmission coefficients at each state are given as and plotted in Fig. 6-(b). The disease prevalence at among states is plotted in Fig. 6-(d). We see that the second state has the highest prevalence , while it has the lowest transmission coefficient .

Observation II. For transfers among states, the state with the highest in-degree or out-degree may not have the highest disease prevalence.

The impacts of transfer rates among states on the distribution of are highly complex and is an interesting topic worthy of further investigation. In the example shown Fig. 6, the state-transfer network is given in Fig. 6-(a) together with the weights. It is clear from the directed graph that the first state (vertex) has the largest number of edges coming in and going out. Its in-degree and out-degree in the weighted directed graph are both the highest among all states. However, the highest disease prevalence occurred at the second state , even when .

Observation III. Proportions of disease incidence among states have a direct influence on distribution of disease prevalence.

From (6.5) we see that has the direct influence on the number of new infections going into state . We expect that distribution of will mimic that of . As shown in Fig. 6-(c) and Fig. 6-(d), this is indeed the case.

Summarizing our observations, the parameter , which is the proportion of new infections entering the state , has the most important impact on deciding which state will have the highest disease prevalence. The biological significance of this observation is that it is important for clinical researchers to measure proportions of newly infected HIV patients who are in different HIV states. It is known that some patients whose HIV infection progresses much faster than normal patients, the so-called ‘fast progressors’, while some patients whose HIV infection progresses much slower, the ‘nonprogressors’, and others can manage to control the level HIV viral load without treatment treatments, the ‘elite controllers’ [28], [29]. Out results show that it is important for the HIV control to estimate the proportions of these ‘abnormal’ patients.

Impacts of transfer among states on the basic reproduction number

The basic reproduction number determines whether an infectious disease can be effectively controlled. Many disease interventions are aiming at reducing . It can be seen in the expression of in (4.12) that reducing the transmission coefficients will lower the value of , which is also biologically intuitive. However, the impacts of changing transfer rates among infection states are not so obvious. The transfer rates can be altered through medical interventions. For instance, antiretroviral therapies for HIV infection can decrease and suppress the level of viral load and allow the immune system to recover. This will result in great health benefits for the infected individual. However, if the viral load is not fully suppressed, the infected individual can still be infectious. A healthier condition may allow an individual to resort to riskier sexual behaviours and may make the individual more infectious. While disease amelioration directly benefits an infected individual, its impacts on and the overall benefits at the community level are far from simple.

When there are two infected states and , when other parameters are fixed in model (6.1), it can be rigorously shown that an increase of or a decrease of can lead to a reduction of if , where , which are the basic reproduction number of the th state when it is isolated from other states.

To illustrate the complexity concerning the impacts of state-transfer on basic reproduction number in model (6.1), we choose the parameters as given in (6.2), (6.6). Then , , , , , thus the basic reproduction number , by (4.12).

Simulation results are shown in Table 1. We can see that an increase in will decease from its baseline value 6.0074; while an increase of will increase from the baseline value.

Table 1

Different effects on of increasing state transfer rates. Up arrows indicate increases from the baseline values and down arrows indicate decreases from the baseline values.

	(δ_21,δ_31,δ_41,δ_51,δ_12,δ_13,δ_14,δ_15)	R_0
Baseline	(0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3)	6.0074 ↓
δ_21↑	(0.25, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3)	6.0216 ↑
δ_31↑	(0.2, 0.25, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3)	6.0009 ↓
δ_41↑	(0.2, 0.2, 0.25, 0.2, 0.3, 0.3, 0.3, 0.3)	5.9568 ↓
δ_51↑	(0.2, 0.2, 0.2, 0.25, 0.3, 0.3, 0.3, 0.3)	5.7936 ↓
δ_12↑	(0.2, 0.2, 0.2, 0.2, 0.35, 0.3, 0.3, 0.3)	5.9931 ↓
δ_13↑	(0.2, 0.2, 0.2, 0.2, 0.3, 0.35, 0.3, 0.3)	6.0110 ↑
δ_14↑	(0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.35, 0.3)	6.0265 ↑
δ_15↑	(0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.35)	6.0674 ↑

These observations strongly suggest that further theoretical investigations of the impacts of state transfers on the overall dynamics are warranted.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.