Effects of Nonmonotonic Functional Responses on a Disease Transmission Model: Modeling and Simulation.

In this article, a novel susceptible-infected-recovered epidemic model with nonmonotonic incidence and treatment rates is proposed and analyzed mathematically. The Monod-Haldane functional response is considered for nonmonotonic behavior of both incidence rate and treatment rate. The model analysis shows that the model has two equilibria which are named as disease-free equilibrium (DFE) and endemic equilibrium (EE). The stability analysis has been performed for the local and global behavior of the DFE and EE. With the help of the basic reproduction number R 0 , we investigate that DFE is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1 . The local stability of DFE at R 0 = 1 has been analyzed, and it is obtained that DFE exhibits a forward transcritical bifurcation. Further, we identify conditions for the existence of EE and show the local stability of EE under certain conditions. Moreover, the global stability behavior of DFE and EE has been investigated. Lastly, numerical simulations have been done in the support of our theoretical findings. © School of Mathematical Sciences, University of Science and Technology of China and Springer-Verlag GmbH Germany, part of Springer Nature 2021.

Introduction

Mathematical modeling is an effective tool to understand the dynamics of infectious diseases/epidemics from over a century. In mathematical epidemiology, many mathematical models have been proposed for disease transmission dynamics such as SIR [15, 7, 19, 10, 28, 8, 2], SIRS [26, 23], SEIR [6], and SVEIR [12] (where S, E, V, I, and R denote the susceptible, exposed, vaccinated, infected, and recovered individuals, respectively). In the disease transmission model, the rate of a new infection mainly depends on the incidence rate. The number of people who become infected per unit of time is defined as the incidence rate [10]. In the pioneering work of Kermack and McKendrick [14], the incidence rate is in the bilinear form , which was unrealistic for a large population because if the number of susceptibles increases, then the incidence rate also increases. So this unreasonable bilinear incidence rate is needed to be modified. Due to this unrealistic situation, many researchers [15, 7, 19, 8, 16, 1, 11, 17] have proposed various nonlinear functional forms such as Holling type II, Beddington–DeAngelis type, and Crowley–Martin type, since these functional forms incorporate inhibitions measures either by susceptibles or by both susceptibles and infectives, thus making the model more realistic for a large population. In 1986, to include the effect of behavioral changes, Liu et al. [20] proposed the nonlinear saturated incidence rate of the form:where measures the infection force of the disease, describes the inhibition effects from the behavioral changes of susceptibles individuals when the number of infected individuals increases, and are positive constants, and is a nonnegative constant [3, 23, 21]. Note that bilinear incidence is a particular case of (1.1) with and . The case when and in (1.1), i.e., , was proposed by Capasso and Serio [1] to describe the crowding effect in modeling the cholera epidemic in Bari in 1973. The function is also known as the Holling type II function. The case when and in (1.1), i.e., , represents a nonmonotonic incidence function which can be used to interpret the “psychological effects” [23]. The function is also known as the simplified Monod–Haldane (M–H) function. In M–H incidence function, for a large number of infectious individuals, the infection force may decrease as the number of infectious individuals increases, because in the presence of a large number of infectious individuals the population may tend to reduce the number of contacts per unit time, as seen with the spread of SARS [15].

It is well known that the treatment rate always plays a crucial role in preventing and controlling the spread of infectious diseases. We know that the treatment resources of any community are limited. Therefore, it is essential to choose an effective treatment rate to control the spread of the disease. In 2004, Wang and Ruan [25] considered an SIR epidemic model with a constant treatment rate given as follows:They performed a stability analysis and showed that their model exhibits various bifurcations. Later, in 2010, Zhang and Sao [27] modified the constant removal rate into the Holling type II:They analyzed the stability behavior of equilibria and investigated that the model exhibits various bifurcations such as Bogdanov–Takens bifurcation, Hopf bifurcation, and Homoclinic bifurcation. In 2013, Dubey et al. [6] proposed an SEIR model with Holling type III and IV treatment rates for disease dynamics. Motivated by these important works, we have considered the simplified Monod–Haldane (M–H) function type treatment rate in the present epidemic model. In the M–H treatment rate, the treatment rate initially rises with the increase in infectives and achieves the greatest, and after that begins decaying. Such a circumstance may emerge due to the limited availability of treatment for a highly infected population. When the availability of treatment is exhausted (due to the fewer medicines, doctors, and medical staff in hospitals, immunization, etc.), at that point, regardless of the high number of infectives, the accessible medicines turn out to be exceptionally low.

