key: cord-0713779-3lzextoc
authors: Ajbar, Abdelhamid; Alqahtani, Rubayyi T.
title: Bifurcation analysis of a SEIR epidemic system with governmental action and individual reaction
date: 2020-10-01
journal: Adv Differ Equ
DOI: 10.1186/s13662-020-02997-z
sha: aefe7da2925c3f241653f40dce8b6ba3afd7e3ff
doc_id: 713779
cord_uid: 3lzextoc

In this paper, the dynamical behavior of a SEIR epidemic system that takes into account governmental action and individual reaction is investigated. The transmission rate takes into account the impact of governmental action modeled as a step function while the decreasing contacts among individuals responding to the severity of the pandemic is modeled as a decreasing exponential function. We show that the proposed model is capable of predicting Hopf bifurcation points for a wide range of physically realistic parameters for the COVID-19 disease. In this regard, the model predicts periodic behavior that emanates from one Hopf point. The model also predicts stable oscillations connecting two Hopf points. The effect of the different model parameters on the existence of such periodic behavior is numerically investigated. Useful diagrams are constructed that delineate the range of periodic behavior predicted by the model.

The spread of infectious diseases is a very complex phenomenon that depends on a large number of factors. Some of them are social, environmental, or economic, which are linked to human activities while other factors have to do with the nature of the pathogen causing the disease. With this complexity, it is not surprising that sometimes simple mathematical models are often used to understand the dynamics of spread of infectious diseases. The use of more complex models will necessarily involve very large number of parameters that are difficult to estimate, making the model predictions weak and uncertain. In this regard, compartmental models are often used to simplify the mathematical modeling of infectious diseases. In these models, the population is divided into compartments, with the assumption that every individual in the same compartment has the same characteristics [1] . The models are usually investigated through deterministic ordinary differential equations. The SEIR mathematical model is an extensively used compartmental epidemic model that is based on the division of the population into four basic compartments; an individual can either be susceptible (S), exposed to the disease but not yet infectious (E), infectious (I), or removed (recovered or deceased) (R). The SEIR model has been used extensively in the literature to model many human infections such as Ebola [2] , H1N1 [3] , influenza [4] and MERS-CoV [5] . Despite being an old model, the SEIR model is quite flexible and many variations and improvements on the model are still possible. Right now, the SEIR model has been applied extensively to analyze the COVID-19 pandemic [6] [7] [8] [9] . Each of these studies includes a variation on the basic SEIR model by either taking into consideration new variables or parameters, ignoring others, selecting different expressions for the transmission rate, or using different methods for parameter estimation. The stability, bifurcation and chaotic behavior of SEIR epidemic models have been investigated over many decades [10] [11] [12] [13] [14] [15] [16] [17] [18] . Tools from nonlinear theory were shown to be useful in revealing conditions for the occurrence of sustained epidemic equilibria. Some early work [10, 11] analyzed the effect of seasonal fluctuations as well as contact rate periodicity in what becomes a forced response problem resulting in harmonic and subharmonic resonances.

Several authors have analyzed the occurrence of periodic solutions in SEIR models due to the presence of time delays and/or nonlinear incidence rates [12] [13] [14] [15] [16] . Chaotic motion has also been studied in such models [17, 18] .

In this paper we examine the bifurcation behavior of a SEIR model when the transmission rate takes explicitly into account governmental action and population response to the severity of the pandemic. Although a number of studies have included the effects of governmental action and individual response in models of COVID-19 [6, 19] , there is few work in the literature on the bifurcation and stability of such models. The recent work of Kwuimy et al. [19] showed the importance of governmental action and social behavior in COVID-19 dynamics.

A final note is to be made about the usefulness of such models in the forecast of the spread of COVID-19 pandemic. Because of the absence of complete data for the disease in Saudi Arabia, the validation of the model could not be carried out. However, the values of model parameters that were used in the mathematical analysis were taken from the literature where a validation of a similar model was carried out for the Wuhan province in China [6] . This provides the proposed model with some credibility. Also, one of the objectives of this paper is to carry out a sensitivity analysis for the effect of model parameters. In this regard, the numerical analysis was carried out for a wide range of values of model parameters, which should also give some credibility to the different behavior predicted by the model. Moreover, it is true that some of the behavior predicted by the proposed model (such as periodic solutions), were not reported so far in the current literature on the COVID-19 disease. However, it is also true that the pandemic is unfortunately not yet contained, with the unfortunate possibility of a second wave or more. Given that the proposed model assumes the disease to be endemic and given the long time period involved, it is not ruled out that new data would report cases of some type of periodic behavior.

