Dynamics of an epidemic model with general incidence rate dependent on a class of disease-related contact functions

1.   Introduction

During disease transmission, most pathogens are transmitted by person-to-person contact. Therefore, establishing mathematical models to study infectious diseases, especially infection links, is a very important research goal. The rate at which new compartments of infected persons are added to a model is usually referred to as the incidence of disease. Different infectious diseases have different transmission routes, so the selection of disease incidence rate in the infectious disease model cannot be studied in the same form. Many infectious diseases break out repeatedly in large numbers of people. Mechanisms that can explain this phenomenon include seasonal change ^[1,2], transmission rate, infection or the time delay in the process of recovery ^[3] and heterostructure of host population ^[4], etc. As research continues, many scholars have focused on no periodicity of bilinear transmission research ^{[5,6,7,8,9,10]}. Through the study of different nonlinear incidence, more abundant dynamics are obtained, which can be used to describe and explain some periodic epidemic phenomena in reality.

The contact form within the population has an important influence on the disease transmission process ^[11]. In ^[12], a mixed proportion of infection incidence was proposed and the number of effective contacts per population is considered to be a fixed constant number. However, with the change of population within the population and the development of society, the contact pattern between populations will also change. It is not very effective to take the number as a constant. We consider it as a function dependent on population size. In a population, it is assumed that the members are evenly distributed and the number of contacts between each member in unit time is the same. Because the population is homogeneous, we take the contact function as $T(N)$, where $N$ is the total number of the population.

For the effective contact function $T(N)$, we propose assumptions (A1)–(A3) described below. Since the function represents the number of effective contacts per unit time of each member of the population, it is non-negative. Thus, the first assumption is as follows:

(A1) $\forall N\in[0, +\infty), 0 = T(0)\leq T(N) < +\infty.$

Considering that the effective contact times of each member of the crowd in unit time will not decrease with the increase of the population, $T(N)$ is a monotone undecreasing function. Thus, we give the second assumption:

(A2) $\forall N\in [0, +\infty), T'(0) > 0, T'(N)\geq 0.$

We make the third assumption of $T(N)$ regarding the growth rate $N$ as follows:

(A3) $\forall N\in [0, +\infty), T'(N)N-T(N)\leq 0.$

Table 1 shows two special cases of the effective contact function $T(N)$ that satisfy assumptions (A1)–(A3). These two types of contact functions correspond exactly to the standard incidence rate and the bilinear incidence rate.

The disease we are concerned with may affect the contact times of the infected person (measures taken by the infected person or the disease itself will cause a behavior change of the infected person). At the same time, we regard susceptible persons $S$ as a reference group, and record the ratio of the number of contacts of infected person $I$ to $S$ as $\eta$, $\eta\in[0, \infty)$. Also, we say that $\eta$ is an intensity factor of disease affecting the behavior of the infected person.

Considering $\eta$ as the observed variable, we have the following description

1. $\eta = 0$. $ This can be regarded as the disease can not be transmitted or infected through regular contact. Infected persons are completely isolated;

2. $\eta\in(0, 1)$. $ It can be seen that when the susceptible person gets sick, the disease will reduce (or the infected person has taken some measures) the number of contacts of the infected person compared to the susceptible person;

3. $\eta = 1$. The disease has no effect on the effective contact function $T(N)$;

4. $\eta\in(1, +\infty)$. $ It can be seen that when the susceptible person gets sick, the disease will increase the number of contacts of the infected person compared to the susceptible person.

For the intensity factor, we propose assumptions (A4)–(A5) described below.

(A4) $\eta\in\left(0, 1 \right].$

(A5) $\eta\in\left(0, 1+\cfrac{\mu}{d} \right].$

Based on the above discussion, we have established the following system (1.1).

where $S(t)$ is the number of susceptible persons in the population. $I(t)$ is the number of infected persons in the population. $\Lambda$ is the number of constant increase (birth and immigration) of the population in unit time. $d$ is migration (death and emigration) rate of population. $\beta$ is the probability that a susceptible person will get sick when an effective contact occurs between a susceptible person and an infected person. $\mu$ is the case fatality rate. $\gamma$ is the recovery rate of the disease. $\eta$ is the intensity factor.

By using $N(t) = S(t)+I(t)$, we can easily obtain a system equivalent to system (1.1).

where $N(t)$ is the total number of the population. The other parameters are the same as the system (1.1).

Let

From a biological perspective, it is clear that the dynamics of system (1.1) is a dynamical system in $\mathcal{X}_{0}$. We obtain the dynamics of (1.1) by studying the dynamics of the equivalent system (1.2) in $\mathcal{Y}_{0}$.

