Skip to main content

ORIGINAL RESEARCH article

Front. Phys., 13 May 2021
Sec. Statistical and Computational Physics

Dynamics of an SIR-Based COVID-19 Model With Linear Incidence Rate, Nonlinear Removal Rate, and Public Awareness

\nAbdelhamid Ajbar
Abdelhamid Ajbar1*Rubayyi T. AlqahtaniRubayyi T. Alqahtani2Mourad BoumazaMourad Boumaza1
  • 1Department of Chemical Engineering, College of Engineering, King Saud University, Riyadh, Saudi Arabia
  • 2Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia

The paper studies the dynamics of the classical susceptible-infectious-removed (SIR) model when applied to the transmission of COVID-19 disease. The model includes the classical linear incidence rate but considers a nonlinear removal rate that depends on the hospital-bed population ratio. The model also includes the effects of media on public awareness. We prove that when the basic reproduction number is less than unity the model can exhibit a number of nonlinear phenomena including saddle-node, backward, and Hopf bifurcations. The model is fitted to COVID-19 data pertinent to Saudi Arabia. Numerical simulations are provided to supplement the theoretical analysis and delineate the effects of hospital-bed population ratio and public awareness on the control of the disease.

1. Introduction

The severity and the global scale of the COVID-19 pandemic have pushed research in many areas including the modeling of the disease dynamics with the goal of using such models to better understand the effects of intervention strategies on the disease control [15]. Several factors are known to affect the disease dynamics including the incidence rate, recovery rate, quarantine strategy, and awareness effects. Classical susceptible-infectious-removed (SIR) epidemic models were used extensively in the literature for modeling infectious diseases. These models are structured on the numbers of susceptible (S), infectious (I), and removed (recovered or deceased) individuals (R). When these classical models are used with the standard linear incidence rate and linear recovery rate they are known to predict at most one endemic static equilibrium and cannot exhibit bistability and periodic behavior. In fact, the model dynamic behavior depends mainly on the basic reproduction number. The disease is controlled and subsequently eliminated if the basic reproduction number is less than one and keeps spreading otherwise [68]. However, in practice infectious diseases can exhibit periodic oscillations and other nonlinear phenomena during the outbreak. Therefore, several studies attempted to model such nonlinear behavior. In this regard, nonlinear incidence rates were shown to play a key role in producing rich dynamics in epidemic models including periodic behavior [913].

The recovery rate is also an important parameter in the SIR model. Many of the studied SIR models in the literature assumed a constant recovery rate with the implicit assumption that health resources are sufficient to tackle the disease. However, the COVID-19 pandemic has shown that health services even in the most developed countries were overwhelmed during the first wave of the pandemic. It is therefore more realistic to consider a nonlinear recovery rate that incorporates factors such as hospitals capacity and efficiency of treatment. In this regard, various studies have analyzed epidemic models that incorporated nonlinear expressions of the recovery rate [14, 15]. Shan and Zhu [14], in particular, showed that the SIR model with linear incidence rate and nonlinear recovery rate can exhibit rich nonlinear phenomena such as saddle-node, backward, Hopf, and Bogdanov-Takens bifurcations.

Another parameter that can also affect the dynamics of the disease is public awareness and the role played by the media in raising such awareness. The media continuously alert the susceptible individuals regarding the severity of the disease and the need to take preventive measures which can lead to a decrease in the disease transmission [16, 17].

In this paper we propose to re-investigate the simple SIR model which is based on the standard incidence rate. We propose two modifications to the model: we assume a nonlinear removal rate that depends on hospital beds availability, and we model the media effects on public awareness in the form of a nonlinear rate with a saturation constant. We show that the proposed model can predict backward and Hopf bifurcations. Numerical simulations are provided to supplement the theoretical proofs. These simulations make use of parameter values fitted to COVID-19 data in Saudi Arabia.

A note should be made about the important issue of model validation. One objective of this work is to identify the parameters of the proposed model using real data pertinent to Saudi Arabia and then use the computed parameter values to analyze the static and dynamic behavior of the model. This parameter estimation is known as the inverse modeling problem. This inverse problem is important for model calibration and to allow public health officials to make reliable predictions about the pandemic. Various approaches were used in the literature for this purpose. These include least squares method that minimizes the sum of squared residuals [2] and variational method applied to time-depended model parameters [18]. Comunian et al. [19] provided a critical analysis of the inverse problem. The inverse problem is known to lack uniqueness of solutions [20] and also to lack the continuous dependence of the parameters to be identified on the data, so that any small errors in the data can lead to large discrepancies in the parameters to be identified. Moreover, the fitting can provide a tool for extracting some key quantities of interest, such as the basic reproduction number, that are not explicit in the epidemiological data streams. Unfortunately, as pointed out rightfully in [21] most epidemiological data streams are not designed for modeling. The challenge for data collection and fitting becomes more severe when the proposed model is highly structured. The balance between model structure and quality of fit was made clear in a number of studies concerning Saudi Arabia, where simple SIR models were shown to provide better fits to COVID-19 data than more complex network-based models [22, 23].