The rest of this paper is organized as follows. In Sect. 2, we propose the SEIR model 

We consider the following SEIR model:

The model equations (Eqs. (1)-(4)) are based on the classical "Susceptible-Exposed-Infectious-Removed" (SEIR) model for a population of size (N ). The model attempts to describe what may unfortunately be an endemic disease (which may persist in a population for long time). Therefore the vital dynamics (births and natural deaths) have to be incorporated in the model. The variable (S) is the fraction of susceptible individuals (those able to contract the disease), (E) the fraction of exposed individuals (those who have been infected but are not yet infectious), (I) the fraction of infective individuals (those capable of transmitting the disease) and (R) is the fraction of removed individuals (those who have recovered or deceased). The model assumes that recovered individuals do not revert to the susceptible class. It is also assumed that all newborns are susceptible with the birth rate set equal to the death rate which is assumed not to be related to the infectious disease. The term βI N represents the force of infection where β is the effective per capita contact rate of infected individuals. The incidence rate is therefore βIS N . The parameter b is the rate of natural birth, α is the rate at which the exposed individuals become infective, so 1 α represents the mean latent period. The term 1 γ represents the mean infectious period. We have added to the classical SEIR model a new equation (Eq. (5)) and a new variable (P) that mimics the public perception of the severity of the pandemic. It can be seen that the dynamics (Eq. (5)) of the public perception of the risk regarding the pandemic is proportional to the number of infected cases (I) with e being the proportion of severe cases and 1 λ the mean duration of public reaction.

Besides the new variable (P), we adopt in this paper a new expression for the transmission rate β that reflects the impact of governmental action and the public perception of the severity of the disease. A number of studies have considered specific forms of the transmission rate. Lin et al. [6] , for instance, adopted the following expression for the transmission rate which was based on the formulation of He et al. [20] :

where D is a state variable representing social behavioral dynamics. The first term in Eq. (7) incorporates the impact of governmental action. It is parameterized by μ and represents all actions, which can be modeled as a step function. The second term in β represents the decreasing contacts among individuals reacting to the severity of the pandemic. The parameter κ represents the intensity of the individual reaction.

Kwuimy et al. [19] , on the other hand, proposed the following expression:

Similarly to Eq. (7), D represents social behavioral dynamics, μ represents the strength of the government action and κ is the strength of public response. Both expressions (Eqs. (7)- (8)) reflect the fact that public reaction would increase when more people get infected, and would naturally diminish over time. In this paper we assume the following expression for the transmission rate:

While keeping the same formulation for the impact of governmental action (all actions which can be modeled as a step function), we have opted for an exponential function to reflect the decreasing contacts among individuals reacting to severity of the disease. Mathematically, the expression exp( -κP N ) can be considered as a good approximation of (1 -P N ) κ (Eq. (7)), especially if the values of P are very small compared to the total population N , as the numerical simulations will show in this paper.

The model at steady state has two equilibria: A trivial one (S = N , E = 0, I = 0, R = 0, P = 0) and a nontrivial one, which satisfies the following transcendental equation:

We have

In the absence of governmental action (μ = 0) and public reaction (κ = 0), the transmission coefficient is constant β = β 0 , and the nontrivial steady state can be solved readily to yield explicit relations for the model state variable I:

The other state variables S, E, R and P can be obtained accordingly.

In the following, we show that the model solutions are positive under non-negative initial conditions. Proof Let x(t) = (S(t), E(t), I(t), R(t), P(t)) be the solution of system under initial conditions x 0 = (S(0), E(0), I(0), R(0), P(0)) = (S 0 , E 0 , I 0 , R 0 , P 0 ) ≥ 0. By continuity of the solution, for all S(t), E(t), I(t), R(t) and P(t) that have a positive initial value at t = 0, we have the existence of an interval (0, t 0 ) such that S(t), E(t), I(t), R(t), P(t) ≥ 0 for 0 < t < t 0 . We will prove that t 0 = ∞.

If S(t 1 ) = 0 for t 1 ≥ 0 and other solutions stay positive at t = t 1 , then

This implies that whenever the solution x(t) touches the S-axis, the derivative of S is nondecreasing and the function S(t) does not cross to negative values. Similarly, when E(t 1 ) = 0 for a t 1 > 0 and the other solutions stay positive,

When I(t 1 ) = 0 for a t 1 > 0 and the other solutions stay positive,

when R(t 1 ) = 0 for a t 1 > 0 and the other solutions stay positive,

Finally, when P(t 1 ) = 0 for a t 1 > 0 and the other solutions stay positive,

