key: cord-0886033-030w3ddr
authors: Jain, Shikha; Kumar, Sachin
title: Dynamical analysis of SEIS model with nonlinear innate immunity and saturated treatment
date: 2021-09-17
journal: Eur Phys J Plus
DOI: 10.1140/epjp/s13360-021-01944-5
sha: 52939a362086b092ce793c442f682320a265e2c1
doc_id: 886033
cord_uid: 030w3ddr

In this paper, we develop an SEIS model with Holling type II function representing the innate immunity as well as the saturated treatment. We obtain the existence and stability criteria for the equilibrium points. We observe that when the reproduction number is less than unity, the disease-free equilibrium always exists and is locally asymptotically stable. The multiple endemic equilibrium points can exist independent of the basic reproduction number, and the system may experience bistability. We find that the system can encounter backward or forward bifurcation at [Formula: see text] , where the contact rate [Formula: see text] is the bifurcation parameter. Therefore, the disease-free equilibrium may not be globally stable. We deduce the criteria for the presence of Hopf bifurcation where the parameter [Formula: see text] acts as the bifurcation parameter and the system is a neutrally stable center. We also observe with the aid of a numerical example that a slight perturbation disrupts the neutral stability and the trajectories become either converging or diverging from the equilibrium point. Numerical simulation is performed with the help of MATLAB to justify the findings. We study the effect of nonlinearity of immunity function and the treatment rate on the dynamics of the disease spread. We find that when both are linear, the reproduction number is the same, but the system has a unique endemic equilibrium point that exists for reproduction number greater than unity. We find that there is neither backward bifurcation nor Hopf bifurcation. We also observe that the saturation in treatment enlarges the domain of backward bifurcation making disease eradication an extremely difficult task. The endemic equilibria in the case of saturated treatment may exist far more to the left of the bifurcation parameter [Formula: see text] . Hence, the nonlinearity of immunity function and treatment function affects the dynamics of an SEIS model highly; therefore, one must be precautious to choose an appropriate function for both while modeling.

Mathematical modeling is the art of transforming real-world problems into a set of mathematical equations whose analysis provides insight and useful information regarding the concerned problem which further helps in simplifying it. It has a significant role in the invesa e-mail: shikhajain01051990@gmail.com (corresponding author) b e-mail: sachinambariya@gmail.com tigation of troubles that influence our day-to-day life. Mathematical models have proved to be a remarkable tool in investigating the dynamics of communicable diseases. Moreover, these models have a pivotal part in policymaking for the prevention and control of an infection spread. The first deemed mathematical model for the study of an infectious disease was given in 1760 by Daniel Bernoulli on inoculation against smallpox [1] . Later, communicable disease modeling became one of the prominent domains of research. The compartmental models are considered the building blocks of epidemiological modeling. The first recognized compartmental model was given by Kermack and McKendrick [2] in 1927 which studied the relationship among susceptible, infected, and recovered individuals. In the past, various researchers have developed different compartmental systems (SEIS, SIR, SEIR, SEITR, SIS, etc.) to analyze the spread of communicable diseases.

Treatment is a term coined for the care given to a patient for an illness or an injury. It is a significant tool in controlling the spread of disease. In early epidemic models, the treatment rate was considered to be linear in infected fraction while in reality, it is highly dependent on access to medical resources like vaccines, medicines, hospital cots, quarantine areas, and treatment capability. Each region has restricted medical assets; hence, it is essential to contemplate a suitable treatment function. Treatment was first introduced in epidemiological modeling by Wang and Ruan [3] . They considered the treatment rate to be constant. Later, Wang [4] modified it as a function where treatment rate is proportionate to the infectious population till the treatment capability reaches its maximum limit. Nowadays, a large domain of researchers is dealing with the infectious disease models with a variety of treatment functions [5] [6] [7] [8] [9] [10] [11] .

Innate immunity is the term given to the non-specific defense mechanisms that come into play immediately or within hours of a pathogen intrusion in the body. The different outcomes among individuals after infection as well as the cause why not each individual in contact with an infectious individual catches infection signifies the importance of innate immunity. Innate immunity includes phagocytosis, production of cytokine, immune identification, and effector processes. Innate immunity is a neglected topic among researchers, and some consider it important while mostly it is ignored.

