Qualitative analysis of generalized multistage epidemic model with immigration

1.   Introduction

Mathematical models have revolutionized our understanding of the spread and control of infectious diseases. By analyzing complex interactions between different factors such as population dynamics, environmental factors and infectious agents, these models have helped public health officials make informed decisions that have saved countless lives. Diseases that progress through multiple stages also present unique challenges in epidemiology. Infectious diseases, such as HIV or tuberculosis, have distinct phases of disease progression, each with its own set of symptoms and transmission characteristics. For example, in HIV infection, the initial acute phase is followed by a long period of asymptomatic infection and, finally, the symptomatic phase leading to AIDS. Similarly, tuberculosis progresses from latent infection to active disease, with varying levels of contagiousness and symptom severity. Accurately identifying the stage of disease progression is critical for the design of appropriate public health interventions, as different stages require different strategies for control and prevention ^[1,2]. Among all of the factors that affect the spread of a disease, migration and immigration are important aspects to consider, especially in diseases where asymptomatic individuals or carriers can transmit it, as their traceability becomes very challenging since most diseases are identified based on symptomatic presentation. Thus, it is crucial to take migration and immigration processes into account in predictive models, as they can result in a more accurate prediction of the dynamics ^[3,4]. By studying the complex dynamics of disease progression in the case of infectious diseases, we can improve our understanding of the underlying biological processes, develop new interventions and ultimately reduce the burden of these diseases on global public health ^[5,6,7]. Therefore, incorporating migration and immigration processes into multistage mathematical models is crucial to comprehensively understanding complex population dynamics. In this paper, we present a mathematical model considering an immigration process that incorporates multiple infection states, including latent state ($E$), which are considered as infectious without showing any symptoms and infectious states ($I_1, I_2, \ldots, I_n$), as well as a recovered state ($R$), with the aim to study the global stability of the equilibrium points. The incidence rate is given by $f(S)g(E)+f(S)\sum_{i = 1}^{n}h_i(I_i)$, where the function $f(S)$ is the contact function and the functions $g(E)$ and $h_i(I_i)$, for $i = 1, \ldots, n$, are the force of infection, while the dead rate for all states of the model are proportional to $\mu_1\sigma_1(S), \mu_2\sigma_2(E), \mu_3\sigma_3(I_1), \ldots, \mu_{n+2}\sigma_{n+2}(I_n)$ and $\mu_{n+3}\sigma_{n+3}(R)$, respectively; and the rates at which each state evolves to an infectious or recovered state are $\alpha_2\sigma(E), \alpha_3\sigma_3(I_1), \ldots, \alpha_{n+1}\sigma_{n+1}(I_{n-1})$, and $\alpha_{n+2}\sigma_{n+2}(I_n)$, respectively.

Therefore, the model is given by the following system of differential equations:

where $\Lambda_1$, $\Lambda_2, \ldots, \Lambda_{n+3}$ represents the constant flow of new members into each compartment, respectively.

Multistage models with a bilinear or nonlinear incidence rate have been studied; for example, in ^[8], a multistage model with amelioration is considered to study the disease progression of HIV/AIDS, as well as the global stability of the equilibrium points. The incidence rate used is the bilinear rate known as the pseudo mass action law of the form $\sum_{m = 1}^{r}c\beta_m\frac{I_m}{N}S$. A generalization of the previous result is given in ^[9], the incidence rate was $f(N)\sum_{i = 1}^{n}\lambda_iI_iS$, where $f$ is a function that depends on the density of the population. In ^[10], the global stability of a general multistage model was proved; the term for the incidence rate was $\sum_{j = 1}^{n}f(N)g_j(S, I_j)$, and the death rate functions and the transfer rate functions were different and permits amelioration or immunity restoration. In ^[11], an extension of a previously proposed model is introduced to investigate global stability in a general cholera model. The model incorporated a nonlinear incidence rate given by the expression $\sum_{j = 1}^{n}D(N)f_j(S, I_j)+\sum_{j = 1}^{m}g_j(S, W_j)$, where $S$ and $I_j$ possess commonly established significance and $W_j$ denotes the number of pathogens shed by individuals. In ^[12], a model with imperfect vaccine and multistage behavior is studied; this model is similar to ^[10]; because the vaccinated people can infect others and they will have an influence on the dynamics of the disease, the incidence rate used was $\sum_{j = 1}^{n}y(N)h_j(S, I_j)+\sum_{j = 1}^{m}y(N)g_j(V, I_j)$. None of the previous multistage models consider immigration terms and nonlinear incidence rates, as our model does. Models that include immigration and, a nonlinear incidence rate have been studied in ^{[13,14,15,16,17,18]} but none of them are multistage. The novelty of our model is that it considers general incidence rates and also takes into account migration processes; besides, it generalizes the work presented in ^[19].