In Section 2, we analyze the dynamics of system (1.1). First, we calculate equilibra and the basic reproduction number. Then, we study the asymptotic stability of equilibria. We show that there exists a limit cycle. In Section 3, we provide some numerical simulation results which show the existence of a Hopf bifucation. In Section 4, we give some conclusions, discuss the effective contact function and propose future research directions.

2.   Dynamics of the model

In this section, we study the dynamics of the model (1.1).

2.1.   Equilibra and basic reproduction number

Simple calculations quickly show that system (1.2) admits a disease free equilibrium (DFE) $E_0 = (\Lambda/d, 0)$.

For system (1.2), we use the next generation matrix method ^[14] to calculate the threshold $\mathcal{R}_0$. Let

By calculating, we have

So, we define the threshold of system (1.2) as

It is simple to get that $E_0$ also is the DFE of the system (1.1), and $\mathcal{R}_0$ also is the threshold of (1.1).

Remark 2.1. $\beta$ is the probability that a susceptible person will get sick during one active contact with an infected person. $\eta T(N_0)$ is the number of effective contacts with susceptible persons in unit time when an infected person is introduced into a DFE state. $1/(d+\mu+\gamma)$ is the average duration of illness for an infected person. Therefore, the threshold $\mathcal{R}_0$ is the basic reproduction number of systems (1.1) and (1.2).

Theorem 2.1. Assume (A1) and (A2) hold. Then we have the following two conclusions:

(Ⅰ) The DFE $E_0$ always exists for system (1.1). Furthermore, if $\mathcal{R}_0\leq 1$, assume (A3) or (A5) holds, then the DFE $E_0$ is the unique equilibrium for system (1.1).

(Ⅱ) If $\mathcal{R}_0 > 1$, then system (1.1) has at least one endemic equilibrium (EE). Furthermore, if $\mathcal{R}_0 > 1$, assume (A3) or (A5) holds, then system (1.1) has a unique EE, $E^*$.

Proof. Obviously, the DFE $E_0$ always exists. We let the endemic equilibria of (1.1), $E^* = (S^*, I^*)$, where $S^* = \Lambda/d-(1+\mu/d)I^*$, and $I^*$ satisfies the following equation

Let $I_m = {\Lambda}/{(d+\mu)}$. For $S^* > 0$, we have $I^*\in[0, I_m)$. In $[0, I_m]$, we let

About the continuous function $H(I)$, we have the following conclusions:

where ${N}_{_{{\eta}}} = \Lambda/d+(\eta-1-\mu/d)I > 0, \; I\in(0, I_m).$

If $\mathcal{R}_0 > 1$, we have $H(0) > 0$. If $\mathcal{R}_0\leq 1$, we have $H(0)\leq 0$. Assume (A3) or (A5) holds. Then, ${\mathrm{d} {H}}/{\mathrm{d} {I}} < 0$.

By the properties of the continuous function $H(I)$, we have the following conclusions:

(a) If $\mathcal{R}_0 > 1$, then $H(I) = 0$ has at least one root in $(0, I_m)$.

(b) $\mathcal{R}_0\leq 1$. If assume (A3) or (A5) holds, then $H(I) = 0$ has no roots in $(0, I_m)$.

(c) $\mathcal{R}_0 > 1$. If assume (A3) or (A5) holds, then $H(I) = 0$ has a unique root, $I^*$, in $(0, I_m)$.

Therefore, (Ⅰ) and (Ⅱ) are proved. □

It is simple to get that if $E^* = (S^*, I^*)$ is the EE of system (1.1), then $\bar{E}^* = (S^*+I^*, I^*)$ is the EE of system (1.2).

2.2.   Local stability of equilibria

We give the local stability of all equilibria of (1.1).

Theorem 2.2. Assume (A1) and (A3) hold. The following statements are valid for system (1.1).

(Ⅰ) If $\mathcal{R}_0 < 1$, then $E_0$ is local asymptotically stable. If $\mathcal{R}_0 > 1$, then $E_0$ is unstable (Saddle).

(Ⅱ) If $\mathcal{R}_0 > 1$, assume (A3) or (A4) holds, then $E^*$ is local asymptotically stable.

Proof. In order to analyze the local stability of the DFE $E_0 = (\Lambda/d, 0)$ of (1.1), by using Eq (2.1), we linearize system (1.2) at the DFE $E_0$,

The Jacobian matrix at $E_0$ is

It is easy get that $\lambda_1 = -d$ and $\lambda_2 = (d+\mu+\gamma)(\mathcal{R}_0-1)$ are the characteristic roots of $J(E_0)$. It follows that all eigenvalues of $J(E_0)$ have negative real parts if and only if $\mathcal{R}_0 < 1.$

Therefore, if $\mathcal{R}_0 < 1$, $E_0$ is local asymptotically stable; if $\mathcal{R}_0 > 1$, then $E_0$ is unstable. The proof of (Ⅰ) is completed.