Shi et al. [12] constructed a sepsis disease model which considers innate as well as acquired immunity. Kumar et al. [13] studied a model for inflammation and presented some suggestions in form of therapeutic strategies resulting from the sensitivity analysis. Heffernan et al. [14] studied the measles infection with a model formulated with the combination of the immunological and epidemiological models. Their detailed study included the relations among disease occurrence, dwindling immunity, and boosting. Mwalili et al. [15] proposed an SEIR epidemic system for COVID-19 infection where innate immunity plays an important role in recuperating from incubating stage to the susceptible compartment. Pigozzo et al. [16] developed the algorithmic system for the study of innate immunity. Banerjee et al. [17] constructed a model by combining the immunological and epidemic models to study the two-stage epidemic spread where an increase in the population of the infected class increases the viral load and inhibits innate immunity. For further literature on the relevance of innate immunity, refer [18] [19] [20] [21] [22] [23] . Jain et al. [24] studied the effect of natural immunity in the SEIS epidemic system which is the motivation of the present manuscript. In this paper, we study the effect of innate immunity when the treatment resources are limited. We consider the immunity and saturated treatment both to be the functions of Holling type II and study the effect of both on the dynamics of an SEIS model. We further compare the disease dynamics to the case when the treatment rate is linear. We observe that the saturation in treatment increases the domain of backward bifurcation, that is, the domain of existence of endemic equilibria is quite large when the basic reproduction number is less than unity. This makes disease eradication a difficult task to attain.

This manuscript is arranged as: In Sect. 1, we formulate an SEIS model having three compartments where the individuals become susceptible after recovering from the exposed compartment with the aid of innate immunity, and the treatment resources are limited. In the next section, we prove the positivity and boundedness of the solutions and calculate the reproduction number. In Sect. 3, we perform the mathematical analysis of the model. We find the conditions for the existence and stability of the equilibrium points and study the possible bifurcations. In the next Section, we execute simulations to justify the mathematical results numerically. Section 5 discusses the effect of nonlinearity of immunity function and saturation of treatment on the different classes. In this section, the results of the previous sections are compared with the cases when treatment and immunity functions are linear. The final section sums up the manuscript in the form of a conclusion.

We formulate the mathematical structure of a generic SEIS epidemiological model. We study the impact of innate immunity and limited medical resources, to make the SEIS model more realistic. The Holling type II functional response is used to represent both. The following assumptions are considered to construct the model:

• The inter-individual interaction is homogeneous, and the population is well mixed. • The total population N (t) is classified into three exclusively disjoint but interacting classes: susceptible, exposed, and infected where S(t) , E(t), and I (t) denote the population of each class, respectively. The total population N (t) is

• We consider the constant recruitment rate for the population growth which is denoted by . The natural death rate is represented by δ, while μ denotes the disease-induced death rate.

• We further assume that the disease has no vertical transmission, that is, the offspring of an infected individual belongs to the susceptible class. • There is no direct transmission from the susceptible class to the infected class, that is, one has to undergo the incubation phase to get infected, and σ represents the rate of incubation. • The rate of contact at which infected individuals transmit the disease to susceptible is represented by β. • The innate immunity and treatment functions are Holling type II functions. The Holling type II functions (H (x)) are continuously differentiable functions as well as monotonically increasing and have following properties:

• The incubating or latent beings move into the class of susceptible with the aid of innate immunity at the rate g(E), assumed to be Holling type II function, defined as the following:

where a/k is the maximum immunity efficiency of the community. The term a E denotes the recovery with the aid of innate immunity, while 1/(1 + k E) denotes the hindrance in recovery from exposed to susceptible compartment owing to weakening of immunity due to diabetes, malnutrition, poor sanitary conditions, and HIV infection, etc.

• We suppose that nobody can recover from infected class without proper treatment, that is, to recover from I class one has to take some sort of medical help. It is assumed that due to the increasing number of infections the treatment rate will eventually saturate rather than grow without limits as a linear rate since every region has limited medical resources rather than unlimited. The infected individuals recover from infected to susceptible class with the recovery rate h(I ) which is assumed to be a Holling type II function, defined as follows:

is the maximal treatment efficiency of the community.

The above assumptions result into the following mathematical model:

All the parameters are non-negative, and the initial conditions corresponding to the model (1.3) are

In the next section, we present the positivity and boundedness of the solutions of the system (1.3).

In this section, we show that the solutions of the system (1.3) are bounded and positive. In the following theorem, we establish the region of attraction for the system (1.3). 

The proof of the above theorem is trivial, hence omitted. The above theorem shows that the system (1.3) is biologically well behaved. Next, we study the elementary characteristics of the system (1.3) like basic reproduction number, the existence of equilibria, stability, and bifurcation analysis.

Every epidemic model has two types of equilibria, disease-free equilibrium (DFE) and endemic equilibrium (EE). The DFE of the system (1.3) is denoted by E 0 and E 0 = ( /δ, 0, 0). Next, we compute the basic reproduction number which signifies the total number of expected secondary infections caused by an infection carrier in a purely susceptible population. We calculate the basic reproduction number with the help of the next-generation matrix method [25] [26] [27] . The transmission matrix T , describing the production of new infections, and the transition matrix describing changes in the state, are as follows:

The basic reproduction number is the dominant eigenvalue of the matrix −T −1 and is computed as follows:

The reproduction number R 0 is independent of the parameters k and b.

The Jacobian matrix of the system (1.3) is as follows:

Next, we study the nature of the roots of the Jacobian matrix with the help of the next lemma.

For > 0, the character of eigenvalues of the Jacobian matrix (3.2) , associated with the EE, is described as

, the matrix has unique real eigenvalue and a complex conjugate pair of eigenvalues where the real part is negative.

3 , the matrix has unique real eigenvalue and a pair of purely imaginary eigenvalues.

3 , the matrix has unique real eigenvalue and a complex conjugate pair of eigenvalues where the real part is positive.

Proof The proof is similar to [24, Lemma 2.2], hence, omitted.

Next, we study the existence of endemic equilibria.

To deduce the existence criterion of the endemic equilibrium points, we solve the third equation of model (1.3) and obtain the value of E at the equilibrium point which is as follows:

Next, we substitute the above value of E * in the second equation of system (1.3) to obtain the value of S * as follows

(3.6)

Substituting the above value in the first equation and on simplification, we get

We can observe that the coefficient A 4 > 0 and nothing can be directly deduced about the sign of coefficients A 3 , A 2 and A 1 . But the coefficient A 0 < 0 when R 0 > 1. Hence, there is at least one sign change among the coefficients of (3.7). Therefore, by Descartes' Rule of Signs [28], equation (3.7) will have at least one positive real root when R 0 > 1. Moreover, uncertainty in the signs of coefficients A 3 , A 2 and A 1 suggests a possibility of the existence of endemic equilibrium points when R 0 < 1 as well as the existence of multiple EE's. We will study this later in detail.

The Jacobian matrix of the linearized system (1.3) on evaluation at DFE gives:

The characteristic equation of the matrix J(E 0 ) is given by

where C 2 = δ,

(3.11)

. Clearly, C 2 and C 1 are positive, and C 0 > 0 when R 0 < 1. Therefore, with the help of Ruth-Hurwitz criterion, the DFE is locally asymptotically stable and the next result follows. 

Next, the Variational matrix is evaluated at the DFE and β = β 0 represented by J(E 0 )| β 0 . The matrix J(E 0 )| β 0 has a zero eigenvalue which is simple, and the rest of the eigenvalues are −δ and −(a + 2δ + γ + σ + μ) which are negative. Hence, center manifold theory can be applied here. The right eigenvector of J(E 0 )| β 0 corresponding to the simple zero eigenvalue represented by r = (r 1 , r 2 , r 3 ) T is as follows: Similarly, the left eigenvector l = (l 1 , l 2 , l 3 ) is calculated as

(3.13)

The local stability in the neighborhood of the bifurcation point β 0 = β is decided by the signs of the two associated constants, denoted by α 1 and α 2 , defined as follows:

with φ = β − β 0 . In order to proceed further, we calculate the partial derivatives at DFE and the required nonzero partial derivatives are

Next, we calculate the values of constants α 1 and α 2 , and we obtain