This paper is organized in the following manner. In Section 2, we impose conditions on the functions included in the model and prove the existence of a single equilibrium point. The global stability of a unique equilibrium point is proved in Section 3. In Section 4, we present numerical simulations to illustrate our main result. Finally, in Section 5, we discuss our results and provide some further extensions of the model.

2.   Model analysis

To investigate the model dynamics, some conditions over the functions are required:

${\rm{I}}).$ $f, \ g, \ h_i, \ \sigma_j$ are strictly increasing functions in $[0, +\infty)$ and $f(0) = g(0) = h_i(0) = \sigma_j(0) = 0$ for $i = 1, \ldots, n$ and $j = 1, \ldots, n+3$.

${\rm{II}}).$ There are positive constants $k_1$, $k_2, \ldots, k_{n+3}$ such that $\sigma_1(S)\geq k_1 S$, $\sigma_2(E)\geq k_2 E$, $\sigma_{i+2}(I_i)\geq k_{i+2} I_i$ for $i = 1, \ldots, n$, $\sigma_{n+3}(R)\geq k_{n+3} R$.

${\rm{III}}).$ $\dfrac{g(E)}{\sigma_2(E)}$ and $\dfrac{h_i(I_i)}{\sigma_{i+2}(I_i)}$ for $i = 1, \ldots, n$ are non-increasing functions in $(0, +\infty)$.

The hypotheses Ⅰ and Ⅲ are necessary conditions for achieving global stability, and Ⅱ is a necessary condition for verifying the existence of an invariant set. Moreover, conditions Ⅰ and Ⅱ allow us to work with positive quantities, i.e., they are biological conditions. Condition Ⅲ can be interpreted as a saturation of the force of infection with respect to the increase in infectious individuals, as this quotient can be a constant, a decreasing function or a combination of these. In fact, $\lim\limits_{I\to \infty} \frac{h(I)}{\sigma(I)} = \lim\limits_{I\to \infty} \frac{h(I)}{I}\frac{I}{\sigma(I)}$ is either a constant or zero. It should be noted that, according to condition Ⅱ, $\frac{I}{\sigma(I)}$ is bounded, which implies that $h(I)$ decreases as $I$ increases. The saturation effect in the force of infection was first discussed in ^[20].

The feasible set for the model (1.1) is given by

where $\Lambda = \Lambda_1+\Lambda_2+\dfrac{\mu_2}{2\alpha_2}\Lambda_3+\sum_{i = 2}^{n}\left(\frac{\Lambda_{i+2}}{2^i}\prod_{j = 2}^{i+2}\frac{\mu_j}{\alpha_j}\right)$, $\delta = \min\bigg\{\mu_1k_1, \frac{\mu_2}{2}k_2, \frac{\mu_3}{2}k_3, \ldots, \frac{\mu_{n+1}}{2}k_{n+1}, \mu_{n+2}k_{n+2}\bigg\}$ and $k_1$, $k_2, \ldots, k_{n+2}$, as defined in condition Ⅱ.

Proposition 1. The set $\Omega$ is positively invariant with respect to system (1.1).

Proof. Let $\Theta(t)$ be the function defined by

where $S$, $E$, $I_1, \ldots, I_n$ are the solution of model (1.1), and let $\left(S(0), E(0), I_1(0), \ldots, I_{n}(0)\right)\in \Omega$ be the initial condition of the system (1.1). Taking the derivative of $\Theta$ with respect to $t$, we have