The rest of the paper is structured as follows. The model is presented in section 2, then in section 3 we study the model positivity, boundedness of its solutions and derive the basic reproduction number. Section 4 discusses the existence of equilibria while section 5 covers the stability of the disease-free equilibrium. In section 6, we study the bifurcations of the model. Numerical simulations are presented in section 7 followed by concluding remarks.

2. The Dimensional Model

The basic SIR model under study is described by the following set of equations:

dSdt=μN-β1SIN-β2SI2(M+I)N-μS    (1)
dIdt=β1SIN+β2SI2(M+I)N-(γ+μ)I    (2)
dRdt=γI-μR    (3)

In these equations, S and I denote the numbers of susceptible and infectious individuals. R, on the other hand, represents the removed (recovered or deceased) individuals. In order to simplify the model we did not include a separate equation for the deceased individuals. We assume a linear incidence rate β1I with β1 representing the transmission rate. The term β2IM+I represents the decrease in the disease transmission as result of increased public awareness due to media. This awareness can lead the public to take preventive measures. The parameter M represents the half-saturation constant. We assume that β1 is always larger than β2.

A rigorous functional description of public awareness is quite challenging. The expression used to quantify public awareness due to media was adopted from [17], but it is not unique. Other expressions were proposed in literature. Liu et al. [12] modeled the psychological effects of media by assuming the transmission rate β1 to be an exponential decreasing function that depends on numbers of reported infectious individuals. Feng et al. [16], on the other hand, used a simple constant parameter to quantify the effects of media on public awareness. The parameter γ represents the removal rate. Since the model considers R to represent both the recovered or deceased cases, γ also accounts for the fatality rate i.e., the rate at which infected people decease. Unlike many epidemic models where γ was assumed constant, we assume in this work that γ depends explicitly on the hospital-bed to population ratio (B) and the number of infectious individuals (I). The parameter (B) is generally considered by health authorities to be a key indicator of the adequacy of public health services [24]. The following relation which was first used in [14] is adopted in this work:

γ(B,I)=γ0+(γ1-γ0)BI+B    (4)

γ1 is the value taken by γ(B, I) when B → ∞ and/or the number of infectious individuals is very small i.e., I → 0. γ1 therefore represents the maximum removal rate when the number of hospital beds is more than sufficient and/or when there are few infectious individuals. γ0 is, on the other hand, the value taken by γ(B, I) when B → 0 and/or the number of infectious individuals is very large i.e., I → ∞. γ0 therefore represents the minimum per capita removal rate that can be sustained in the face of a large number of infections. Both values of γ0 and γ1 also depend on the inherent characteristics of the specific disease (COVID-19). In fact, when B is large, then the health system can provide a good recovery rate, so that the removed individuals R includes mostly recovered people. On the other hand, when B is small, then the health system cannot provide care to infected persons, and removed persons are mostly deceased persons. In other words, γ0 should be close to the “natural” fatality rate.

3. Model Analysis

The model is rendered partially dimensionless by defining:

s=SN,  i=IN,  r=RN,  m=MN,  b=BN

Equations (1–3) become:

dsdt=μ-β1si-β2si2(m+i)-μs    (5)
didt=β1si+β2si2(m+i)-(γ+μ)i    (6)
drdt=γi-μr    (7)

with

γ=γ0+(γ1-γ0)bi+b.    (8)

The initial conditions of the system (5–8) are given by: s(0) = s0 > 0, i(0) = i0 ≥ 0, r(0) = r0 ≥ 0.

In the following we prove that the proposed model is well-posed by showing that it satisfies the positivity condition.

It can be noted that by adding Equations (5–7) we can deduce that the time derivative of n(t): = s(t) + i(t) + r(t) is always zero and therefore n(t) = 1 for every t and the solutions are therefore always bounded.

3.1. Positivity of Solutions

We start by showing that the model solutions are always positive for non-negative initial conditions.

Theorem 1. Let s0, i0, r0 ≥ 0. The solution of Equations (5–8) with (s(0), i(0), r(0)) = (s0, i0, r0) is non-negative, that is s(t), i(t), r(t) ≥ 0, for t > 0.

Proof: The first equation (Equation 5) can be written as

ds(t)dt=μ-s(t)(β1i(t)-β2i(t)2(m+i(t))-μ),ds(t)dt=μ-s(t)y(t)    (9)

where

y(t)=β1i(t)-β2i(t)2(m+i(t))-μ.

Integrating Equation (9) we obtain the following expression:

s(t)=(0tμe0τy(u)dudτ+s0)e-0ty(τ)dτ>0.

This implies that s(t) is positive for all t. From Equation (6) we have:

di(t)dt-(γ+μ)i(t).    (10)

Integrating Equation (10) we obtain:

i(t)=i0e-0t(γ+μ)dτ>0.

It can also be shown in the same way that

dr(t)dt-μr(t)

which gives

r(t)=r0e-0tμdτ>0.

We conclude therefore that the model solutions are always positive.

4. Existence and Classification of Equilibria

In this section we study the number of real and positive equilibrium solutions of the system (Equations 5–8). The model has always a disease-free equilibrium given by E0(s, i, r) = (1, 0, 0). The non-trivial equilibria are obtained by setting the right sides of equations (Equations 5–7) to zero. The equilibrium condition of (7) is given by r=γμi and Equations (5) and (6) are independent of r, so that they can be treated as an independent sub-system. From Equation (6) we have:

s=((γ0+μ)i+b(γ1+μ))(m+i)(i+b)((β1+β2)i+β1m)    (11)

Substituting Equation (11) into Equation (5) leads to the following cubic polynomial:

E1(i)=a3i3+a2i2+a1i+a0=0    (12)

where a0, a1, a2 and a3 are defined by

a3=-(γ0+μ)(β1+β2),a2=-(γ1+μ)(β1+β2)b-(γ0+μ)β1m     +(β1+β2-γ0-μ)μ,a1=(-β1(γ1+μ)m+μ(β1+β2-γ1-μ))b     +μ(β1-γ0-μ)m,a0=μbm(R0-1)

where

R0=β1γ1+μ    (13)

is the basic reproduction number derived following the techniques presented in [25].

We can see from Equation (12) that the coefficient a3 is always negative. The coefficient a0, on the other hand, is positive when R0>1 and negative when R0<1. The possible number of positive roots can be obtained by the Descartes rule of signs as shown in Table 1.

TABLE 1
www.frontiersin.org

Table 1. Number of positive roots of Equation (12).

Theorem 2. The system (5-7)

1. has a unique steady state solution if R0>1 and whenever cases 1, 5, and 7 are satisfied;

2. can have more than one steady state solution if R0>1 and case 3 is satisfied;

3. can have two steady state solutions if R0<1 and cases 2, 4, and 6 are satisfied.

The results of Table 1 can be described qualitatively using bifurcation diagrams. It can be seen from the table that for R0<1 [cases (2,4,6,8)] there can be either two solutions or no solution. For R0>1 there can be only one solution. The possibility of three solutions (case 3) can be ruled out since it can be shown that the cubic equation (Equation 12) for case 3 can not predict three positive solutions. Therefore, only two expected behavior can be found as far as the bifurcation diagram (R0-i) is concerned. Figure 1A shows the case of two endemic equilibria for R0<1 and one endemic equilibrium for R0>1. Figure 1B shows, on the other hand, the case of the existence for R0<1 of only the disease-free equilibrium (no endemic solution) and one endemic equilibrium for R0>1.

FIGURE 1
www.frontiersin.org

Figure 1. Qualitative diagrams showing the two possible bifurcation behavior expected in the model. (A) backward bifurcation; (B) forward bifurcation (The stability nature of the branches is ignored).

The steady state multiplicity for R0<1 (Figure 1A) indicates the possible existence of a backward bifurcation where a stable endemic equilibrium coexists with the disease-free solution.

Before proceeding to the stability analysis, another useful relation can be obtained for the case of R0<1. The analysis of the coefficients ai(i = 0, 3) (Equation 12) shows that the following inequalities are always satisfied:

β1-β2-γ1-μ<β1-β2-γ0-μ<β1-γ0-μ                          since γ1>γ0  and  β2>0.    (14)

Therefore, if the last term in Equation (14) is negative i.e., β1 ≤ μ + γ0 the other terms of Equation (14) are also negative and therefore a3 and a2 are always negative which means that no positive solution can exist under this condition.

5. Stability of the Disease-Free Equilibrium

In this section we study the stability of the disease-free equilibrium point E0.

Theorem 3. For the nonlinear dynamical system (Equations 5–6)

 1. if R0<1 then E0 is an attracting node;

 2. if R0>1 then E0 is a saddle;

 3. if R0=1 then E0 is non-hyperbolic, and

(a) if a11=(β2m+γ1-γ0b-β12μ)0, E0 is a saddle-node.

(b) if a11 = 0 and a22=(β2γ1μm+β1(γ1-γ0)μb-β12γ1μ2)<0 then E0 is an attracting node.

Proof: The Jacobian matrix J(E0) of the system (Equations 5–6) at E0 is given by

J(E0)=[-μ-β10β1-γ1-μ].    (15)

