key: cord-0056781-m7ilquyh
authors: Kumar, Abhishek; Nilam
title: Effects of Nonmonotonic Functional Responses on a Disease Transmission Model: Modeling and Simulation
date: 2021-03-02
journal: Commun Math Stat
DOI: 10.1007/s40304-020-00217-4
sha: 20f841fd6e75d698ddc826fc0bf655239df004c0
doc_id: 56781
cord_uid: m7ilquyh

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 [Formula: see text] , we investigate that DFE is locally asymptotically stable when [Formula: see text] and unstable when [Formula: see text] . The local stability of DFE at [Formula: see text] 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.

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 (β S I ), 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 β I l S measures the infection force of the disease, 1 1+α I h describes the inhibition effects from the behavioral changes of susceptibles individuals when the number of infected individuals increases, β, l, and h are positive constants, and α is a nonnegative constant [3, 23, 21] . Note that bilinear incidence (β S I ) is a particular case of (1.1) with l 1 and α 0. The case when l 1 and h 1 in (1.1), i.e., f (I ) β I 1+α I , was proposed by Capasso and Serio [1] to describe the crowding effect in modeling the cholera epidemic in Bari in 1973 . The function f (I ) β I S 1+α I is also known as the Holling type II function. The case when l 1 and h 2 in (1.1), i.e., f (I ) β I 1+α I 2 , represents a nonmonotonic incidence function which can be used to interpret the "psychological effects" [23] . The function f (I ) β I S 1+α I 2 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.

In this section, we propose a mathematical disease transmission model with an assumption that the total population size (P) is constant. We divide the total population into three subpopulations or classes, which are named as: susceptible S(t), infected I (t), and recovered R(t). These subpopulations may vary with respect to the time t, but ultimately the sum of subpopulations is fixed, i.e., S(t) + I (t) + R(t) P constant. 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 A 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, d be the death rate due to the disease, δ be the recovery rate of infected individuals, a be the cure rate to disease, and b be the limitation rate in treatment availability. These assumptions lead to the following nonlinear system of ordinary differential equations:

In model (2.1), the term G(S, I ) β S I 1+α I 2 , 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, β I S measures the infection force of the disease and 1 1+α I 2 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 h (I ) a I 1+bI 2 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 .

For biological reasons, it is assumed that all state variables (S, I , R) and parameters (A, β, α, μ, d, δ, a, b) of model (2.1) are nonnegative. We now show that the model solutions always lie in a compact region D as follows:

Therefore,

It follows that all the solutions of model (2.1) lie in the invariant set D (S, I , R) ∈ R 3 + : S + I + R ≤ A μ . It is noticed that the right-hand side of model (2.1) is completely continuous and locally Lipschitzian on D. 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 P is constant, model (2.1) may be reduced to a two-dimensional system. The third equation (equation of R) of model (2.1) is traditionally omitted; this equation is decoupled from the first two equations of model (2.1). The condition P S(t) + I (t) + R(t) Constant may be used to find R(t). Therefore, the following reduced system is considered for further mathematical analysis:

(3.1)

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

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

is the matrix of new infection term, and V (x) is the matrix of transfer terms into compartments and out of compartments. The Jacobian matrices of Z (x) and V (x) at Q A μ , 0 are given by

By using the concept of [5] , the basic reproduction number is the spectral radius of the next-generation matrix Z V −1 . Therefore, the basic reproduction number (R 0 ) for our model is

Hence, by the above basic reproduction number and theorem 2 of [5] , we propose the following theorem:

Biologically, Theorem 4.1 says 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 R 0 > 1.

In this section, we analyze the local stability behavior of disease-free equilibrium when R 0 is equal to one. For this, first, we redefine S x 1 and I x 2 . Then, system (3.1) can be rewritten as

Now, the linearized matrix of system (4.1) at Q A μ , 0 and at the bifurcation parameter β β * μ(μ+d+δ+a) A is given by