The present article aims to study the transmission and control of infectious diseases when psychological effects and nonmonotonic treatment of infected individuals exist simultaneously in the community. Therefore, to model this situation, we propose an SIR epidemic model with M–H type incidence rate and treatment rate. For the dynamics of the proposed model, we derive the basic reproduction number [5] of the model and study the dynamical behavior through the stability analysis. We give a precise indication of the bifurcation phenomenon using center manifold theory, which ensures the occurrence of the forward bifurcation. We prove the uniform persistence of the system and show the local and global stability of the equilibria under certain conditions. The existence of the limit cycle has also been shown, which shows the occurrence of Hopf bifurcation, regarding the transmission rate as a bifurcation parameter.

The rest of this paper is structured as follows: in Sect. 2, an SIR mathematical model is formulated. In Sect. 3, the basic properties of the model are studied. The local and global stability behaviors of equilibria are discussed in Sect. 4. The model is numerically simulated in Sect. 5. Finally, a brief discussion and conclusion of the work are presented in Sect. 6.

Mathematical Model

In this section, we propose a mathematical disease transmission model with an assumption that the total population size is constant. We divide the total population into three subpopulations or classes, which are named as: susceptible , infected , and recovered . These subpopulations may vary with respect to the time , but ultimately the sum of subpopulations is fixed, i.e., . Further, it is assumed that a disease can spread due to the direct contact between susceptibles and infectives only. It is also assumed that recovered individuals are immune to their entire life, and they will not reinfect in their rest of life. Let be the constant recruitment rate of susceptibles, be the natural death rate of the population, be the transmission rate of a disease, be the rate of psychological effects or behavioral changes due to the susceptibles, be the death rate due to the disease, be the recovery rate of infected individuals, be the cure rate to disease, and be the limitation rate in treatment availability. These assumptions lead to the following nonlinear system of ordinary differential equations:whereIn model (2.1), the term , represents the Monod–Haldane (M–H) functional type incidence rate, which describes a nonmonotonic behavior of incidence rate due to the psychological effect or behavioral changes of susceptibles in case of the high density of infectives in the community. Here, measures the infection force of the disease and describes the inhibition effects of the behavioral changes of susceptible individuals when the number of infected individuals increases. The rate of inhibitory effect may be affected by government interference. The term represents the M–H functional type treatment rate, which describes the nonmonotonic behavior of the treatment rate due to the limitation in the treatment availability of infectives. A summary of the symbols used in the proposed model is given in Table 1.

Table 1

Summary of symbols

Symbol	Description
S	Susceptible population
I	Infected population
R	Recovered population
A	The constant recruitment rate of susceptibles
μ	Natural death rate
β	Transmission rate
α	Rate of inhibitory or psychological effect
d	Disease-induced death rate
δ	Recovery rate
a	Cure rate
b	Limitation rate in treatment availability
R_0	Basic reproduction number
P	Total constant population
QAμ,0	Disease-free equilibrium
Q^∗S^∗,I^∗	Endemic equilibrium
β^∗	Bifurcation parameter
GS,I	M–H type incidence rate
hI	M–H type treatment rate

Basic Properties of the Model

For biological reasons, it is assumed that all state variables and parameters of model (2.1) are nonnegative. We now show that the model solutions always lie in a compact region as follows:

Theorem 3.1

All solutions of model (2.1) starting in with initial conditions (2.2) are bounded and enter the positively invariant region for all time

Proof

Note thatTherefore,It follows that all the solutions of model (2.1) lie in the invariant set

It is noticed that the right-hand side of model (2.1) is completely continuous and locally Lipschitzian on . Thus, it follows that the solution of model (2.1) exists and is unique. It completes the proof of the theorem.□

Under the assumption that the total population size is constant, model (2.1) may be reduced to a two-dimensional system. The third equation (equation of ) of model (2.1) is traditionally omitted; this equation is decoupled from the first two equations of model (2.1). The condition may be used to find . Therefore, the following reduced system is considered for further mathematical analysis:

Equilibria and Their Stability Analysis

System (3.1) has two nonnegative equilibria, obtained by equating the derivatives of the system to zero. These are given as follows:

Disease-free equilibrium (DFE) .

Endemic equilibrium (EE) .

Local Stability Analysis of Disease-Free Equilibrium

To investigate the local stability behavior of disease-free equilibrium, we need to compute the basic reproduction number [7, 5]. Therefore, using the next-generation matrix method [5, 4], the basic reproduction number is computed as follows.

Letwhere , is the matrix of new infection term, and is the matrix of transfer terms into compartments and out of compartments. The Jacobian matrices of and at are given byNow,By using the concept of [5], the basic reproduction number is the spectral radius of the next-generation matrix . Therefore, the basic reproduction number for our model isHence, by the above basic reproduction number and theorem 2 of [5], we propose the following theorem:

Theorem 4.1

The disease-free equilibrium is locally asymptotically stable (LAS) when and unstable when .

Biologically, Theorem 4.1says that the disease will die out from the society whenever the basic reproduction number is less than unity, and a positive equilibrium will exist if .

Analysis of

In this section, we analyze the local stability behavior of disease-free equilibrium when is equal to one. For this, first, we redefine . Then, system (3.1) can be rewritten asNow, the linearized matrix of system (4.1) at and at the bifurcation parameter is given byThe matrix has a simple zero eigenvalue at , and another eigenvalue of the matrix has a negative real part. Thus, the linearization technique fails to determine the behavior of system (4.1) [7, 8]. Therefore, we adopt the center manifold theory [22] to analyze the behavior of disease-free equilibrium at . Using Theorem 4.1 of [2], the bifurcation constants and can be computed asandwhere and are the right and the left eigenvectors of the matrix corresponding to a zero eigenvalue, respectively. Then, we haveThe nonzero partial derivatives of the functions of system (4.1) at and areTherefore,andIt can be seen that is negative and is positive. Hence, bifurcation is forward. Therefore, we propose the following theorem:

Theorem 4.2

System (3.1) exhibits a forward transcritical bifurcation at and bifurcation parameter .

This theorem implies that a positive equilibrium will always exist whenever slightly crosses the unity.

The forward transcritical bifurcation is illustrated in Fig. 1.

Fig. 1

Forward transcritical bifurcation graph in support of theorem (3.1) with numerical data as given in Table 2

Table 2

Parameters and their numerical values

Parameter	P	A	μ	β	α	d	δ	a	b
Numerical value	100	2	0.02	0.003	0.0005	0.02	0.002	0.02	0.005

Existence of Endemic Equilibrium

For the existence of an endemic equilibrium , system (3.1) is rearranged to get and , which gives:where is given by the equationwhereNow, for , using the Descartes’ rule of signs [24], the biquadratic Eq. (4.2) has a unique positive real root if the following condition holds:After getting the value of , we can obtain the value of . Hence, there exists a unique endemic equilibrium if the above condition holds. Hence, we state the following theorem:

Theorem 4.3

If , system (3.1) admits a unique endemic equilibrium .

Now, we show the uniform persistence of system (3.1). Let denote the interior of and denote the boundary of . Epidemiologically, persistence means that the subpopulation always exists and will not lead to elimination if initially, they exist. For this, the following theorem has been proved:

Theorem 4.4

If , system (3.1) is uniformly persistent, which means that there exists a positive constant such that every solution of system (3.1) with the initial data satisfieswhere is independent of initial data in .

Proof

From Theorem 4.3, for , there exists a unique endemic equilibrium . From Theorem 4.1, we know that implies that the DFE is unstable. By theorem 4.3 in [9], the instability of , together with , implies the uniform persistence of the state variables of system (3.1). Therefore, there exists a positive constant such that every solution of system (3.1) with the initial data satisfieswhere is independent of initial data in .□

The uniform persistence, together with the boundedness of , is equivalent to the existence of a compact set in the interior of which is absorbing for system (3.1) [13]. So, we have

Theorem 4.5

If , then there exists a compact absorbing set .

Local Stability Analysis of Endemic Equilibrium