The eigenvalues of the matrix Equation (15) are λ1=β1-γ1-μ=(γ1+μ)(R0-1) and λ2 = −μ. Therefore, we can see that if R0 < 1 then E0 is an attracting node and is a saddle if R0>1.

In the case of R0=1 the first eigenvalue is zero (λ1 = 0) and E0 is a non-hyperbolic point. In order to examine the nature of E0 we make use of the center manifold theory. We first transform the model (Equations 5–6) into the origin. Let u = s−1 and v = i then we have:

dudt=μ-β1(u+1)v-β2(u+1)v2m+v-μ(u+1)    (16)
dvdt=β1(u+1)v+β2(u+1)v2m+v-(γ0+(γ1-γ0)bv+b+μ)v.    (17)

Expanding the system (Equations 16–17) in a second order Taylor series around (u, v) = (0, 0) yields

dudt=-μu-β1v-β1uv-β2v2m    (18)
dvdt=(β1-γ1-μ)v+β1uv+(β2m+γ1-γ0b)v2.    (19)

When R0=1 then β1 − γ1 − μ = 0 and the system (Equations 18–19) becomes:

dudt=-μu-β1v-β1uv-β2v2m    (20)
dvdt=β1uv+(β2m+γ1-γ0b)v2.    (21)

The Jacobian matrix of the system is

J=[-μ-β100]    (22)

with eigenvalues

[0-μ].    (23)

We can therefore rewrite the system (Equations 20–21) as follows:

dxdt=a11x2+β1xy,    (24)
dydt=-μy+a22x2+a33xy    (25)

with

a11=(β2m+γ1-γ0b-β12μ),a22=(β2γ1μm+β1(γ1-γ0)μb-β12γ1μ2),a33=γ1β1μ.

If a11 ≠ 0 i.e., b(β12m-β2μ)μm(γ1-γ0) then it can be seen from (Equations 24–25) that E0 is a saddle-node.

If a11 = 0 i.e., b=(β12m-β2μ)μm(γ1-γ0) then we restrict the system (Equations 24–25) to the center manifold parameterized by y=h(x)=c1x2+c2x3 where c1 and c2 are parameters that need to be determined. The dynamics on the center manifold can be shown to be described by:

x=β1a22μx3+higher ordersy=-μy+higher orders.

Therefore, E0 is stable if a11 = 0 and a22 < 0. Substituting the condition a11 = 0 into a22 yields

a22=β12-β2μm.    (26)

The condition a22 < 0 requires that

β12<β2μm,    (27)

otherwise point E0 is unstable.

6. Bifurcation

6.1. Backward Bifurcation

In this section we set up the conditions under which the model exhibits a backward bifurcation.

Theorem 4. The system (Equations 5–6) exhibits a backward bifurcation whenever b<bcr=γ1-γ0(γ1+μ)2μ-β2m and no backward bifurcation otherwise.

Proof: We examine the existence of a backward bifurcation using the methodology described in [26]. First, we make the transformation of variables as follows: x1 = s, x2 = i. The model (Equations 5–6) is rewritten as follows:

dx1dt =f1=μ-β1x1x2-β2x1x22(m+x2)-μx1    (28)
dx2dt =f2=β1x1x2+β2x1x22(m+x2)-(γ+μ)x2.    (29)

We assume the bifurcation parameter to be β1. The condition R0=1 corresponds to β1=γ1+μ:=β1*. For β1=β1*, the Jacobian matrix of the system (Equations 28–29) at the disease-free equilibrium E0 = (1, 0) is given by

J=[-μ-β10β1-β1*].    (30)

The right eigenvector of this matrix at the zero eigenvalue is W=(w1,w2)T=(-β1*μ,1). Similarly, the left eigenvector is given by V = (v1, v2) = (0, 1). The stability parameters (a) and (b) that determine the conditions for the existence of a backward bifurcation can be found by the methodology explained in [26]. The coefficient (a) is given by:

a=k,i,j = 12vkwiwj2fkxixj,    (31)

which is reduced (since v1 = 0 and v2 = 1) to

a=w1w12f2x12+2w1w22f2x1x2+w2w22f2x22.    (32)

The different derivatives evaluated at E0 are given by

    2f2x12=02f2x1x2=β1*    2f2x22=2(γ1-γ0)b+2β2m.

Substituting yields

a= 2(γ1-γ0)b+2β2m-2β1(μ+γ1)μ    (33)

On the other hand, the bifurcation parameter b is given by

b= k,i = 12vkwi2fkxiβ1    (34)

which is reduced to

b=w22f2x1β1+w22f2x2β1,    (35)

and can be shown to be equal to 1.

We conclude therefore that since b = 1 is always positive then the system exhibits a backward bifurcation when a > 0. The condition a > 0 is equivalent to

γ1-γ0b>(γ1+μ)2μ-β2m    (36)