By the condition Ⅱ, we have that $-\mu_1\sigma_1(S)\leq -\mu_1k_1S$, $-\frac{\mu_2}{2}\sigma_2(E)\leq -\frac{\mu_2}{2}k_2E, \ldots$, $-\dfrac{\mu_2\mu_3\cdots\mu_{n+1}}{2^n\alpha_2\alpha_3\cdots\alpha_{n+1}}\mu_{n+2}\sigma_{n+2}(I_{n})$$\leq -\dfrac{\mu_2\mu_3\cdots\mu_{n+1}}{2^n\alpha_2\alpha_3\cdots\alpha_{n+1}}\mu_{n+2}k_{n+2}I_{n}$, which implies that

By taking $\delta$ of the form $\delta = \min\bigg\{\mu_1k_1, \frac{\mu_2}{2}k_2, \frac{\mu_3}{2}k_3, \ldots, \frac{\mu_{n+1}}{2}k_{n+1}, \mu_{n+2}k_{n+2}\bigg\}$, we obtain

It follows that $\frac{d \Theta}{d t}\leq 0$ if $\Theta \geq \frac{\Lambda}{\delta}$. Besides, we have

In particular, $\Theta\leq \frac{\Lambda}{\delta}$ if $\Theta(0)\leq \frac{\Lambda}{\delta}$. Therefore, the set $\Omega$ is positively invariant. In addition, if $\Theta(0) > \frac{\Lambda}{\delta}$, then either the solutions enters into the set $\Omega$ infinite times or $\Theta (t)$ approaches $\frac{\Lambda}{\delta}$ asymptotically. Hence, the set $\Omega$ attracts all solutions in $\mathbb{R}_+^{n+2}$.

□

Proposition 2. There exists an endemic equilibrium $(S^*, E^*, I_1^*, I_2^*, \ldots, I_n^*)$ of the system (1.1).

Proof. To obtain the equilibrium points of the model (1.1), we need to solve the following system of equations:

From the equations (2.2) and (2.3), we have

in this way,

and from the equation (2.4), we obtain

Continuing in this way with the equation (2.5), for $i = 2, \ldots, n$, we get

Since $\sigma_1, \sigma_2, \ldots, \sigma_{n+2}$ are strictly increasing functions, we can solve the equations (2.6)-(2.8) in the form $E = \phi_0(S), I_1 = \phi_1(S), \ldots, I_n = \phi_{n+2}(S)$. Now, let $\varphi$ be the function defined by

which depends just on $S$. If we solve the equation $\varphi(S) = 0$ for some $S^{*}$, we get that $E^* = \phi_0(S^*), I_1^* = \phi_1(S^*), \ldots, I_n = \phi_{n+2}(S^*)$, and in this way, we obtain the equilibrium point.

We notice that, in the equation (2.6), when

we have that $\sigma_2(E) = 0$. Now remembering that $\sigma_1(0) = 0$ and $\sigma_1$ is a strictly increasing function, there exists $\bar{S}$ such that $\sigma_1 (\bar{S}) = (\Lambda_1+\Lambda_2)/\mu_1$, and, for the same reason, we have that $\sigma_2(E) = 0 \Leftrightarrow E = 0$ when $S = \bar{S}$. Similarly, from equation (2.7), there exists $\bar{I_1}$ such that $\sigma_3(\bar{I_1}) = \Lambda_3/\mu_3$ when $S = \bar{S}$, and from equations (2.8), there exists $\bar{I_i}$ such that $\sigma_{i+2}(\bar{I_i}) = \frac{\Lambda_{i+2}}{\mu_{i+2}} + \sum_{k = 3}^{i+1}\left(\frac{\Lambda_k}{\mu_k}\prod_{j = k}^{i+1}\frac{\alpha_j}{\mu_{j+1}}\right)$ when $S = \bar{S}$. Therefore, we want to find a root of equation (2.9) in the interval $(0, \bar{S})$. To do this, observe that $\varphi(0) = \Lambda_1 > 0$ and

