Delayed Dynamics of SIR Model for COVID-19


This paper presents a new modified SIR model which incorporates appropriate delay parameters leading to a more precise prediction of COVID-19 real time data. The efficacy of the newly developed SIR model is proven by comparing its predictions to real data obtained from four counties namely Germany, Italy, Kuwait, and Oman. Two included delay periods for incubation and recovery within the SIR model produce a sensible and more accurate representation of the real time data. In the absence of the two-delay period () the dynamical behavior of the model will not correspond to today’s picture and lag the detection of the epidemic peak. The reproductive number R0 is defined for the model for values of recovery time delay of the infective case. The effect of recovery time may produce second wave, and/or an oscillation which could destabilize the behavior of the system and a periodic oscillation can arise due to Hopf bifurcation phenomenon.

Share and Cite:

Ebraheem, H. , Alkhateeb, N. , Badran, H. and Sultan, E. (2021) Delayed Dynamics of SIR Model for COVID-19. Open Journal of Modelling and Simulation, 9, 146-158. doi: 10.4236/ojmsi.2021.92010.

1. Introduction

The announcement of the global pandemic outbreak of coronavirus COVID-19 by the world health organization has become one of human’s most concern due to its enormous infectious diseases, both in terms of medical, and economics. The outbreak of the coronavirus COVID-19 has stimulated search for drugs and intensified scientists to search and understand the dynamics of spreading such pandemic. Mathematical modelling is an essential tool to understand the mechanism of spread of a disease such as COVID-19 in the human population. These models generate insights into the transmission dynamics of infectious diseases and assist health officials and policymakers to control its extensive spread. Mathematical models of infectious disease dynamics have been extensively analyzed for the past century [1] [2] [3] [4]. These models are based on SIR model. The controlling of the spread of a virus is based on dimensionless constant known as the reproduction value “R0” which measures the transmission potentials of a virus. The value for R0differs from one virus to another (ex., R0 value for influenza spread ≈ 3 [5]. When R0 < 1, an infectious infects less than one individual and the speed of the disease is expected to stop. However, R0 alone provides limited information about the transmission potential of infectious diseases [6] [7] [8], further theoretical work is needed to connect the patterns of an epidemic. Recent studies on COVID-19 [9] [10] [11] [12] [13] have been helpful to gain insights into the transmission dynamics and potential role of different intervention strategies such as mitigation and suppression to slow down the epidemic spread, reduce the peak health care to protect those who are most at risk from infections, and reducing the number of infective cases to the low level, as well as enforcing lockdown to region of highly infective cases, home isolation of suspect cases, home quarantine of those living in the same household, and implementing social distancing among individuals. COVID-19 is a newly emergent virus, and much remains to be understood about its transmission. Therefore, understanding the incubation and recovery period are very important for health authorities as it allows introducing more effective quarantine systems for individual suspected of carrying the virus.

The epidemic model presented in this document is not novel, perhaps similar to models presented for different diseases. To develop an exact dynamical model is extremely difficult if not impossible to COVID-19 specification because of the different prevention measures taken by different countries (forced isolation, lock-down, social distancing ...) to control and contain the spread of the disease. In the present study, an investigation is presented which shows that a modified SIR model that includes latency time constants resulted in a closely accurate estimate of the real time data for several countries. The accuracy of the suggested modified SIR model is shown to better forecast the number of invectives and recoveries of COVID-19, as well as provides more accurate indicators of how fast the infection is propagating. The benefit of accurately estimating the recovery/infectious rates is to predict a possible slowdown or growth of the infection numbers and allow public health policymakers to determine which containment measures are more effective decisions to take by the Government in combating the spread of the COVID-19 pandemic.

2. Construction the Model: The SIR Framework

In this section, we describe the mechanistic transmission model that enables us to understand and forecast the spread of an epidemic. The model is based on dividing the population into three distinct compartments. Each compartment is based on the infectious status of individuals in the population [1] [2] [3] [4]. Therefore, population is divided into three classes termed as: Susceptible “S”, Infected “I”, and Recovered “R”.

Susceptible individuals are assumed to never have been infected before; however, they can become infected through contacts with infectious individuals at rate proportional to constant β. Meanwhile, some infectious individuals recover and become immune at rate proportional to constant γ. Finally, the rest of the infected individuals migrate due to epidemic at rate proportional to constant α as shown in the schematic diagram Figure 1.

The model assumes that the population size is fixed (no birth), the recovered individuals receive total immunity, and the death rate is very small compared to the total invectives.

The infectious process is modeled by a set of ODE’s:

d S d t = β S ( t ) I ( t ) N (1)

d I d t = β S ( t ) I ( t ) N γ I ( t ) α I ( t ) (2)

d R d t = γ I ( t ) (3)

The total population size is assumed to remains constant over the period of the epidemic such that N = S + I + R is constant ( R = N S I ). Therefore,

d N d t = d S d t + d I d t + d R d t = 0

2.1. The Reproductive Ratio R0, and Equilibrium Points

The reproductive number, R0, is the most important quantity in epidemiology. It is the number of invectives produced by a primary infective in a fully susceptible population “virgin population”. It tells us about the initial rate of increase of the disease over a generation. In other words, Ro depends on the transmission dynamics of the disease and measures the growth potential of an infectious disease [14] [15]. The value of R0 remains invariant during the early epidemic growth dynamics. Hence, if R0> 1, the disease will increase and become eventually epidemic, and if R0< 1, the disease will die out and tend to zero, and if R0= 1, the disease is self-sustaining and the number of infected will remain constant. In the model we assume the population at time t = 0 has one infected individual and correspond to S ( 0 ) = N 1 , I ( 0 ) = 1 , R ( 0 ) = 0 . If one infected individual appears in the population, there will be an epidemic if and only if d I d t > 0 , otherwise, the rate of infective is decreasing d I d t < 0 . Using the initial values for S with N 1 , and I set to 1 in the SIR model represented by Equation (2) at t = 0 yields

Figure 1. Schematic diagram of the basic SIR model. Boxes represent compartment, and arrows indicates the rate of transfer between compartments.

β N 1 N ( γ + α ) < 0 .

Rearrange the inequality and shift the threshold to 1 as for the reproductive number gives

β γ + α N 1 N < 1 .

For huge populations, the term N 1 N approaches 1 which implies the disease will dissipate. On the other hand, when β γ + α > 1 the disease is expected to become a pandemic, the quantity β γ + α is defined as the reproductive number “R0”. To obtain the equilibrium points of the SIR model (1 - 3) are set to zero.

The equilibrium points are where the variables do not change with time. i.e. d S d t = d I d t = d R d t = 0 ; such that

d S d t = 0 = β S I N

d I d t = 0 = β S I N γ I α I

This gives two equilibria: ( S = S = N , I = 0 ) and ( S = 0 , I = 0 ), meaning that there are no infected humans. Shifting the variables so that the origin is at equilibrium ( N , 0 ) ( 0 , 0 ) :

S * = N S S = N S *

I * = I 0 I * = I

d S * d t = d N d t d S d t = β I * N ( N S * )

d I * d t = d I d t = β I * N ( N S * ) ( γ + α ) I *

Considering small deviations from equilibrium, so that S * and I * are small, and higher powers of S * and I * are neglected.

d S * d t β I *

d I * d t β I * ( γ α ) I *

The autonomous linear differential equations have solution of the form A e λ t . Substituting into the above equation reduces it to a standard exponential form as

d S * d t = S 0 λ e λ t = β I 0 e λ t .

d I * d t = I 0 λ e λ t = ( β γ α ) I 0 e λ t

Transforming the above exponential equations into a standard matrix form as follows

λ ( S 0 I 0 ) = [ 0 β 0 β γ α ] ( S 0 I 0 ) [ λ β 0 λ β + γ + α ] ( S 0 I 0 ) = 0 .

Setting the determinant of the matrix to zero, yields two eigenvalues λ = 0 , and λ = β γ α . If λ = β γ α > 0 , the solutions grow away from equilibrium. The equilibrium is unstable, and there is an epidemic ( R 0 = β γ α > 1 ).

If λ = β γ α < 0 , the solution contract back toward equilibrium. The equilibrium is stable, and there is no epidemic ( R 0 = β γ α < 1 ).

Define E 0 = ( N , 0 ) as the disease-free equilibrium points of system (1 - 3) [16]. The disease-free equilibrium of the system is locally asymptotically stable if the eigen values of the characteristic equation are negative. Hence, the equilibrium point is locally stable if R 0 < 1 , and unstable if R 0 > 1 .

2.2. Incubation and Recovery Time

A novel human coronavirus spread in December 2019 in Wuhan, China. It is an enigmatic and confusing illness, wrapped with uncertainty because there have not been enough scientific studies on how long an individual might have the symptoms or be contagious or totally immune. Despite the uncertainty, ranges have been identified and were found to vary from one individual to another. The current understanding of the incubation and recovery periods are limited. The mean incubation period observed to be 3 days (0 - 24 days) conducted on 1324 cases [17]. The recovery period tends to be 1 - 3 weeks depending on how mild or severe the disease is. However, these are rough guidelines. For example, the symptoms of mild illness could extend to 3 weeks [18].

3. Delayed SIR Model

This part of the paper is devoted to construct the dynamical model for our proposed problem, the disease is assumed to have an incubation period of the virus τ 1 > 0 , and recovery period τ 2 > 0 . The incubation period represents the delay time from exposure to the development of symptoms of the virus. The bilinear transmission incidence will be a function of ( t τ 1 ). The recovery period represents the delay time from being infected to getting totally immune and move to the susceptible compartment and will be a function of ( t τ 2 ). The process model dynamics can be described as

d S d t = β S ( t τ 1 ) I ( t τ 1 ) (4)

d I d t = β S ( t τ 1 ) I ( t τ 1 ) γ I ( t τ 2 ) α I ( t ) (5)

d R d t = γ I ( t τ 2 ) (6)

In the system of equation, a susceptible individual is assumed to interact with an infective individual and does not move to the infected compartment until after certain time “incubation period” as of the case of COVID-19. The incubation period “ τ 1 ” is only when moving from the susceptible compartment to the infected compartment. Similarly, “ τ 2 ” is the period for an infected individual moving from infected compartment to the recovery compartment. Our modification of the SIR model considers the delay constants, while a detailed description can be obtained by using a system of cities connected by traffic streams.

If the pandemic duration is long enough (over a year), the structure behavior of the modified SIR proposed model presented in (4 - 6) reduces to the classic SIR model presented in (1 - 3). That is the incubation and recovery periods “ τ ” are very small compared to the long duration time scale. It can be proven using tailor series expansion for small τ around the parameter t and ignoring the terms ( O τ , O τ 2 ) .

4. Numerical Simulation and Discussion

Simulations for dynamical system the classical SIR model (1 - 3) and the proposed time delayed SIR model (4 - 6) are compared to real data collected by the official site of World Health Organization (WHO) [19]. The data were collected for up to July 7th, 2020 for four countries (Germany, Italy, Kuwait, and Oman) that have tried different strategies toward controlling, mitigation and suppressing the disease. The data used are the original time series data which shows significant daily fluctuations. There are three stages for the spreading of the disease; the early stage of the spread, where the transmission of the disease is slow among the individuals, the middle stage, where the infection propagates very aggressively toward its maximum daily infective numbers, and the last stage, the epidemic starts to slow down and decay until eventually dies out. In what follows we compare both theoretical models using simulations including incubation and recovery durations with the obtained daily data for the infectives and cumulative recovery from WHO [19]. In the simulations performed, it is assumed that S(0) and I(0) are initial susceptible and infective individuals, respectively. Noting that S S ( 0 ) = N 2 , and I ( 0 ) = 2 N , are two infected individuals, and N is the population of the country. With initial conditions appropriately set, the model demonstrates self-similarities with the real time data. The parameters β, γ, and α are varied for optimal curve fitting and the results shown in figures with corresponding captions of each country.

4.1. Germany

Figure 2 simulate both models to real time infective daily cases and cumulative recovery. Figure 2(a) demonstrates a sensible smooth fit to the real data by

Figure 2. Top: Germany’s SIR curve compared with the modified time-delayed SIR model estimates with β = 0.17 , γ = 0.03 , α = 0.01 , R 0 = 4.25 , τ 1 = 3 and τ 2 = 9 vs. lower curve representing basic SIR model with τ 1 = τ 2 = 0 .

introducing two-time delays, and reproduction value R 0 = 3.7 . Where τ 1 is the time period effect for an individual moving from the susceptible compartment to the infected compartment and τ 2 is the time period effect for an individual moving from an infected compartment to the recovery compartment. The parameters β, γ, and α values are shown in the corresponding captions of each figure. They are varied to get the best fit in curvature with the real time data curve.

Figure 2(b) does not correspond to today’s picture, predicts a false alarm detection of the epidemic due to considering the incubation and recovery time ( τ 1 = τ 2 = 0 ), which makes it difficult to cope with the spreading of the disease. In the figure, different phenomena are observed significant spreading (sharp rise) as well as large daily fluctuations in the infective cases. Germany had its own epidemic plans focusing on containment, protection, and mitigation [20]. Despite Germany’s strong health care system and early progress on developing a tracing application to alert individuals had encounter people tested positive, imposing lock down, and the different precautionary measures to interrupt transmission of the disease during the severe spreading events (20/3-10/4). The situation persists and remained uncontrolled even after locked down and social distancing in attempt to reduce the reproduction number R 0 = 3.5 . The main cause for the high spread among the youth is that they did not abide by the safety protocols and advice provided by the local health authority [21].

4.2. Italy

It took three weeks for Italy to go from discovering the 1st case to the closure of all non-essential business, and all movement of individuals within the whole territory. Unfortunately, the government underestimated how fast the virus spread and how quickly it could push their health care system to the verge of collapse. The infection grew very fast and was highly lethal overcoming other major countries in the number of infected people as early as March 23rd reaching over 50,000 total infections [22]. Figure 3 shows numerical solution to both models respectively and compared to real time data for daily invective cases and cumulative recovery patients. The result shown in Figure 3(a) fit well with real data, and the curves are well constrained. Although the approximated reproductive number R 0 = 2.25 considered to be relatively low, the Italian outbreak was the worst in Europe with the actual reproductive number in the range R 0 = 2 - 3.4 . Mitigation strategies (social distancing, quarantine measures, lock down) are essential to contrast further epidemic spreading, especially in countries experiencing a lag time behind the Italian outbreak such as Kuwait and Oman. Figure 3(b) shows a delayed detection time of the epidemic due to not considering the incubation and recovery time ( τ 1 = τ 2 = 0 ).

Figure 3. Top: Italy’s SIR curve compared with the modified time-delayed SIR model estimates with β = 0.18 , γ = 0.03 , α = 0.05 , R 0 = 2.25 , τ 1 = 3 and τ 2 = 8 vs. lower curve representing basic SIR model with τ 1 = τ 2 = 0 .

4.3. Kuwait and Oman

The global Epidemic outbreak alarmed Gulf region with a lag time compared to Europe, causing the government to enforce control measures, and implement several strategies to reduce, interrupt, and/or control spreading of the disease (mitigation and suppression). Measurements were taken in early March to suspend non-essential governmental agencies, businesses, studies, and all flights. By mid of March, the government announced the country entered a transmission stage of coronavirus which led to enforcing curfews, lockdowns, and city isolation, as well as enforcing punishments and fines for violators of guidelines including imprisonment. Health authorities have relied on random testing (PCR) to identify the infected individuals to treat in the hospital, and the contacted ones are quarantined, and monitored. In Kuwait, the aim is to maintain the health care system by reducing the number of cases to low level or eliminate human to human transmission until a vaccine is available. In April 78% of spreading was in migrant cities, the government imposed zonal quarantine to five densely populated migrant worker cities to minimize the ongoing spread [23]. The reproductive number varied R 0 = 1.08 - 3.8 depending on type of action implemented by the government in the country. Figure 4 represent a numerical simulation for both models when the reproductive number R 0 = 1.82 . Figure 4(a), shows a simulation solution of model Equations (4)-(6) compared to real time data for daily infective cases and cumulative daily recovery. The inconsistency in the fitting to interpret the actual evolution of the infection is due to the complex procedure taken by authorities and the inhomogeneous daily number of tests that led to a biased estimation for COVID-19. Figure 4(b), shows a delayed detection of the epidemic due to setting τ 1 = τ 2 = 0 .

Figure 4. Top: Kuwait’s SIR curve compared with the modified time-delayed SIR model estimates with β = 0.208 , γ = 0.02 , α = 0.01 , R 0 = 1.8 , τ 1 = 3 and τ 2 = 10 vs. lower curve representing basic SIR model with τ 1 = τ 2 = 0 .

Figure 5 simulates both models to the real time data for daily infective cases and the cumulative daily recovery in Oman. Figure 5(a) demonstrate a smooth fit of the modified model with real time data and a reproduction number R 0 = 1.82 , while Figure 5(b) shows a delay epidemic when τ 1 = τ 2 = 0 . Oman reached the highest daily infective cases later than Kuwait (mid of July). Oman implemented similar strategies but at different timing, locked down, restricted movement and isolation were imposed to a densely populated province causing 45% of the total infective cases. Most of the cases 62% were in the expatriates [24].

Figure 5. Top: Oman’s SIR curve compared with the modified time-delayed SIR model estimates with β = 0.255 , γ = 0.05 , α = 0.09 , R 0 = 1.82 , τ 1 = 3 and τ 2 = 9 vs. lower curve representing basic SIR model with τ 1 = τ 2 = 0 .

4.4. Extending the Recovery Period

For all the above cases extending recovery period τ 2 will stretch the duration of the pandemic with another wave of disease and/or perhaps become an endemic.

Figure 6 shows a numerical simulation of model (2) with different values of the parameters given in the corresponding caption. In the figure, the negative region of infectives is ignored since it cannot be negative or perhaps immunity due to prior exposure. Although the model shows large epidemics occurring; nonetheless, the level of the disease reaches a plateau around a constant level. Figure 6(b) depicts the Hopf bifurcation of the periodic outbreak of the disease which destabilize the system meaning d I d t > 0 . The oscillation and the periodic behavior in Figure 6 do rely exclusively on the parameters of the model and not


Figure 6. The delayed SIR model undergoes (a) damping oscillation with a parameter β = 0.6 , γ = 0.05 , α = 0.08 , τ 1 = 3 , τ 2 = 37 . (b) Hopf bifurcation with a parameter β = 0.255 , γ = 0.05 , α = 0.01 , τ 1 = 3 , τ 2 = 37 .

necessarily on higher reproductive number. In Figure 6(a), a value of R 0 = 4.6154 produced a decaying oscillation, while varying the parameters of the system to reproduce a lower reproductive number R 0 = 4.25 resulted in a periodic waves of disease.

5. Conclusion

This paper demonstrates the effect of incorporating two delay periods within the SIR model as compared to real data provided by [19]. The delay periods correspond to the duration of the incubational and recovery periods as it appears in COVID-19. The values of the two-delay period were picked in accordance to recent literatures [17] [18]. The parameter variations in the model β , γ , α , τ 1 and τ 2 make up different cases corresponding to different situations. The variations are introduced based on what is being observed in the literature to predict what possible actions of the future are. The model is well fitted on epidemic real time data for four countries. Simulation results are consistent with data and generated curves are well constrained. The parameters are varied for producing a reproduction number within range of the countries. Despite the real time data shown in the figures are up to July 7th, 2020 updating the daily infective cases would still fit the model because the number of daily infectives is gradually decreasing. Therefore, extending the recovery period will result in double exposure and/or an epidemic, to avoid oscillations, the model should not exceed a critical threshold at the endemic equilibrium. Although the mortality rate is a minor fraction of the total infective cases, it is important to include it for a more accurate reproduction number. Furthermore, the future is murky because it depends on human actions, both individual and collective, which remains a challenge to be more accurately predicted and modeled.


The authors would like to thank Ms. Nourah Ebrahim for collecting the real time data and Dr. Ali Hajjiah for his valuable comments.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.


[1] Kermack, W. and McKendrick, A. (1927) A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society A, 115, 700-721.
[2] Hethcote, H.W. (2000) The Mathematics of Infectious Diseases. SIAM Review, 42, 599-653.
[3] Anderson, R.M. (1982) Population Dynamics of Infectious Diseases. Chapman and Hall, London.
[4] Chowell, G., et al. (2016) Mathematical Models to Characterize Early Epidemic Growth: A Review. Physics of Life Reviews, 18, 66-97.
[5] Anderson, R.M. and May, R.M. (1991) Infectious Diseases of Humans Dynamics and Control. Oxford University Press, Oxford.
[6] Davis, S., et al. (2008) The Abundance Threshold for Plague as a Critical Percolation Phenomenon. Nature, 454, 634-637.
[7] Salkeld, D.J., et al. (2010) Plague Outbreaks in Prairie Dog Populations Explained by Percolation Thresholds of Alternate Host Abundance. Proceedings of the National Academy of Sciences, 107, 14247-14250.
[8] Cross, P.C., et al. (2007) Utility of R0 as a Predictor of Disease Invasion in Structured Populations. Journal of the Royal Society Interface, 4, 315-324.
[9] Maleewong, M. (2020) Time Delay Epidemic Model for COVID-19.
[10] Menéndez, J. (2020) Elementary Time-Delay Dynamics of COVID-19 Disease.
[11] Ivorra, B., et al. (2020) Mathematical Modeling of the Spread of the Coronavirus Disease 2019 (COVID-19) Taking into Account the Undetected Infections. The Case of China. Communications in Nonlinear Science and Numerical Simulation, 88, Article ID: 105303.
[12] Arino, J. and Porteta, S. (2020) A Simple Model for COVID-19. Infectious Disease Modelling, 5, 309-315.
[13] Wang, N., et al. (2020) An Evaluation of Mathematical Models for the Outbreak of COVID-19. Precision Clinical Medicine, 3, 85-93.
[14] Heesterbeek, J.A.P. (2002) A Brief History of R0 and a Recipe for Its Calculation. Acta Biotheoretica, 50, 189-204.
[15] Blackwood, J.C. and Childs, L.M. (2018) An Introduction to Compartmental Modeling for the Budding Infectious Disease Modeler. Letters in Biomathematics, 5, 195-221.
[16] Porwal, P., Shrivastava, P. and Tiwari, S.K. (2015) Study of Simple SIR Epidemic Model. Advances in Applied Science Research, 6, 1-4.
[17] Guan, W.J., Ni, Z.Y., Hu, Y., et al. (2020) Clinical Characteristics of 2019 Novel Coronavirus Infection in China.
[18] Bai, Y., Yao, L., Wei, T., et al. (2020) Presumed Asymptomatic Carrier Transmission of COVID-19. JAMA, 323, 1406-1407.
[19] COVID-19 Coronavirus Pandemic.
[20] an der Heiden, M. and Buchholz, U. (2020) Modellierung von Beispielszenarien der SARS-CoV-2-Epidemie 2020 in Deutschland. Report, Robert Koch Institute, Berlin.
[21] Coronavirus SARS-CoV-2: Chronik der bisherigen Maßnahmen [Coronavirus SARS-CoV-2: Chronology of Measures Taken]. German Federal Ministry of Health.
[22] De Natale, G., Ricciardi, V., et al. (2020) The Covid-19 Infection in Italy: A Statistical Study of an Abnormally Severe Disease. Journal of Clinical Medicine, 9, 1564.
[23] Alkhamis, M.A., Al Youha, S., et al. (2020) Spatiotemporal Dynamics of COVID-19 Epidemic in the State of Kuwait. International Journal of Infectious Diseases, 98, 153-160.
[24] Zia, K. and Farooq, U. (2020) COVID-19 Outbreak in Oman: Model-Driven Impact Analysis and Challenges.
[25] Radha, M. and Balamuralitharan, S. (2020) A Study on COVID-19 Transmission Dynamics: Stability Analysis of SEIR Model with Hopf Bifurcation for Effect of Time Delay. Advances in Difference Equations, 2020, Article No. 523.

Copyright © 2021 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.