The matrix J has a simple zero eigenvalue at R 0 1, 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 R 0 1. Using Theorem 4.1 of [2] , the bifurcation constants a 1 and b 1 can be computed as

where w (w 1 , w 2 ) T and u (u 1 , u 2 ) are the right and the left eigenvectors of the matrix J corresponding to a zero eigenvalue, respectively. Then, we have

The nonzero partial derivatives of the functions f i s (i 1, 2) of system (4.1) at R 0 1 and β β * are

Therefore,

It can be seen that a 1 is negative and b 1 is positive. Hence, bifurcation is forward. Therefore, we propose the following theorem: This theorem implies that a positive equilibrium will always exist whenever R 0 slightly crosses the unity.

The forward transcritical bifurcation is illustrated in Fig. 1 .

For the existence of an endemic equilibrium Q * (S * , I * ), system (3.1) is rearranged to get S * and I * , which gives:

where I * is given by the equation

where C 1 bμα(μ + d + δ),

Now, for R 0 > 1, using the Descartes' rule of signs [24] , the biquadratic Eq. (4.2) has a unique positive real root I * if the following condition holds:

After getting the value of I * , we can obtain the value of S * . Hence, there exists a unique endemic equilibrium Q * (S * , I * ) if the above condition holds. Hence, we state the following theorem: and ∂ D r denote the boundary of D r . 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: Table 2 123

where K is independent of initial data in D 1 .