Therefore, whenever x(t) touches any of the axes S = 0, E = 0, I = 0, R = 0, P = 0, it never crosses them. Now, let N(t) = S(t) + E(t) + I(t) + R(t), we can see that

so the value of N is constant.

Next, we study the stability of the trivial solution. The Jacobian matrix evaluated at X 0 (S, E, I, R, P) = (N, 0, 0, 0, 0) is

The eigenvalues of the Jacobian matrix J(S, E, I, R, P) evaluated at the disease-free equilibrium are

where

It follows that the eigenvalues element (Eq. (19)) are negative if It should be noted that the results of Eq. (20) can also be derived by obtaining the expression of the basic reproduction number R 0 , which is an important parameter for the monitoring of the spread of the disease but also for the stability of the disease-fee solution. It is known that values of R 0 > 1 indicate an unstable equilibrium [21] . The derivation of the expression for R 0 is carried out in the appendix and yields

The disease-free equilibrium is stable provided that R 0 < 1, which is equivalent to Eq. (20) .

In this section we study the conditions of the occurrence of Hopf points in our fivedimensional model. We recall the conditions for a five-dimensional system to exhibit Hopf points [22] . We address the following characteristic equation:

This polynomial has exactly one pair of imaginary roots, χ 1,2 = ± √ θ , if and only if one of the following sets of conditions is satisfied:

The Jacobean J of the model is given by The term in the first Hopf condition C 1 (Eq. (23)) can be shown to be a, a 3 , a 4 , I, b, y) 

The obtained expressions for b i (i = 1, 5) and are quite complicated and are not amenable to analytical manipulation; therefore numerical simulations will be carried out in the next section. However, we can obtain an important result for the case when public reaction is absent (κ = 0). In this case it can be seen from Eq. (28) that ϒ = 0 and therefore b 4 > 0 (28), which contradicts the Hopf conditions of C 2 and C 3 (Eqs. (24)-(25)). As to the first Hopf condition C 1 (Eq. (23)), it can be seen that for (κ = 0) the term (Eq. (30)) is reduced to c 1 c 2 c 3 c 4 which is always strictly positive, which contradicts the Hopf condition.

The numerical analysis of the model is carried out using standard bifurcation techniques [23] with the help of the software AUTO [24] . The nominal values of the model parameters are listed in Table 1 . These parameters correspond to physically realistic values pertinent to COVID-19 disease. Figure 1(a) shows a typical behavior using the transmission rate β 0 as the main bifurcation parameter. There are two static solutions in the diagram: the 11)). The enlargement of the figure shown in Fig. 1(b) indicates that only values of β 0 greater than the value of Eq. (20) will lead to nontrivial solution. The diagram ( Fig. 1(a) ) is also characterized by the existence of one Hopf point. The periodic branch emanating from the Hopf point can be seen to terminate as it collides with the static branch. For β 0 larger than the value of Eq. (20) and up to the HB point the system will settle on the nontrivial solution, but for β 0 larger than the HB point, periodic behavior is expected in the model for a wide range of β 0 . An example of a limit cycle is shown in Fig. 2 for β 0 = 70. At this point, it is useful to show the effect of the different model parameters on the existence of periodic behavior. Each curve of Fig. 3 shows the locus of the Hopf point. Figure 3(a) shows the effect of governmental action. It can be seen that the model cannot exhibit a Hopf point in the absence of governmental action (i.e. μ = 0). Moreover, an increase in the strength of the governmental action (i.e. increase in the value of μ), decreases the range of periodic behavior, as the Hopf point moves to larger values of β 0 . The effect of the public response (κ) on the Hopf point is shown in Fig. 3(b) . An increase in the strength of the public reaction (larger κ) increases the range of periodic behavior as the Hopf point occurs at small values of β 0 .

The rest of the graphs of Fig. 3 shows that the range of periodic behavior decreases with smaller values of the latent period of the disease (i.e. larger α) or smaller values of the mean duration of public reaction (i.e. larger λ). On the other hand, the periodic behavior increases with the increase in the birth rate (larger b) or we have an increase in the rate of recovery (larger γ ).