It can be seen that in the absence of media effects i.e., β2 = 0 the condition Equation (36) can be satisfied for some values of hospital-bed to population ratio (b). However, in the presence of media effects β2 ≠ 0 and if the right hand side of the inequality is negative i.e.,

β2m>(γ1+μ)2μ,    (37)

then the inequality (36) is always satisfied for any values of b. This means that media awareness also plays a critical role in the occurrence/absence of a backward bifurcation. If, on other hand, the inequality (37) is not satisfied then a backward bifurcation occurs when

b<γ1-γ0(γ1+μ)2μ-β2m    (38)

6.2. Hopf Bifurcation

In this section we study the Hopf bifurcation of Equations (5−−6). The Jacobian matrix is

J(s,i)=[J11J12J21J22],    (39)

where

J11=-[β1i+β2i2m+i+μ],J21=-[β1s+β2si(2m+i)(m+i)2],J21=β1i+β2i2m+i,J22=β1s+2β2sim+i-β2si2(m+i)2+(γ1-γ0)bi(i+b)2-γ0       -(γ1-γ0)bi+b-μ.    (40)

We use the fact that along the steady state solution we have

β1i+β2i2m+i+μ=μs,β1s+β2sim+i-γ0-(γ1-γ0)bI+b-μ=0.

Thus, we have

                                   J22=β2sim(m+i)2-(γ0-γ1)bi(i+b)2,Trace(J)=J11+J22=-(β1i+β2i2m+i+μ)+β2sim(m+i)2                                           -(γ0-γ1)bi(i+b)2.    (41)

We can write this equation as follows:

Trace(J)=-Φ(i)(m+i)2(i+b)2μ    (42)

where

Φ(i)= B5i5+B4i4+B3i3+B2i2+B1i+B0    (43)

with

B5=μ(β1+β2),B4=μ2+2(β1+β2)(b+m)μ+β2γ0m,B3=μ(β1+β2)b2+β1μm2+mμ(2μ-β2)         +((4β1μ+γ0β2+γ1β2+4β2μ)m         +μ(γ0-γ1+2μ))b,B2=m12μ2+((2β1μ+γ1β2+2β2μ)m+μ2)b2         +(2β1μm2+2μ(γ0-γ1+2μ-β2)m)bB1=(β1μm2+μ(2μ-β2)m)b2+μ(γ0-γ1+2μ)m2b1B0=b2m2μ2.    (44)

The coefficients of Equation (44) are always positive when the following condition holds

β2γ0-γ1+2μ    (45)

Simple analysis shows that the maximum value of γ1 when the right side of the condition (45) is positive satisfies

γ1,max2μ+γ0.    (46)

The other condition of Hopf point Det(J) > 0 can also be readily derived. However, both Hopf conditions are not easily amenable to analytical manipulations and therefore numerical simulations are used to illustrate the Hopf-induced dynamics of the model.

7. Numerical Simulations

In order to carry out numerical simulations we needed first to choose carefully the values of model parameters that reflect the disease parameters in Saudi Arabia, a country that has recorded more than 362,000 cases as on December 25, 2020 [27]. Only two model parameters were assumed constant: the natural death rate μ which is 4.21 × 10−5 [28] and the hospital-bed to population ratio which is 22 per 10,000 [29]. The rest of model parameters (m, β1, β2, γ0, γ1) were fitted to COVID-19 data in the country. The fitting was carried out with the help of the optimization function (fmincon) of MATLAB. The COVID-19 data used for fitting consist of infected cases and removed (recovered and deceased) cases. The data was collected for a period of 90 days starting from March 2 (the day of the detection of the first COVID-19 case in the country) until June 2. The data is publicly accessible via the COVID-19 dashboard of the Saudi Ministry of Health [27]. The reason we restricted the data to this period is to be able to compare the obtained reproduction number with similar results available in the literature [30, 31].

Figure 2 shows the fitting results. Both the infected cases (Figure 2A) and removed cases (Figure 2B) show a reasonable quality of fit. Since the optimization task led to local optima, we were guided by the results reported in [22, 30, 31] in selecting initial conditions and bounds on the fitted model parameters. The optimum values of model parameters are shown in Table 2.

FIGURE 2
www.frontiersin.org

Figure 2. (A,B) Results of fitting the model to COVID-19 data in Saudi Arabia for a period of 90 days starting from March 2, 2020. Blue solid line (model predictions); red points (actual data).

TABLE 2
www.frontiersin.org

Table 2. Model parameters.