.

(3. 16) We notice that α 2 > 0; therefore, the sign of α 1 is the determining factor for the direction of bifurcation [29, Theorem 4.1]. When α 1 > 0, the system undergoes backward bifurcation, and α 1 < 0 leads to forward bifurcation. This shows that the endemic equilibrium points can exist when R 0 < 1 as we have already observed in (3.8) . Moreover, we perceive that the endemic equilibrium points can exist when R 0 < 1, the DFE may not be globally stable. The next result concludes these observations.

For α 1 < 0, there exists a unique endemic equilibrium E * 1 whenever R 0 > 1, and as the sign of R 0 − 1 changes from negative to positive, E * 1 switches its stability from unstable to locally asymptotically stable and the forward bifurcation occurs at R 0 = 1, and β 0 works as the bifurcation parameter. When R 0 < 1, endemic equilibrium points exist for α 1 > 0, leading to backward bifurcation with β 0 as the bifurcation parameter and α 1 is described as the following

We also observe that nonlinearity of innate immunity and recovery rate, that is, parameters k and b play crucial roles in the occurrence of backward bifurcation. 

We claim the presence of Hopf bifurcation by establishing the existence of limit cycle surrounding the EE E * . We follow the following steps to reach this goal [24, 30] (1) Determine B 1 , Y 1 and Y 2 in the terms of the bifurcation parameter γ .

(2) Simplify the equation Y 1 + Y 2 = −2B 1 /3 and calculate the critical value γ = γ * of the bifurcation parameter. (3) Determine the endemic equilibria related to the particular set of parameter values. (4) Next, the real part of complex eigenvalues is differentiated with respect to γ . We further evaluate it at γ * to verify the transversality condition which is

(3.17)

The remaining eigenvalue λ 1 should be negative at γ = γ * . Next, we establish the presence of Hopf bifurcation with the help of Hopf Bifurcation Theorem [31] . (3.4) .

Proof The characteristic equation of the Variational matrix (3.2) is

, the characteristic equation at the bifurcation parameter γ = γ * reduces to

where B i s = 1, 2 are given by (3.4) . The eigenvalues of the (3.19) are −B 1 and ±ι √ B 2 . For some ε > 0 and γ ∈ (γ * − ε, γ * + ε), the eigenvalues are (3.20) where j i s, i = 1, 2 ∈ R. Now, we verify the transversality condition which is as follows

Next, we substitute λ 2 (γ ) from (3.20) in the characteristic equation (3.18) of the Jacobian matrix (3.2), and we get

Differentiating with respect to γ , we get 3( j 1 (γ ) + ιj 2 (γ )) 2 (j 1 (γ ) + ιj 2 (γ )) + 2B 1 ( j 1 (γ ) + ιj 2 (γ ))(j 1 (γ ) + ιj 2 (γ )) +Ḃ 1 ( j 1 (γ ) + ιj 2 (γ )) 2 + B 2 (j 1 (γ ) + ιj 2 (γ )) +Ḃ 2 ( j 1 (γ ) + ιj 2 (γ )) +Ḃ 3 = 0.

(3.23)

Comparing the real and imaginary parts, we get

where

Now, solving the equations in (3.23) forj 1 , we geṫ

(3.26)

Now, at γ = γ * , Case 1. When

Therefore,

and λ 1 (γ * ) = −B 1 (γ * ) < 0. Hence, the Hopf Bifurcation Theorem [31] proves the existence of limit cycle. Therefore, Hopf bifurcation occurs at the bifurcation parameter γ = γ * .

For a = 0.7, β = 0.2, b = 0.8, = 0.9, , δ = 0.01, μ = 0.02, k = 0.01, σ = 0.03 and γ = 1.35052, the Hopf bifurcation occurs, and γ = γ * is the bifurcation parameter. Next, the transversality condition results

Therefore, Hopf bifurcation takes place at γ = γ * around EE E 2 . Moreover, we note the following observations while varying the value of γ *