Another type of bifurcation behavior predicted by the model can be shown in Fig. 4 where α (the inverse of the mean latent period) is chosen to be the main bifurcation parameter with the rest of model parameters set at their values of Table 1 with β 0 = 50. In the diagram the appearance of two Hopf points can be seen that are connected by a stable periodic branch. Unlike the case of Fig. 1 , the periodic behavior in this case is confined to a range of values of α. Figure 5 shows an example of limit cycle for α = 0.25. The effect of the different model parameters on the existence of such periodic behavior is shown in Fig. 6 where the two curves show the loci of the two Hopf points. It can be seen that an increase in the strength of the government action (μ) decreases the range between the two Hopf points and therefore decreases periodic behavior. The Hopf points will not exist beyond a critical point. On the contrary, increasing the strength of the public reaction (κ) Fig. 6 shows that an increase in the basic transmission rate (β 0 ) would increase the range of periodic behavior, and no periodic behavior can be found below a critical value of β 0 . The effect of the mean duration of public reaction shows a closed loop, which means that periodic behavior is confined within two critical values of λ. The same can be said about the effect of the birth rate b. Periodic behavior is expected only within a specific range of birth rate. Finally, when γ increases (i.e. the mean infection periodic decreases) the range of periodic behavior increases. The periodic behavior cannot exist if γ falls below a critical value.

This paper has proposed and analyzed the stability of a SEIR model structured upon susceptible, exposed, infected and removed cases, in addition to a new social behavior variable that mimics the public perception of risk regarding the severity of the pandemic. We also proposed a new expression of the transmission rate that modeled the impact of governmental action as a step function, and the individual reaction as a decreasing exponential function. The model was shown to predict one and two Hopf points. There is a fundamental difference between the two predicted periodic behaviors. While the former would exist for a wide range of model parameters, the latter is generally confined between some critical values. Regardless of the type of periodic behavior (emanating from one or two Hopf points) it was found that periodic behavior will increase in range if the disease has a large latent period, or if the mean duration of the public reaction increases, or the birth rate is high or the rate of recovery increases. Both the governmental action and public reaction have strong effects on the periodic behavior. A periodic behavior would not exist if no governmental action is taken or if there is no individual reaction. The range of periodic behavior would increase with a decrease in the strength of the governmental action or an increase in the strength of the public reaction.

The computation of the basic reproduction [21] number starts by identifying the infected compartments, i.e. E and I,

where F denotes the rate of appearance of new infections and V denotes the rate of transfer of individuals between compartments. Calculating the derivative of F and V with respect to x = (E, I), respectively, then substituting for initial values (S 0 , E 0 , I 0 , R 0 , P 0 ) yields

R 0 is the largest eigenvalue of matrix FV -1 i.e. R 0 = ρ(FV -1 ). Substituting for β = β 0 (1μ) exp( -κP N ) yields

Understanding Infectious Disease Dynamics

A modified SEIR model for the spread of Ebola in Western Africa and metrics for resource allocation

Extension and verification of the SEIR model on the 2009 influenza A (H1N1) pandemic in Japan

An SEIR model of influenza A virus infection and reinfection within a farrow-to-finish swine farm

Estimation of basic reproduction number of the Middle East respiratory syndrome coronavirus (MERS-CoV) during the outbreak in South Korea

A conceptual model for the coronavirus disease 2019 (COVID-19) outbreak in Wuhan, China with individual reaction and governmental action

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

An updated estimation of the risk of transmission of the novel coronavirus (2019-nCov)

Modified SEIR and AI prediction of the epidemics trend of COVID-19 in China under public health interventions

Oscillatory phenomena in a model of infectious diseases. Theor

Analysis and control of an SEIR epidemic system with nonlinear transmission rate

An Introduction to Mathematical Epidemiology, 1st edn

Some epidemiological models with nonlinear incidence

The Hopf bifurcation analysis and optimal control of a delayed sir epidemic model

Hopf bifurcation of a delayed epidemic model with information variable and limited medical resources

Seasonality in epidemic models: a literature review

Oscillations and chaos in epidemics: a nonlinear dynamic study of six childhood diseases in Copenhagen

Persistence, chaos and synchrony in ecology and epidemiology

Nonlinear dynamic analysis of an epidemiological model for COVID-19 including public behavior and government action

Inferring the causes of the three waves of the 1918 influenza pandemic in England and Wales

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

Complete coefficient criteria for five-dimensional Hopf bifurcations, with an application to economic dynamics

Dynamics of the Chemostat: A Bifurcation Theory Approach

Auto: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations

The authors would like to thank the anonymous referees for their valuable suggestions, which have greatly helped in improving the presentation of this paper. The first author would like to thank the Research Center of the College of Engineering at King Saud University, Riyadh, Saudi Arabia, for generous support.

Without loss of generality we can take N = 1. The initial values are S 0 = 1 and P 0 = 1 for the disease-free equilibrium. This yieldsFunding No funding available.

Not applicable.

The authors declare that they have no competing interests.

The two authors worked together in the derivation of the mathematical results. All authors read and approved the final manuscript. 

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 20 July 2020 Accepted: 22 September 2020