The $E^* = (S^*, I^*)$ satisfies the following equation

We take $T^{*} = T(S^*+\eta I^*)$, $T^{'*} = T'(S^*+\eta I^*)$. In order to analyze the local stability of $E^*$ of (1.1), by using $E^*$ and (2.3), we linearize system (1.1) at $E^*$,

where

We obtain that the Jacobian matrix at $E^*$ is

Thus,

As $T(N)$ required in (A1)–(A3), we may take

Therefore, we have $a_{21} > 0$ and

.

If (A1)–(A3) hold, then we have $\mathrm{tr}(J(E^*)) < 0$ and $\det(J(E^*)) > 0$. According to the Routh-Hurwitz criteria, we get that, if the conditions (A1)–(A3) hold, $E^*$ is local asymptotically stable.

Of course, if we replace condition (A3) with (A4), we have

Then we easily get $\mathrm{tr}(J(E^*)) < 0$ and $\det(J(E^*)) > 0$. Thus, $E^*$ is local asymptotically stable. The proof of (Ⅱ) is completed. □

2.3.   Global stability of equilibria

Lemma 2.3. Assume (A1), (A3) hold. If $\mathcal{R}_0 > 1$, assume (A3) or (A4) holds, then $\bar{E}^*$ is globally asymptotically stable for system (1.2) in $\mathcal{Y}_{0}\backslash N$-axis.

Proof. For system (1.2), we take the Dulac function $B(N, I)$ as

Assume (A1), (A2) hold, we have

If $\mathcal{R}_0 > 1$, assume (A3) or (A4) holds, then $\bar{E}^*$ is the unique EE and local asymptotically stable. Also, we obtain that

By Theorems 2.1 and 2.2, system (1.2) has only one equilibrium, $E_0$, in $\mathcal{Y}_{0}\cap N\mbox{-axis}$. We have $E_0$ is unstable. According to the Bendixson-Dulac criteria, the conclusion in the theorem holds. □

So, we have the condition for the stability of equilibria of (1.1).

Theorem 2.4. Assume (A1), (A3) hold. The undermentioned statements are valid for system (1.1).

(I) The DFE $E_0$ always exists.

(II) If $\mathcal{R}_0 < 1$, then the DFE $E_0$ is locally asymptotically stable.

(III) If $\mathcal{R}_0 < 1$, assume (A3) or (A5) holds, then the DFE $E_0$ is the unique equilibrium and globally asymptotically stable in $\mathcal{X}_{0}$.

(IV) If $\mathcal{R}_0 > 1$, then the DFE $E_0$ is unstable.

(V) If $\mathcal{R}_0 > 1$, then system (1.1) has at least one EE.

(VI) If $\mathcal{R}_0 > 1$, assume (A3) holds, then the EE $E^*$ is the unique EE and globally asymptotically stable in $\mathcal{X}_{0} \backslash S$-axis.

(VII) If $\mathcal{R}_0 > 1$, assume (A5) holds, then the EE $E^*$ is the unique EE. Furthermore, assume (A4) holds, then $E^*$ is globally asymptotically stable in $\mathcal{X}_{0} \backslash S$-axis.

2.4.   Hopf bifurcation

From Theorem 2.4, we can find that complex dynamics in system (1.1) are only possible when the assumption (A3) of the contact function is not valid. We construct a series of simple polynomial functions to satisfy assumptions (A1) and (A2), but assumption (A3) does not hold. Then, we conduct numerical simulations to see if periodic solutions can be found. In the next section, we provide a contact function and corresponding parameter values that enable system (1.1) to discover periodic solutions. In order to analyze how parameter affects the direction of Hopf bifurcation, we may choose $\eta$ as a bifurcation parameter.

From Theorem 2.2, we obtain that the Jacobian matrix at $E^*$ is

The characteristic equation is

where $\mathrm{tr}(J(E^*))$ is as shown in (2.5), $\det(J(E^*))$ is as shown in (2.6). If there is an $E^*$ that makes $\mathrm{tr}(J(E^*)) = 0, \; \det(J(E^*)) > 0$, then there exists a periodic solution around the EE $E^*$ for $\Gamma_1 > 0$.

For $\eta$, the Hopf bifurcation point is a solution that holds at the following conditions

The Lyapunov number can be computed by^[15]

where

Now, we are able to state the results for the Hopf bifurcation.

Theorem 2.5. If $\mathcal{R}_0 > 1$, $\mathrm{tr}(J(E^*)) = 0, \; \det(J(E^*)) > 0$, then system (1.1) has a subcritical (supercritical) Hopf bifurcation at $E^*$ if $\Gamma_2 > 0$ ($\Gamma_2 < 0$).

3.   Numeric simulations

In this section, the influence of stability of equilibria of the model on the basic reproduction number $\mathcal{R}_0$, the contact function $T(N)$ and the intensity factor $\eta$ is investigated. The numerical simulations are performed using MATLAB (R2018a) software. The set of parameter values given in Table 2 and the initial conditions $(S(0), I(0)) = (1, 0.5)$ are used for the simulations.

From Theorems 2.1, 2.2 and 2.4, we find that the relationship between the threshold $\mathcal{R}_0$ and $1$ determines the equilibria of the system, but the stability of the equilibriia is also influenced by the selection of contact functions and parameter $\eta$. About contact function $T(N)$, we choose $T(N)$ as

where this $T(N)$ satisfies assumptions (A1) and (A2), and does not satisfy assumption (A3). Because these functions are too broad, the contact function (3.1) is just an example of the complex dynamic behavior where the system (1.1) will exhibit periodic solutions when the conditions in Theorem 2.4 are not met.

Figures 1 and 2 show the simulation of susceptible person $S(t)$ for the connect function $T(N) = 2$, $T(N) = 2N$, with $\eta$ being a variable, the given initial conditions and other parameter values in Table 2. From the figure, when $\mathcal{R}_0 < 1$, susceptible individuals ultimately reach their maximum scale. When $\mathcal{R}_0 > 1$, susceptible individuals ultimately reach a stable value. As $\mathcal{R}_0$ increases ($\mathcal{R}_0 > 1$), the final size of susceptible individuals also gradually decreases. These results are consistent with the conclusions in Theorem 2.4.

Figure 3 shows the simulation of susceptible person $S(t)$ for the connect function $T(N)$ in (3.1), with $\eta$ being a variable, the given initial conditions and other parameter values in Table 2. From the figure, it can be seen that as $\eta$ increases, the basic regeneration number $\mathcal{R}_0$ linearly increases. and the disease progresses from extinction to epidemic, while there may be periodic epidemics. And as $\eta$ increases, the final size of susceptible individuals continues to decline, and the scale of disease outbreaks is increasing.

Figure 4 is a bifurcation diagram drawn using Matcont, where $\eta = 5.5273$ and $\eta = 6.0396$ are Hopf bifurcation points, and their first Lyapunov exponents are $22.6683$ and $-23.1353$ respectively. It can be seen that when system (1.1) is at $\eta \in (5.5273, 6.0396),$ other parameter values in Table 2 and the connect function $T(N)$ using (3.1), the equilibrium of endemic disease disease is unstable, and there is a stable limit cycle, that is, the disease will eventually circulate periodically.

4.   Conclusions and discussion

In this paper, we have studied an SI model that has the general contact function and the impact of diseases on infected persons. We have provided some sufficient conditions for the ultimate extinction or epidemic of the disease. This model admits two types of equilibria and the global dynamics are determined by the basic reproduction number $\mathcal{R}_0$ and the intensity factor $\eta$. It must be pointed out that the results obtained in this paper only apply to simple cases where the contact function $T(N)$ satisfies the assumptions (A1) and (A2), as well as partially satisfies (A3).

From a biological perspective, assumptions (A1) and (A2) are usually valid. From Theorems 2.1, 2.2 and 2.4, if assumption (A3) holds, system (1.1) will only have two types of the equilibria, $E_0$ and $E^*$, and their global dynamics are determined by the basic reproduction number $\mathcal{R}_0$. If assumption (A3) does not hold, assumption (A5) and the relationship between the threshold $\mathcal{R}_0$ and $1$ are sufficient conditions for Theorem 2.1 (Theorems 2.2 and 2.4) to hold. The necessary and sufficient conditions for the theorems to hold true regarding the range of $\eta$ are work that needs to be further studied.

Numerical simulations were conducted to show that if the contact function $T(N)$ satisfies assumptions (A1) and (A2), and does not satisfy assumption (A3), system (1.1) may have a Hopf bifurcation, resulting in a stable periodic solution. Whether there will be other bifurcations is still a question that needs to be studied. This reveals that the epidemic of the disease cycle may be related to the mode of disease transmission and the impact of the disease on patient behavior patterns. When a disease is prevalent, behavioral constraints on potential carriers may be more efficient public health policies.

$T(N)$ describes the effective number of contacts between individuals within a population, and the population activity habits at different times and regions are usually not consistent. Therefore, when studying practical problems, it is necessary to analyze them in a specific manner and not use the same contact function $T(N)$ without thinking. The dynamic behaviors of models where the contact functions do not satisfy assumption (A3) require further study. Our main research in the future is to find a type of contact function that makes Hopf branches widely exist.

Use of AI tools declaration

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

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant No. 12071382). We thank the associate editor and the reviewers for their useful feedback that improved this paper.

Conflict of interest

The authors declared that they have no conflicts of interest to this work.