(1) When γ = γ * = 1.35052, the system has three equilibrium points, E 0 = (90, 0, 0), E 1 = (40.3553, 41.4098, 2.74497) and E 2 = (26.2739, 49.3865, 4 .77987). The system undergoes Hopf bifurcation around the EE E 2 , and γ * is the critical value of bifurcation parameter. The system is a neutrally stable center with concentric limit cycles around the EE E 2 . Figure 1a and b show the presence of concentric limit cycles around E 2 . The different colors represent the trajectories obtained with different initial conditions. (2) When γ = 1.345 < γ * , the system is no more a neutrally stable center, and the EE becomes locally asymptotically stable as well as the limit cycles convert into attracting spirals toward E 2 . Figure 1c shows the existence of attracting spirals. (3) When γ = 1.352 > γ * , the EE becomes unstable and the trajectories are diverging away from E 2 in the form of unstable spirals. Figure 1d shows the diverging trajectories from the EE E 2 .

In this section, we justify the theoretical results with the help of some numerical examples. We use MATLAB for carrying out the simulations with the aid of several sets of parameters which are demonstrated in the next examples: Next, we discuss the significance of the nonlinearity of the immunity function and the treatment rate as well as analyze the effect of variation of parameters k and b on the different population classes.

• Effect of variation in b: Fig. 3 represents the effect of variation in the value of parameter b on the distinct population compartments. It can be observed that when b = 0, the dynamics completely change as the curves in Fig. 3 for b = 0 are quite distinguishable. Hence, it is very easy to observe that the saturated treatment changes the dynamics of disease spread to a high extent. • Effect of variation in k: Similarly, Fig. 4 represents the effect of variation in the value of k on the distinct population compartments. It can be observed that in the neighborhood of k = 0.5, a slight variation may cause convergence to another equilibrium point.

Hence, we deduce that though the basic reproduction number is independent of b and k, these parameters still have a high influence on the spread of disease. 

We have already observed that the parameters b and k do not effect the basic reproduction number. However, (3.16) shows that there is no backward bifurcation when both of them are zero. When both the parameters k and b are zero, the system (1.3) reduces to the following:

The basic reproduction number for the system (5.1) is same as for system (1.3) and given by (3.1). The system (5.1) has unique DFE, E 0 = ( /δ, 0, 0) and unique endemic,

The endemic equilibrium exists only when R 0 > 1 and is unique. We further observe that the nonlinearity of the immunity function and treatment rate forbids the possibility of multiple endemic equilibria as well as backward bifurcation. Hence, the DFE is globally stable while for model (1.3), the DFE may not be globally stable. Hence, the nonlinearity of immunity and treatment function changes the dynamics of infection spread highly distinct from when they are linear. We have proved the existence of Hopf bifurcation in the model (1.3), and next, we show that when parameters b and k are zero, the limit cycle does not exist. We define the new variables

Using (5.3), the system (5.1) transforms to the following: 

(5.7)

Suppose g = g 1 + g 2 + g 3 where

We observe thatf · g = 0 in D, wheref = (f 1 ,f 2 ,f 3 ). We further calculate

where n = (1, 1, 1 + μ 1 ) is the normal vector. Therefore, by [32, Corollary 4.2] , there is no periodic solution of (5.4) and ultimately no periodic solution of system (5.1). Hence, there is no Hopf-bifurcation and the unique endemic equilibrium E * 1 is globally stable. We conclude that Hopf bifurcation cannot exist whenever treatment and immunity functions are linear. We have observed throughout this section how much nonlinearity of immunity function and saturation in treatment affect the dynamics of disease spread.

Next, we compare the results obtained in this manuscript with the work done by Jain et al. [24] . We observe that the system (1.3) reduces to the model given by [24] when we consider b = 0. The constant which determines the direction of transcritical bifurcation at R 0 = 1, α 1 for the system (1.3) is given by (3.16) and for b = 0, as calculated in [24] is α 1 = 2ak(δ + γ + μ) 2 σ (a + δ + σ ) − β 0 (δ(γ + δ + σ ) + μ(δ + σ )) δ(a + δ + σ ) . (5.10) We can see that they differ by the term 2bγ . This implies that when the treatment is saturated or the medical facilities are limited, the domain of backward bifurcation is enlarged. Hence, the endemic equilibria in case of saturated treatment may exist far more to the left of bifurcation parameter β 0 or simply the domain of existence of endemic equilibria is quite large when R 0 < 1. This makes disease eradication quite a difficult task. Hence, the saturated treatment makes disease eradication a more difficult task to attain.