For bifurcation studies the basic reproduction number R0 was taken to be the main bifurcation parameter so the transmission rate β1 can be calculated by β1=R0(γ1+μ) (Equation 13). The bifurcation diagrams were obtained using the continuation software AUTO [32]. Using the obtained values of model parameters, the condition (Equation 38) for a backward bifurcation to occur can be shown to be b < 2.01 × 10−4. Figure 3 shows an example of bifurcation diagrams in the (R0-i) plane for b = 6.322 × 10−7. This is the actual value of the dimensionless hospital-bed to population ratio, and can be seen to be much smaller than the aforementioned critical value. Saddle-node bifurcation can be seen to occur in this case. A static limit point LP appears at R0=0.803. The diagram is also characterized by the presence of a Hopf point occurring at R0=0.896. An unstable limit cycle can be seen to emerge from the Hopf point and disappears from homoclinic bifurcation. Here we have a situation where the disease-free equilibrium coexists with the endemic equilibrium for the range of basic reproduction number extending from HB (R0=0.896) to (R0=1).

FIGURE 3
www.frontiersin.org

Figure 3. Bifurcation diagram for b = 6.322 × 10−7 showing a backward bifurcation. solid line (stable branch); dashed line (unstable branch); circles (unstable periodic branches); LP (static limit point); HB (Hopf point); blue line (endemic equilibrium); red line (disease free equilibrium).

Figures 4A,B show numerical simulations for the Hopf point and the limit cycle. Figure 4A shows the birth of a limit cycle close to the Hopf point while Figure 4B shows an example of time trace of an unstable limit cycle at (R0=0.899). The oscillations grow in magnitude until they terminate as they collide with the static branch (S = 1).

FIGURE 4
www.frontiersin.org

Figure 4. (A) Time trace close to the Hopf point R0=0.896 of Figure 3. (B) Time trace for R0=0.899 of Figure 3 showing unstable oscillations growing in magnitude until they collide with the static branch.

From an epidemiological perspective, the existence of a Hopf point in the model and the nature of the Hopf bifurcation (subcritical) has two effects on the dynamics of the disease transmission: (1) The endemic equilibrium coexists now with the disease-free solution for values of R0 ranging from the Hopf point to R0=1. This means that even if the basic reproduction number is reduced below unity, it is still possible for the pandemic to persist unless the reproduction number is reduced below the point associated with the Hopf point (HB, R0=0.896). (2) The existence of a subcritical Hopf bifurcation (as opposed to supercritical when a stable limit cycle surrounds an unstable equilibrium point) means that the pandemic cannot have a stable periodic behavior and this is consistent with the current situation in the country which showed a steady decline in the disease active cases. Even if initial conditions (S, I) would fall close to the limit cycle, the resulting unstable oscillations would increase in magnitude but they will ultimately collapse and lead to the disease suppression (Figure 4B).

Figure 5 shows, on the other hand, the bifurcation diagram for b = 2.5 × 10−4, a value just larger than the aforementioned critical point. The backward bifurcation disappears and we have the occurrence of a forward bifurcation where the endemic equilibrium is the only stable point. A forward bifurcation is of course the outcome to aim for in a pandemic since if the basic reproduction number is reduced to a value below unity, then the disease outbreak is suppressed. Figures 3, 5 therefore showed the important role the number of hospital beds plays in the control of the pandemic. However from a practical point of view, the current hospital-bed to population ratio (22 per 10,000 people) is much smaller than the critical value. This means that a forward bifurcation is unlikely.

FIGURE 5
www.frontiersin.org

Figure 5. Bifurcation diagram for b = 2.5 × 10−4 showing forward bifurcation; solid line (stable branch); dashed line (unstable branch); blue line (endemic equilibrium); red line (disease free equilibrium).

Next, we examine the effects of variations in model parameters on the occurrence of the backward bifurcation. This can be suitably done by showing the location of the Hopf point when these parameters are varied, as shown in Figures 6A–E obtained also using AUTO. These figures are quite useful since they illustrate the range of the backward bifurcation, that is the range in term of R0 of the stable coexistence of the endemic equilibrium with the disease-free equilibrium. This sensitivity analysis is also necessary given the unavoidable degree of uncertainty in the fitted values of model parameters.

FIGURE 6
www.frontiersin.org

Figure 6. (A–E) Locus of the Hopf point of Figure 3 in terms of the model parameters.

We start with the effects of media shown in Figures 6A,B. The media effects are represented by two parameters: the awareness rate β2 and the half-saturation constant m. It can be seen from Figure 6A that the range of backward bifurcation decreased in term of R0 since the Hopf point moves to large values. The effect of m (Figure 6B) shows that as m increases the Hopf point moves to larger values of R0 which means again a smaller range of coexistence with the disease-free equilibrium. We can therefore conclude that the media also have an important role in attenuating the backward bifurcation.