Proof From Theorem 4.3, for R 0 > 1, there exists a unique endemic equilibrium Q * . From Theorem 4.1, we know that R 0 > 1 implies that the DFE (Q) is unstable. By theorem 4.3 in [9] , the instability of Q, together with Q ∈ ∂ D r , implies the uniform persistence of the state variables of system (3.1). Therefore, there exists a positive constant K such that every solution (S, I ) of system (3.1) with the initial data (S(0),

where K is independent of initial data in D 1 .

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

The local stability of the endemic equilibrium Q * (S * , I * ) is explored as follows.

The variational matrix J Q * corresponding to Q * (S * , I * ) of system (3.1) is given by

The characteristic equation of the variational matrix J Q * is given as:

where

where la a 1 + I * 2 α 2 + 1 + bI * 2 2 d 1 + I * 2 α 2 + δ + 2μ + I * (1 + I * (S * + I * )α)β + I * α 2 + I * 2 α (δ + 2μ) 1 + bI * 2 2 1 + I * 2 α 2 , lb ab I * + I * 3 α 2 + S * 1 + bI * 2 2 β ma a 1 + I * 2 α (μ + I * (β + I * αμ)) + 1 + bI * 2 2 I * 1 + I * 2 α β(d + δ) + d 1 + I * 2 α 2 + I * (1 + I * (S * + I * )α)β + 1 + I * 2 α 2 δ μ + μ + I * 2 αμ 2 1 + bI * 2 2 1 + I * 2 α 2 , mb S * 1 + bI * 2 2 βμ + abI * 2 1 + I * 2 α I * β + μ + I * 2 αμ 1 + bI * 2 2 1 + I * 2 α 2 .

The eigenvalues of the matrix J Q * have a negative real part if and only if p 1 > 0 and p 2 > 0. Thus, by Routh-Hurwitz criterion, the following results are proposed in the form of theorems given as follows: 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. Proof The condition l a l b implies that p 1 0 in Eq. (10) , and the condition m a > m b implies that p 2 > 0. Thus, Eq. (4.3) has purely imaginary roots. From Theorem 4.6(i) and (ii), it follows that the equilibrium Q * (S * , I * ) changes its behavior from stable to instability as the parameter β crosses the critical value β β * , where

0.

Hence, system (3.1) shows a Hopf bifurcation near Q * (S * , I * ) 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:

Proof We take a real-valued function in the interior of the S-I plane given as follows:

Let us consider,

It can be seen that the above expression can never be equal to zero when α ≥ b, 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 D 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. 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.

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

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

Theorem 4.10 The disease-free equilibrium Q A μ , 0 is globally asymptotically stable when R 0 ≤ 1 and the condition ab ≤ α(μ + d + δ) hold simultaneously.

Proof Let L 1 be the Lyapunov function defined as:

Since all parameters of the model are positive, it follows that dL 1 dt < 0 if R 0 ≤ 1 and [22, 18] , DFE is globally asymptotically stable.

To prove the global stability of an endemic equilibrium Q * (S * , I * ), we consider the following hypothesis: Differentiating L 2 with respect to t, we get

Note that

Hence, by H (1), we see that dL 2 dt ≤ 0 for all S, I ≥ 0. It is easy to verify that the largest invariant in (S, I ) ∈ R 2 + : dL 2 dt 0 is singleton {Q * }. By the Lyapunov-LaSalle asymptotic stability theorem [22, 18] , Q * (S * , I * ) is globally asymptotically stable.

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 Q * (S * , I * ) (25.4721, 26.2010) . Therefore, R * of system (1) can be calculated as R * P − S * − I * 48.3269. 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 Q * with respect to the time. 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 123 A. Kumar, Nilam psychological or inhibitory effect. In both figures, the infected population is initially increasing, and as time passes, the infected population approaches their steady state. Figure 5 shows the infected population with different values of the initially infected individuals I (0). It can be observed from the figure that the infected population approaches the same steady state at any value of I (0) with respect to the time. 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 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.

To obtain the numerical result in support of Theorem 4.7, we consider the following values of the parameters (Table 3) . 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 E * (47.87, 1.29).

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 (R 0 ) is less than unity. For R 0 1, we Fig. 7 The behavior of the infected population for different incidence rates obtained that disease-free equilibrium exhibits a forward transcritical bifurcation which implies that a positive equilibrium will always exist whenever R 0 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 R 0 ≤ 1 and endemic equilibrium is globally asymptotically stable under the hypothesis H (1), when R 0 > 1.

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.

A generalization of the Kermack-McKendrick deterministic epidemic model

Dynamical models of tuberculosis and their applications

Complex dynamics of an SIR epidemic model with nonlinear saturate incidence and recovery rate

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

Reproduction numbers and sub-threshold endemic equilibria for compartment models of disease transmission

Modelling and analysis of a SEIR model with different types of nonlinear treatment rates

Dynamics of n SIR model with nonlinear incidence and treatment rate

An SIR model with nonlinear incidence rate and Holling type III treatment rate

Uniform persistence and flows near a closed positive invariant set

Nilam: A mathematical and numerical study of a SIR epidemic model with time delay, nonlinear incidence and treatment rates

Nilam: A deterministic time-delayed SVIRS epidemic model with incidences and saturated treatment

An SVEIR model for assessing the potential impact of an imperfect anti-SARS vaccine

Permanence and the dynamics of biological systems

A contribution to the mathematical theory of epidemics

Nilam: Mathematical analysis of a delayed epidemic model with nonlinear incidence and treatment rates

A short study of an SIR model with inclusion of an alert class, two explicit nonlinear incidence rates and saturated treatment rate

Nilam: A deterministic time-delayed SIR epidemic model: mathematical modeling and analysis

The Stability of Dynamical Systems

Dynamic behavior of a modified SIR model in epidemic diseases using nonlinear incidence and recovery rates

Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological model

Dynamical behavior of an epidemic model with a nonlinear incidence rate

Analysis, Stability and Control

Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate

A simple proof of Descartes's rule of signs

Bifurcation in an epidemic model with constant removal rates of the infectives

Stability of a delayed SIRS epidemic model with a nonlinear incidence rate

Qualitative analysis of a SIR epidemic model with saturated treatment rate

Dynamics of a SIR epidemic model with limited medical resources revisited

The authors are gratefully acknowledged to Delhi Technological University for providing monetary support for this research. The authors are also thankful to the handling editor and anonymous reviewers for their constructive comments and suggestions which enhance the quality of the paper.

Conflict of interest The authors declare that they have no conflict of interest.