We have observed that the saturation in treatment and innate immunity both change the dynamics of disease spread in an SEIS model in comparison to the model where both are considered linear. We find that linearity of both results in a system where a unique endemic exists, that too for R 0 > 1. There is neither backward bifurcation nor Hopf bifurcation.

As there cannot exist multiple endemic equilibria, bistability cannot occur. Moreover, the saturation in treatment enlarges the domain for backward bifurcation and makes disease eradication a difficult task. Hence, we conclude that one must consider a proper form of immunity and treatment function while modeling as a wrong choice can deviate the system from a real scenario to a great extent. This study can be quite advantageous for modeling the diseases such as TB, dengue, flu, Zika, SARS (Severe acute respiratory syndrome), and the current COVID-19 pandemic. The application of natural immunity with limited treatment resources can be further broadened to distinct compartment models such as SEIR, SVEIS, SEITRS, SEITS. In this paper, we have formulated a simple system for the study of innate immunity in an environment with restricted medical assets. This study can be further diversified with various constraints such as nonlinear incidence, herd immunity and hence, can be very beneficial in epidemiological modeling.

The mathematics of infectious diseases

A contribution to the mathematical theory of epidemics

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

Backward bifurcation of an epidemic model with treatment

Dynamics of an epidemic model with quadratic treatment

Chaos detection in SIR model with modified Beddington-DeAngelis type incidence rate and saturated treatment

Assessing the effects of treatment in HIV-TB co-infection model

Progressive dynamics of a stochastic epidemic model with logistic growth and saturated treatment

Backward bifurcation of an epidemic model with saturated treatment function

Dynamics of an SIR epidemic model with limited medical resources revisited

Bifurcation analysis of an SIS epidemic model with saturated incidence rate and saturated treatment function

Mathematical model of innate and adaptive immunity of sepsis: a Modeling and simulation study of infectious dDisease

The dynamics of acute inflammation

Implications of vaccination and waning immunity

SEIR model for COVID-19 dynamics incorporating the environment and social distancing

On the computational modeling of the innate immune system

Immuno-epidemiological model of two-stage epidemic growth

A reduced mathematical model of the acute inflammatory response II. Capturing scenarios of repeated endotoxin administration

Cellular innate immunity: an old game with new players

Mathematical modelling of immune response in tissues

silico models of acute inflammation in animals. Shock

Global stability of a virus dynamics model with Beddington-De Angelis incidence rate and CTL immune response

Analyzing global stability of a viral model with general incidence rate and cytotoxic T lymphocytes immune response

Dynamic analysis of the role of innate immunity in SEIS epidemic model

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

The construction of next-generation matrices for compartmental epidemic models

A simple proof of Descartes's rule of signs

Dynamical models of tuberculosis and their applications

Chaos detection and optimal control in a cannibalistic prey-predator system with harvesting

Mathematical biology

Analysis of a disease transmission model in a population with varying size

The authors would like to thank the anonymous referees for their extensive comments on the revision of the manuscript which improved the quality of the paper. The author, Shikha Jain, is financially supported by the University Grant Commission (UGC), Government of India (Sr. No. 2061440971). She gratefully acknowledges the support for the research work.

 Page 15 of 17 952 Consider D = (N ,Ē,Ī ) ∈¯ :N = 1 − μ 1Ī = (S,Ē,Ī ) ∈¯ :S +Ē + (1 + μ 1 )Ī = 1 . We observe that dN /dτ = 0 in D, and D is an attractive invariant set in¯ as all the solutions of (5.4) in¯ tend to D monotonically when τ tends to infinity. Hence, all ω-limit sets of (5.4) are contained in D. Therefore, it is sufficient to show that there is no periodic solution in D.Suppose f i , i = 1, 2, 3 denote the left hand sides of (5.4) maintaining the sequence, andf i denote the restriction of f i on D. Next, we computef k (Ē,Ī ),f j (S,Ī ), andf l (S,Ē) using the relationS +Ē + (1 + μ 1 )Ī = 1 where k = 2, 3; l = 1, 2 and j = 1, 3. We obtain the followingf