This means that, for the continuity of $\varphi$, there exists $S^{*} \in (0, \bar{S})$ such that $\varphi(S^{*}) = 0$. In conclusion, we show that there is a single equilibrium point $(S^*, E^*, I_1^*, I_2^*, \ldots, I_n^*)$. □

3.   Global stability

In this section, we prove the stability of the unique equilibrium point by using a Lyapunov function. First, we prove a proposition that will guarantee that the Lyapunov function is positive and only zero at the endemic equilibrium point.

Proposition 3. If $\phi$ is a continuous increasing function in $(0, \infty)$, the function

is positive for $x > 0$ and $\varPsi(x) = 0$ just for $x = x^*$.

Proof. First, we start showing that, in $x^* > 0$, the function $\varPsi$ has a minimum. In fact, $\varPsi'(x) = 1-\frac{\phi(x^*)}{\phi(x)}$; since $\phi$ is an increasing function, we have that $\varPsi'(x)$ is negative if $x < x^*$ and positive if $x > x^*$. Thus, $\varPsi$ attains its minimum value at $x^*$. Finally, we can see that $\varPsi(x^*) = 0$; thus, $\varPsi(x) > 0$ for $x > 0$ ($x\neq x^*$). □

For the proof of the next theorem, we omit the equation for recovery state R because it does not appear in the other equations.

Theorem 4. The equilibrium point $(S^*, E^*, I_1^*, I_2^*, \ldots, I_n^*)$ is globally asymptotically stable.

Proof. Let $L$ be the Lyapunov function defined by

where $a_1 = \frac{f(S^*)\sum_{i = 1}^{n}h_i(I_i^*)}{\alpha_2\sigma_2(E^*)}$, $a_i = \frac{f(S^*)\sum_{j = i}^{n}h_j(I_j^*)}{\alpha_{i+1}\sigma_{i+1}(I_{i-1}^*)}$. By Proposition 3, we have that $L > 0$ and $L = 0$ just in $(S^*, E^*, I_1^*, I_2^*, \ldots, I_n^*)$. This Lyapunov function was first proposed in ^[21]. The orbital derivative of $L$ is

From the equilibrium equations, we have

Puting these equations in the orbital derivative, we obtain

rearranging and grouping terms, we have that

Adding the terms that have in common $f(S^{*})g(E^{*})$ and those that have $f(S^{*})\sum_{i = 1}^{n}h_i(I_i^{*})$, we get

From the equations for $a_1$ and $a_i$, we have that

□

Now, adding the like terms, we obtain

After respectively adding and subtracting the terms $f(S^{*})g(E^{*})$, $\frac{\sigma_2(E)}{\sigma_2(E^{*})}\frac{g(E^{*})}{g(E)}$ and $f(S^{*})h_i(I_i^{*})$, $\frac{\sigma_{i+2}(I_i)}{\sigma_{i+2}(I_i^*)}\frac{h_i(I_i^{*})}{h_i(I_i)}$ for each $i = 1, \ldots, n$ and regrouping some expressions, we get

Now, we see that

for $i = 1\ldots, n.$ From this, we can write the last equations for $\dot{L}$ as

We know that the geometric mean is lower than or equal to the arithmetic mean, i.e., $\sqrt[n]{x_1x_2\cdots x_n}\leq \frac{x_1+x_2+\cdots+x_n}{n}$ or $n\sqrt[n]{x_1x_2\cdots x_n}\leq x_1+x_2+\cdots+x_n$. Hence, we obtain

By the conditions Ⅰ and Ⅲ, and for $i = 1, \ldots, n$, we have

Therefore, $\dot{L} \leq 0$ for all $(S, E, I_1, I_2, \ldots, I_n)$ and $\dot{L} = 0$ if and only if $(S, E, I_1, I_2, \ldots, I_n) = (S^*, E^*, I_1^*, I_2^*, \ldots, I_n^*)$. So, the equilibrium point is globally asymptotically stable.

4.   Numerical simulation