Figure 6C shows, on the other hand, the effects of the minimum removal rate γ0. It can be seen that larger values of γ0 decrease the range of the coexistence since the Hopf point moves to larger values and closer to 1. The effect of the maximum removal rate γ1 is shown in Figure 6D. Unlike γ0, it can be seen that as γ1 decreases, the Hopf point moves to larger values of R0. While it is true that the values of γ0 and γ1 depend on the inherent properties of the virus causing the COVID-19 disease, they also depend on health care and basic clinical resources [14]. Their role is therefore important in attenuating the backward bifurcation. Finally, Figure 6E shows that as the hospital-bed to population ratio b increases, the location of the Hopf point moves to larger values of R0. This confirms the already obtained results that the backward bifurcation is reduced with the increase in available hospital beds.

8. Conclusions

The paper examined the dynamics of a SIR model applied to COVID-19 disease with two modifications: (1) the removal rate was assumed to be nonlinear and function of both the number of infectious individuals and hospital-bed to population ratio; (2) the model included public awareness due to media. The model was validated using COVID-19 data pertinent to Saudi Arabia. The model was shown to exhibit backward and Hopf bifurcations. The Hopf point is subcritical and unstable limit cycles terminate homoclinically. Therefore no stable periodic behavior occurred in the model for the chosen set of parameters. The main behavior observed in the model is a static coexistence of the disease-free equilibrium with the endemic point for a range of the basic reproduction number extending from the Hopf point to R0=1. The analysis also showed the important role played by health care resources in term of hospital-bed to populations ratio as well as the minimum and maximum removal rates. In fact with the country's current health resources, a backward bifurcation is the most likely outcome and can be attenuated by the increase in the media awareness.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors extend their appreciation to the Deanship of Scientific Research at King Saud University for funding this work through research group no (RG-1441-188).

References

1. Chen TM, Rui J, Wang QP, Zhao ZY, Cui JA, Yin L. A mathematical model for simulation the phase-based transmissibility of novel coronavirus. Infect Dis Pov. (2020) 9:24. doi: 10.1186/s40249-020-00640-3

CrossRef Full Text | Google Scholar

2. Bentout S, Chekroun A, Kuniya T. Parameter estimation and prediction for corona virus disease outbreak 2019 (COVID-19) in Algeria. AIMS Public Health. (2020) 7:306–18. doi: 10.3934/publichealth.2020026

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Belgaid Y, Helal M, Venturino E. Analysis of a model for Corona virus spread. Mathematics. (2020) 8:820. doi: 10.3390/math8050820

CrossRef Full Text | Google Scholar

4. Mishra AM, Purohit SD, Owolabi KM, Sharma YD. A nonlinear epidemiological model considering asymptotic and quarantine classes for SARS CoV-2 virus. Chaos Soliton Fract. (2020) 138:109953. doi: 10.1016/j.chaos.2020.109953

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kwuimy CAK, Nazari F, Jiao X, Rohani P, Nataraj C. Nonlinear dynamic analysis of an epidemiological model for COVID-19 including public behavior and government action. Nonlinear Dyn. (2020) 101:1545–59. doi: 10.1007/s11071-020-05815-z

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Brauer F, Castillo-Chavez C. Mathematical Models in Population Biology and Epidemiology. New York, NY: Springer-Verlag (2001).

Google Scholar

7. Diekmann O, Heesterbeek JAP. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. New York, NY: John Wiley & Sons Ltd. (2000).

Google Scholar

8. Hethcote HW. The mathematics of infectious disease. SIAM Rev. (2000) 42:599–653. doi: 10.1137/S0036144500371907

CrossRef Full Text | Google Scholar

9. Ruan S, Wang W. Dynamical behavior of an epidemic model with a nonlinear incidence rate. J Differ Equat. (2003) 188:135–63. doi: 10.1016/S0022-0396(02)00089-X

CrossRef Full Text | Google Scholar

10. Tang Y, Huang D, Ruan S, Zhang W. Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate. SIAM J Appl Math. (2008) 69:621–39. doi: 10.1137/070700966

CrossRef Full Text | Google Scholar

11. Hu Z, Ma W, Ruan S. Analysis of SIR epidemic models with nonlinear incidence rate and treatment. Math Biosci. (2012) 238:12–20. doi: 10.1016/j.mbs.2012.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Liu R, Wu J, Zhu H. Media/psychological impact on multiple outbreaks of emerging infectious diseases. Comput Math Methods Med. (2007) 8:612372. doi: 10.1080/17486700701425870

CrossRef Full Text

13. Wang W. Backward bifurcation of an epidemic model with treatment. Math Biosci. (2006) 201:58–71. doi: 10.1016/j.mbs.2005.12.022

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Shan C, Zhu H. Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds. J Differ Equat. (2014) 257:1662–88. doi: 10.1016/j.jde.2014.05.030

CrossRef Full Text | Google Scholar