The local stability of the endemic equilibrium is explored as follows.

The variational matrix corresponding to of system (3.1) is given byThe characteristic equation of the variational matrix is given as:wherewhereThe eigenvalues of the matrix have a negative real part if and only if and . Thus, by Routh–Hurwitz criterion, the following results are proposed in the form of theorems given as follows:

Theorem 4.6

Endemic equilibrium is

locally asymptotically stable if the conditions and hold simultaneously.

a saddle point if the conditions and hold simultaneously.

unstable if and .

Biologically, Theorem 4.6 says that when disease persists in society, it can be controlled under condition (i); otherwise, the disease will continue to spread, under conditions (ii) and (iii), respectively.

Theorem 4.7

System (3.1) exhibits a Hopf bifurcation near to the endemic equilibrium if the conditions and hold simultaneously.

The condition implies that in Eq. (10), and the condition implies that . Thus, Eq. (4.3) has purely imaginary roots. From Theorem 4.6(i) and (ii), it follows that the equilibrium changes its behavior from stable to instability as the parameter crosses the critical value , whereAgain, we haveHence, system (3.1) shows a Hopf bifurcation near when .

Theorem 4.7 says that if model (3.1) undergoes a Hopf bifurcation under the conditions mentioned above, then (from the Poincare and Bendixson theorem) every solution initiating the limit cycle must approach a stable limit cycle.□

Next, we derive the condition for the nonexistence of a periodic solution in the interior of the positive quadrant of the S–I plane. For this, we prove the following theorem:

Theorem 4.8

If , then system (3.1) does not have any periodic solution in the interior of the positive quadrant of S–I plane.

We take a real-valued function in the interior of the S–I plane given as follows:Let us consider,Then,It can be seen that the above expression can never be equal to zero when , and also, its sign will be unchanged in the positive quadrant of the S–I plane. Thus, using Dulac’s criterion [22], it can be said that system (3.1) does not have any periodic solution in the interior of the positive quadrant of the S–I plane. Epidemiologically, this theorem says that if the given condition holds, the disease will not reappear in the society.

Since the set defined in Theorem 3.1 is a positively invariant set, the following theorem is a direct result of the Poincare–Bendixson theorem [22] showing the existence of a limit cycle in the interior of the positive quadrant of the S–I plane.□

Theorem 4.9

If either the conditions or are satisfied simultaneously, system (3.1) has at least one limit cycle in the interior of the positive quadrant of the S–I plane.

Epidemiologically, the above theorem implies that if an endemic equilibrium is a saddle point or unstable, then disease may re-occur in society in the future.

Global Stability Analysis

In this section, we discuss the global stability behavior of disease-free and endemic equilibria by Lyapunov functions.

Global Stability of Disease-Free Equilibrium

In this subsection, the global stability behavior of disease-free equilibrium is discussed by Lyapunov function given as follows:

Theorem 4.10

The disease-free equilibrium is globally asymptotically stable when and the condition hold simultaneously.

Proof

Let be the Lyapunov function defined as:where

Differentiating with respect to , we getSince all parameters of the model are positive, it follows that if and hold simultaneously, and if and . Hence, is a Lyapunov function on .

It implies that the largest compact invariant set in is the singleton set . From LaSalle’s invariance principle [22, 18], DFE is globally asymptotically stable.□

Global Stability of Endemic Equilibrium

To prove the global stability of an endemic equilibrium , we consider the following hypothesis:

H(1)

for all

Theorem 4.11

Under hypothesis H(1), endemic equilibrium is globally asymptotically stable.

Let be the Lyapunov function defined as:Differentiating with respect to , we getNote thatTherefore,Hence, by , we see that for all It is easy to verify that the largest invariant in is singleton . By the Lyapunov–LaSalle asymptotic stability theorem [22, 18], is globally asymptotically stable.□

Numerical Simulations

In this section, system (3.1) is simulated numerically using Mathematica 11 (Table 2). The following numerical experimental values of parameters are taken for the simulation:

The initial values of populations are as follows:At the above parameters values, the endemic equilibrium point of system (3.1) is calculated as Therefore, of system (1) can be calculated as .