Here, we introduce an example of HIV/AIDS transmission dynamics to show the theoretical results. The population is divided into four stages of disease progression, susceptible to HIV infection ($S$), HIV-positive individuals in the acute HIV infection stage ($E$), HIV-positive individuals in the chronic HIV infection stage ($I$) and individuals with full-blown AIDS ($A$). We divided the chronic HIV infection stage into two groups ($I_1$ and $I_2$). Besides, we assume that people in compartment $I_2$ are more infectious than those in $I_1$.

Susceptible individuals can be infected through contact with HIV-positive individuals. The susceptible individuals that were infected go to the acute HIV infection stage ($E$). After a period $\alpha_2^{-1}$, they progress to the first chronic stage of the disease ($I_1$). Individuals in the first chronic stage progress to the second chronic phase after a period $\alpha_3^{-1}$. The individuals in the $I_2$ compartment progress to full-blown AIDS after a period $\alpha_4^{-1}$. Here, all compartments have a recruitment rate $\Lambda$ and natural mortality rate $\mu$. Also, individuals with full-blown AIDS have an additional mortality rate due to the disease.

Therefore, we obtain the following system of differential equations:

where $h_1 (I_1) = m_1 I_1 / (1+I_1)$ and $h_2(I_2) = m_2 I_2 / (1 + I_2)$. Table 1 shows the model parameters and their description. Figures 1 and 2 show the dynamics of infectious individuals for the scenarios described before.

To simulate some epidemiological scenarios, we assumed four scenarios with different initial conditions for the infected individuals. In scenario 1, we assumed that $E(0) = I_1(0) = I_2(0) = 1$ and $A(0) = 10$. For the scenario 2, $E(0) = I_1(0) = I_2(0) = 10$ and $A(0) = 100$; in scenario 3, we suppose that $E(0) = I_1(0) = I_2(0) = 50$ and $A(0) = 500$; finally for the last scenario, we have that $E(0) = I_1(0) = I_2(0) = 100$ and $A(0) = 1000$. For susceptible individuals, we set it as 100 in all scenarios. Figure 2 shows the dynamics of individuals with full-blown AIDS according to the scenarios of the initial conditions.

As expected from the global stability, all solutions converge to the respective equilibrium point coordinate.

5.   Discussion

In this paper, we analyzed a multistage mathematical model that includes a general incidence function, death rate functions and immigration in all stages of the model. Our model has only one equilibrium point due to the constant rate of immigration in all populations that transmit the disease. We prove that this equilibrium point is globally asymptotically stable by using an appropriate Lyapunov function and considering sufficient conditions for the functions involved in the model. Our results provide a foundation for creating models that can help us understand how different incidence rates affect the spread of diseases that have multiple stages in the presence of migration and immigration. Additionally, our allows us to explore what might happen in scenarios that do not match our initial assumptions. In the context of the model, migration or immigration terms can be interpreted as vertical transmission. There is also the possibility of investigating non-constant migration terms; we believe that the global stability theorem presented here would need little change, that the only thing left to prove would be the existence of an endemic equilibrium. One limitation of the present model is that it considers the process of immigration/migration as a constant influx, but rarely is this true in a real context. Generally, this process occurs in a discontinuous way. Another limitation is the absence of $R_0$, which results in only one endemic equilibrium point. If we set the migration terms to zero, we can recover the disease-free equilibrium point and, therefore, the possibility of finding the $R_0$ threshold that allows for the eradication of the disease, as mentioned and discussed in ^[26]. This suggests that stopping migration or immigration is an effective measure to try to stop the transmission of a disease. This can also be interpreted as quarantine measures taken to halt a disease.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Miller Cerón Gómez acknowledges the financial support of grant 2472 from Vicerrectoria de Investigaciones e Interacción Social (Viis-Udenar) acuerdo 21 from march 8, de 2022. Felipe Rubio acknowledges support from the grant CEX2018-000806-S funded by MCIN/AEI/ 10.13039/501100011033, as well as support from the Generalitat de Catalunya through the CERCA Program.

Conflict of interest

The authors declare that there is no conflict of interest.