15. Cui Q, Qiu Z, Liu W, Hu ZY. Complex dynamics of an SIR epidemic model with nonlinear saturate incidence and recovery rate. Entropy. (2017) 19:35. doi: 10.3390/e19070305

CrossRef Full Text | Google Scholar

16. Feng L, Jing S, Hu S, Wang DF, Huo HF. Modelling the effects of media coverage and quarantine on the COVID-19 infections in the UK. Math Biosci Eng. (2020) 17:3618–36. doi: 10.3934/mbe.2020204

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Mohsen A, Al-Husseiny HF, Zhou X, Hattaf K. Global stability of COVID-19 model involving the quarantine strategy and media coverage effects. AIMS Public Health. (2020) 7:587–605. doi: 10.3934/publichealth.2020047

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Marinov CT, Marinova RS. Dynamics of COVID-19 using inverse problem for coefficient identification in SIR epidemic models. Chaos Soliton Fract. (2020) 5:5100041. doi: 10.1016/j.csfx.2020.100041

CrossRef Full Text | Google Scholar

19. Comunian A, Gaburro R, Giudici M. Inversion of a SIR-based model: a critical analysis about the application to COVID-19 epidemic. Phys D. (2020) 413:132674. doi: 10.1016/j.physd.2020.132674

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Giordano G, Blanchini F, Bruno R, Colaneri P, Di Filippo A, Di Matteo A, et al. Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy. Nat Med. (2020) 26:855–60. doi: 10.1038/s41591-020-0883-7

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Keeling MJ, Dyson L, Guyver-Fletcher G, Holmes A, Semple MG, Tildesley MJ, et al. Fitting to the UK COVID-19 outbreak, short term forecasts and estimating the reproductive number. medRxiv. (2020). doi: 10.1101/2020.08.04.20163782

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Alharbi Y, Alqahtani A, Albalawi O, Bakouri M. Epidemiological modeling of COVID-19 in Saudi Arabia: spread projection, awareness, and impact of treatment. Appl Sci. (2020) 10:5895. doi: 10.3390/app10175895

CrossRef Full Text | Google Scholar

23. Alrasheed H, Althnian A, Kurdi H, Al-Mgren H, Alharbi S. COVID-19 Spread in Saudi Arabia: modeling, simulation and analysis. Int J Environ Res Public Health. (2020) 17:7744. doi: 10.3390/ijerph17217744

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Boaden R, Proudlove N, Wilson M. An exploratory study of bed management. J Manage Med. (1999) 13:234–50. doi: 10.1108/02689239910292945

CrossRef Full Text | Google Scholar

25. Van den Driessche P, Watmough J. Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Castillo-Chavez C, Song B. Dynamical models of tuberculosis and their applications. Math Biosci Eng. (2004) 1:361–404. doi: 10.3934/mbe.2004.1.361

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Saudi Arabia COVID-19 Dashboard. Available online at: https://covid19.moh.gov.sa/ (accessed November 2, 2020).

28. General Authority for Statistics. Available online at: https://www.stats.gov.sa/en/5305 (accessed November 2, 2020).

29. Saudi Ministry of Health. Available online at: https://www.moh.gov.sa/en/Pages/Default.aspx (accessed November 2, 2020).

Google Scholar

30. Aletreby W, Alharthy A, Faqihi F, Mady AF, Ramadan OE, Huwait BM, et al. Dynamics of SARS-CoV-2 outbreak in the Kingdom of Saudi Arabia: a predictive model, Saudi. Crit Care J. (2020) 4:79–83. doi: 10.4103/sccj.sccj_19_20

CrossRef Full Text | Google Scholar

31. Alshammari F. A Mathematical model to investigate the transmission of COVID-19 in the Kingdom of Saudi Arabia. Comput Math Method M. (2020) 2020:9136157. doi: 10.1155/2020/9136157

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Doedel EJ, Kernevez JP. Auto: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations. Pasadena, CA: Institute of Technology (1999).

Google Scholar

Keywords: SIR model, variable recovery rate, hospital beds, backward bifurcation, Hopf bifurcation

Citation: Ajbar A, Alqahtani RT and Boumaza M (2021) Dynamics of an SIR-Based COVID-19 Model With Linear Incidence Rate, Nonlinear Removal Rate, and Public Awareness. Front. Phys. 9:634251. doi: 10.3389/fphy.2021.634251

Received: 27 November 2020; Accepted: 06 April 2021;
Published: 13 May 2021.

Edited by:

Carla M. A. Pinto, Instituto Superior de Engenharia do Porto, Portugal

Reviewed by:

Dawei Ding, Anhui University, China
Mauro Giudici, University of Milan, Italy
Miguel A. F. Sanjuan, Rey Juan Carlos University, Spain

Copyright © 2021 Ajbar, Alqahtani and Boumaza. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Abdelhamid Ajbar, aajbar@ksu.edu.sa

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.