Figure 2 shows the combined population of susceptible and infected individuals. From this figure, it can be observed that as time passes, the number of susceptibles is decreasing, the number of infectives is increasing, and finally, both populations are approaching the endemic equilibrium with respect to the time.

Fig. 2

Combined population of susceptibles and infectives

Figures 3 and 4 show the infected population at various values of transmission rate and psychological effects , respectively. Figure 3 shows that the number of infected individuals increases with the increase in transmission rate, and Fig. 3 shows that the number of infected individuals decreases with the increase in the psychological or inhibitory effect. In both figures, the infected population is initially increasing, and as time passes, the infected population approaches their steady state.

Fig. 3

Infected population at increased values of transmission rate

Fig. 4

Infected population at increased values of psychological effects

Figure 5 shows the infected population with different values of the initially infected individuals . It can be observed from the figure that the infected population approaches the same steady state at any value of with respect to the time.

Fig. 5

Infected population with different values of the initially infected individuals

Figure 6 shows the infected populations with and without M–H treatment rate. It can be observed from the figure that when the treatment to the infected population is given according to M–H treatment rate, the total number of infected individuals is less, while the total number of infected individuals is high when M–H treatment rate is not considered. Thus, the M–H treatment rate plays a vital role in the control of the spread of infection in society.

Fig. 6

Infected individuals with and without M–H treatment rate

Figure 7 shows the difference in the number of infected individuals when different types of incidence rates, namely Holling type I (H-I), Holling type II (H-II), and Monod–Haldane (M–H), are considered. It shows the effect of incidence rates on the number of infected individuals. It can be seen that, in comparison with the H-I and H-II incidence rates, the number of infectives is less when we consider the M–H type incidence rate for disease dynamics, which has an essential role in disease prevention. Reducing incidence is one of several possible strategies to prevent the spread of infectious disease. Thus, the M–H incidence rate has a vital role in the study of disease dynamics.

Fig. 7

The behavior of the infected population for different incidence rates

To obtain the numerical result in support of Theorem 4.7, we consider the following values of the parameters (Table 3).

Table 3

Parameters and their numerical values

Parameter	P	A	μ	β	α	d	δ	a	b
Numerical value	100	2	0.02	0.043	0.0005	0.02	0.002	2	0.005

Figure 8 shows the occurrence of the limit cycle and confirms the existence of Hopf bifurcation. Clearly, this figure shows the stable limit cycle with an equilibrium point .

Fig. 8

Limit cycle in S–I plane

Discussion and Conclusion

In this study, we proposed and analyzed a mathematical SIR disease transmission model to investigate the role of psychological effects and the limitation in treatment availability. We considered the nonmonotonic Monod–Haldane functional as the incidence rate and the treatment rate. From the mathematical analysis, we obtained that our model has two equilibria, namely disease-free and endemic. Further, we investigated that the disease-free equilibrium is locally asymptotical stable when the basic reproduction number is less than unity. For , we obtained that disease-free equilibrium exhibits a forward transcritical bifurcation which implies that a positive equilibrium will always exist whenever slightly crosses the unity. We checked the existence of the endemic equilibria and obtained that the endemic equilibrium is locally asymptotically stable, a saddle point, and unstable when the conditions in Theorem 4.6 satisfied, respectively. Further, we investigated that system (3.1) exhibits a Hopf bifurcation under the conditions, as stated in Theorem 4.7. We also investigated that system (3.1) is uniformly persistent under the conditions, stated in Theorem 4.4. Moreover, global stability behavior of disease-free and endemic equilibria has been analyzed and obtained that disease-free equilibrium is globally asymptotically stable when and endemic equilibrium is globally asymptotically stable under the hypothesis , when .

From numerical simulations, we observed that the number of infected individuals increases when the transmission rate increases (Fig. 3), while the number of infected individuals decreases when the psychological effect increases (Fig. 4). Further, we observed that the number of infected individuals is low when the treatment rate of the infected individuals is managed according to the M–H treatment rate (Fig. 6). We also portray the various bifurcations through numerical simulation, such as forward transcritical and Hopf bifurcations.

In conclusion, we find that the SIR model with a combination of Monod–Haldane incidence rate and treatment rate is capable of successfully addressing the dynamics of such diseases where the psychological effects and limitation rate in treatment availability are significantly available.