From the early attempts [1–4] to recent studies, epidemic modeling which is applicable in a wide range of fields from informatics [5, 6] to chemistry [7–9] has drawn the attention of researchers in various disciplines. Since the basic compartmental model Susceptible-Infected-Removed (SIR) which is commonly used to model diseases for which the infection confers permanent immunity was introduced by Kermack and McKendrick in 1927 [4]; other compartmental models [10–12] have been developed to model diseases with different structures and dynamics. Especially in recent years, major outbreaks such as avian flu in 2005, swine influenza in 2006, and H1N1 influenza in 2009 have highlighted the need for more effective and reliable models to control the spread of disease and to provide a better knowledge for the prediction of future threats and for the development of stronger containment strategies. In Refs. [13–18], some results on the modeling of different types of epidemic diseases, the solution form of these models, the observation of global stability, and the determination of the final size of the epidemic are obtained. In Ref. [15], global stability criteria are derived for the SEIS model which can be regarded as a model with no immunity, and different SIR models are examined in Refs. [16, 17]. Moreover, the works [17, 18] provide some useful results on the final size of the epidemic for SIR models. With the same motivation of these works but rather a different contribution to literature, we use a multistage model [14] in this article to explain the time shift observed in several surveys such as Spanish flu (1918–1919) [19], SARS (2002–2004), the 2009 H1N1 in Istanbul, and recently, COVID-19 [20] (see Figures 1–3).

In the literature, available observed data that are used for the ordinary differential equation system representing the classical SIR epidemic model is based on the curve of removed individuals. Usually, this curve is obtained by taking into consideration only the fatality data of the epidemic disease, whereas in some research studies, not only the fatality data but also the hospitalization data for the epidemic are taken into account in the modeling process. In Ref. [21], it is shown that there exists a delay between the peak of the hospitalization (infectious) curve and the inflection point of the fatality (removed) curve based on the data collected. The original contribution of this article to the literature is that we explain this time shift by the multidimensional form of SIR and SEIR models and also provide numerical evidence that the expected delay is approximately half of the infectious period of the epidemic disease for both of the multistage systems.

In the second section of this work, we give a brief summary of the classical epidemic SIR and SEIR models and define the multistage form of these models that will be used in further analysis. The graphs obtained by the numerical evaluations of the classical SIR and SEIR models and their multistage forms are also given. Analysis of these graphs confirms that the multistage SIR and SEIR models explain the time shift observed in several surveys. In the third section, the evaluation of delay for different epidemic parameters is presented by using the numerical evaluations of these multistage models. In the fourth section, the distance between the points where successive stages and hence any two stages assume their maximum is found approximately. The last section includes a summary of the results obtained in the previous sections as well as motivations for future analysis.

The SIR model is commonly used to model diseases for which the removed individuals are assumed to be immune to reinfection. In addition, the total population with constant size is divided into three distinct compartments the size of which change with time t. These compartments are called the susceptible class S, the infective class I, and the removed class R. Healthy individuals with no immunity are members of class S until they are infected with a pathogen and become capable of transmitting the disease to others. They move from the class S into the class I once they are infected and then from I to R once they recover or die. Childhood illnesses like measles or rubella are good examples for the SIR model.

The SIR epidemic model without vital dynamics, that is, the recruitment of new susceptible through birth or immigration as well as the loss through mortality or emigration are ignored, is defined by the following system of nonlinear ordinary differential equations: S I R :         S ′ = − β S I , I ′ = β S I − γ I , R ′ = γ I , (1) where the coefficient β refers to the disease transmission rate and 1 / γ represents the duration of infection period. Note that since S ′ + I ′ + R ′ = 0 , we may assume S + I + R = 1 by the use of appropriate normalization.

The standard SIR model ignores a latent phase which is the delay between the time of the acquisition of infection and the onset of infectiousness. In order to define this latent phase, the introduction to the SIR model of an exposed class E whose members are individuals who have been infected with a pathogen but are not yet infectious due to the incubation period of pathogen yields the SEIR model. Chicken pox is suitable for the SEIR model, which is defined by the following system of ordinary differential equations S E I R :         S ′ = − β S I , E ′ = β S I − ϵ E , I ′ = ϵ E − γ I , R ′ = γ I , (2) where 1 / ϵ represents the mean exposed period. Note that we may assume S + E + I + R = 1 by the use of appropriate normalization.

These two models are suitable for mathematical modeling of seasonal diseases, but they fail to reproduce the time shift that was observed in the modeling of the 2009 H1N1 epidemic in Istanbul, Turkey [21], as shown in Figure 2. For COVID-19, a similar time shift is observed between the curves of total cases and removed individuals in China, South Korea, Iran, Turkey, Germany, and Brazil as shown in Figure 3. Publicly accessible data that have been released by the state offices of each country are used for this analysis.

In the literature, a similar time shift is also observed for Spanish flu [19] and SARS epidemics. The graphs for the data of these epidemics are shown in Figure 1. Weekly case and fatality reports for Copenhagen, Denmark, in 1918–1919 are displayed in Figure 1A. In this figure, it can clearly be seen that there is a time shift between the peak of the infectious cases and the peak of fatalities. Similarly, in Figure 1A–B, the time shift is also observable between cumulative cases and cumulative fatalities.

In this study, we use multiple infectivity periods [13, 14] to explain this delay that is unforeseen in the standard SIR model. The approach of multiple infectious stages consists of replacing the single infectious stage I with N + 1 substages denoted by I i , which is the density of individuals in the ith infectious stage. Unlike the model in Ref. [14], each of these stages may have different infectivity β i and a variable infectious period 1 / γ i . In order to compare the solution curves with the ones for the standard SIR and SEIR models, we set 1 γ 0 + 1 γ 1 + ⋯ + 1 γ N = 1 γ . The multistage SIR and SEIR epidemic models are defined by the following systems: M u l t i s t a g e     S I R :         S ′ = − S ( β 0 I 0 + β 1 I 1 + ⋯ + β N I n ) , I ′ 0 = S ( β 0 I 0 + β 1 I 1 + ⋯ + β n I N ) − γ 0 I 0 , I ′ 1 = γ 0 I 0 − γ 1 I 1 , ⋯   ⋯ ⋯ ⋯   ⋯ ⋯ I ′ N = γ N − 1 I N − 1 − γ N I N , R ′ = γ N I N , (3) M u l t i s t a g e    S E I R :         S ′ = − S ( β 1 I 1 + ⋯ + β N I N ) , E ′ = S ( β 1 I 1 + ⋯ + β N I N ) − ϵ E , I ′ 1 = ϵ E − γ 1 I 1 , I ′ 2 = γ 1 I 1 − γ 2 I 2 , ⋯   ⋯ ⋯ ⋯   ⋯ ⋯ I ′ N = γ N − 1 I N − 1 − γ N I N , R ′ = γ N I N . (4) The multistage SIR and SEIR systems with β 0 = … = β n and γ 0 = … = γ n correspond to the choice of gamma-distributed “Infection Period Distribution” (IPD) in the integral equation formulation of the SIR model [14].

The linear parts of apparently different infectious stages for i ≥ 1 in the multistage SIR model and i ≥ 2 in the multistage SEIR model have a similar structure. We write the linear parts of each equation above as a system and then rearrange and rename as follows to keep the models as clear and simple as possible S J R : J = I 0 + β 1 β 0 I 1 + ⋯ + β N β 0 I N , S ′ = − β 0 S J , I ′ 0 = β 0 S J − γ 0 I 0 , I ′ i = γ i − 1 I i − 1 − γ i I i ,   f o r   i = 1 , … N , R ′ = γ N I N , (5) S E J R : J = I 1 + β 2 β 1 I 2 ⋯ + β N β 1 I N , S ′ = − β 1 S J , E ′ = β 1 S J − ϵ E , I ′ 1 = ϵ E − γ 1 I 1 , I ′ i = γ i − 1 I i − 1 − γ i I i ,   f o r   i = 2 , … N , R ′ = γ N I N . (6) Subsequently, the numerical evaluations of these two systems, SJR and SEJR, defined above will be used for some of the structural comparisons of the classical models and multistage models. A MATLAB^® ODE45 solver is used for all the numerical evaluations. In typical cases, initial conditions are chosen so that S ( 0 ) is close to 1, whereas E ( 0 ) and I i ( 0 ) are close to zero, and R ( 0 ) is chosen so that the sum of all variables is 1, that is, S ( 0 ) = 1 − 10 − 4 ,   I 0 ( 0 ) = 10 − 4 ,     I i ( 0 ) = 0 ,     i = 1 , … . , N ,     R ( 0 ) = 0 ; (7) and for the SEJR model, S ( 0 ) = 1 − 10 − 4 ,   E ( 0 ) = 5 × 10 − 5 ,   I 1 ( 0 ) = 5 × 10 − 5 ,   I i ( 0 ) = 0 ,   i = 2 , … . , N ,   R ( 0 ) = 0. (8) Numerical evaluations of the classical SIR and SEIR models are made for specific epidemic parameters, and the results are given in Figure 4. For the evaluations of the classical models, γ and ϵ are fixed as 1 / 5 and 1 / 3 , respectively, whereas the value of R 0 is chosen to be 2.5 and 10, respectively. In all cases, the maximum of the infectious stage and the inflection point of the removed stage occur at the same point in time; hence, there is no time shift in these classical models.

Numerical evaluations of the multistage SIR model are repeated for various stage numbers. First, R 0 is set as 2.5, while the total number of stages, N, is given the values 5,10,30 , and 60, and the corresponding graphs of the solutions are shown in Figure 5. In these evaluations, γ i is chosen to be 0.2 × N , for i = 0,1 , … , N . In Figure 5, the time shift between the maximum point of the curve, representing the sum of the infectious stages, J ( t ) and the inflection point of R ( t ) can clearly be seen. As predicted, the system given by Eq. 5 confirms the delay and therefore seems adequate to explain the time shift between infectious and removed stages observed in Istanbul data [21]. Note that the time shifts for N = 10,30,60 seem to be equal but the one for N = 5 is smaller.

Similarly, numerical evaluations of the multistage S E I R model are obtained for various stage numbers. First, R 0 is set as five, while the stage number N is given the values 5,10,30 , and 60, and the corresponding graphs of the solutions are shown in Figure 6. In these evaluations, ϵ and γ i are chosen to be 1 / 3 and 0.2 × N , respectively, for i = 1 , … , N . Figure 6 illustrates that just like it is seen in the SJR system, there exists a time shift between the maximum point of J and the inflection point of R in the SEJR model, with characteristics similar to the ones for the SIR model.

In this section, we investigate the dependency of the delay on the epidemic parameters by numerical evaluations of system Eq. 5. Graphs I ( i ( t ) for N = 1,2,3,4 for the SIR and SEIR models are given in Figures 7 and 8, respectively. To analyze the nature of the time shift for relatively large values of the stage number, the graphs of the solutions for the same epidemic parameters are obtained for N ranging from 2 to 16 in steps of two and the corresponding graphs are given in Figure 9. As it can be seen from Figure 9, the solution curves start to resemble as the stage number N increases. Finally, in Figure 10, we present the dependency of the time shift on the number of stages N, for 1 ≤ N ≤ 150 . Similar computations are repeated for the SEIR model and the results are presented in Figures 11, 12.

For the SIR model, to investigate the effect of the basic reproduction number R 0 and the infectious period 1 / γ on the infectious dynamics and the resulting delay, system (5) is solved with the initial conditions given by Eq. 7 for some parameter values. To this end, the pair ( R 0 , 1 / γ ) is chosen ( 2.5 , 5 ) , ( 5,5 ) , ( 10,5 ) , ( 2.5 , 10 ) , ( 5,10 ) , ( 10,10 ) , ( 2.5 , 20 ) , ( 5,20 ) , and ( 10,20 ) , respectively, and the numerical evaluations for various infectious stages N are shown in Figure 10. It can be observed from this figure that the delay is almost half of the infectious period. This fact can also be seen in Figure 9 where the time value of the maximum of J / N (normalized J) in time is located at the middle of time values of the maximum of the first infectious stage and the maximum of the last infectious stage. Comparison of panels of Figure 10 shows that the change in the reproduction number R 0 for a fixed infection period has no effect on the delay. On the other hand, the delay depends on the infection period; in fact, it is approximately half of it for large N.

The same analysis is repeated for the multistage SEIR model. The system defined by Eq. 6 is solved for the same epidemic parameters and initial conditions as above and with ε = 1 / 3 , and the resulting graphs are given in Figure 11. As for the SIR model, the solution curves start to resemble for large N and the peak of J ( t ) is located at the midpoint of the delay interval.

To illustrate the dependency of the delay on the system parameters of the S E I R model, the pairs ( ε , γ ) are chosen as ( 1 / 3,1 / 3 ) ( 1 / 3 , ( 1 / 5 ) , ( 1 / 5,1 / 3 ) , ( 1 / 5,1 / 5 ) , and variations of the delay for each of these cases are presented in Figures 12A–D, for R 0 = 5 and R 0 = 7.5 . An analysis of these graphs yields that as N increases, the delay converges to a value. Moreover, it could easily be observed that the delay is independent of ϵ and R 0 , but yet it is influenced by 1 / γ (i.e., infectious period). Note that the delay for the multistage SEIR model is shorter than the delay observed in the multistage SIR model.

In this section, we use a quadratic approximation to each I i ( t ) around the peak of I i − 1 ( t ) to show that within the validity of the quadratic approximation, the delays between successive peaks are 1 / γ i .

Let t i be the time where each I i assumes its maximum and 1 / γ i be the corresponding infectious period for i = 0,1 , … N . To determine the distance between the points t i where each substage I i assumes its maximum value, we use quadratic approximation of the Taylor series expansion of I i at the point t = t i − 1 where, I i ( t ) = I i ( t i − 1 ) + I i ′ ( t i − 1 ) ( t − t i − 1 ) + 1 2 I i ′ ′ ( t i − 1 ) ( t − t i − 1 ) 2 , (9) for i ≥ 1 . Differentiating I i ′ ( t ) = I i ′ ( t i − 1 ) + I i ′ ′ ( t i − 1 ) ( t − t i − 1 ) (10) and then substituting t = t i in Eq. 10 and using the fact that I i ′ ( t i ) = 0 since I i reaches its maximum at t i , one obtains t i − t i − 1 = − I i ′ ( t i − 1 ) I i ′ ′ ( t i − 1 ) . (11) The multistage SIR model defined by the equations in Eq. 5 suggests that for i ≥ 1 , I i ′ ( t ) = γ i − 1 I i − 1 ( t ) − γ i I i ( t ) . (12) Differentiating Eq. 12 yields I i ′ ′ ( t ) = γ i − 1 I i − 1 ′ ( t ) − γ i I i ′ ( t ) . (13) By considering the fact that I i − 1 ′ ( t i − 1 ) = 0 since I i − 1 assumes its maximum value at t i − 1 , one gets the following by replacing t = t i − 1 in equation I i − 1 ′ ′ ( t i − 1 ) = − γ i I i ′ ( t i − 1 ) . Substitution of the equation above in Eq. 11 gives the approximate distance formula as follows: t i − t i − 1 = 1 γ i . (14) Therefore, the distance between any t i is t i − t j = ∑ k = j + 1 i 1 γ k .

Results obtained by the numerical evaluations are compatible with Eq. 14. To observe the distance between the maximum points of the independent infectious stages, solutions of the multistage SIR model with respect to various infectious periods are chosen. In this respect, the basic reproduction numbers R 0 and γ i (i.e., duration is 5) are set as 2.5 and 0.2 × N , respectively, and the related graphs of the solutions for various stage numbers ( N = 1,2,3,4 ) are given in Figure 7. Comparison of graphs in Figure 7 reveals that the value of the difference of the points where successive stages reach their maximum is approximately 1 / γ i .

It should be emphasized that Eq. 14 is also valid for the multistage SEIR model. The distance between the maximum points of the independent infectious stages including the E stage is approximately 1 / γ i . Since the proof is the same as in the SIR case, we do not repeat the derivation of Eq. 14 again to avoid repetition. However, to observe the distance between the maximum points of the infectious stages numerically, the basic reproduction numbers R 0 , ϵ, and γ i (i.e., duration is 5) are set as 5, 1 / 3 , and 0.2 × N . Then, the SEIR model is solved for N = 1,2,3 , and 4, and the related graphs are given in Figure 8. As in the case of the SIR model, it is observed that the distance between the maxima is found to be approximately 1 / γ i , too.

Epidemic data display a time shift between the peaks of infectious cases and fatalities. This time shift is not foreseen by the ordinary differential equations for the SIR and SEIR models since in both of them, the derivative of R ( t ) is proportional to I ( t ) . Nevertheless, this can be remediated by using gamma distributions instead of exponential distributions for the infectious period distribution (IPD) in the original SIR model of Refs. [4, 10] given in terms of integral equations [14], leading to a multistage model.

In this article, we propose a generalization of these multistage models by allowing the parameters to be unequal in different stages. We showed that within the validity of a quadratic approximation to I i ( t ) , the distance between the points where each infectious stage reaches its maximum is approximately 1 / γ i .

We solved the multistage models for a range of epidemic parameters, and we have seen that the solution curves reveal the time shift. While the delay varies for relatively small stage numbers, it is observed that the delay becomes nearly stable as the number of stages increases, and it is independent of the basic reproduction number. This fact supports the validity of the quadratic approximation and shows that the delay phenomenon observed in the infectious diseases defined by the epidemic models SIR and SEIR can be successfully explained by the multistage forms of these models.

In addition to a theoretical contribution, the existence and the estimation of the time shift between the progression of the infectious cases and the fatalities have a practical importance, in the sense that, in order to take timely actions, the severity of the epidemic should be measured in terms of the increase in the number of infectious cases.

Finally, we note that the importance of the effects of quarantine is realized during the COVID-19 pandemic. In the literature, there are two basic approaches, one of which adds a compartment to the model as in Ref. [22] and the other adds a function β / ( γ + Q ( t ) ) , where Q ( t ) is a time-dependent rate as in Ref. [23] and has the effect of shortening the duration of the infection period. This may explain the relatively shorter delays in China (compared with South Korea) where strict quarantine was in effect.

The datasets analysed for this study can be found in the webpages: http://www.worldometers.info/coronavirus/ (last access: 05 November 2020), http://epikhas.khas.edu.tr/ (last access: 05 November 2020) and https://www.who.int/csr/sars/country/en/ (last access: 05 November 2020).

APD, AHB and SA performed theoretical results; AD, APD, AHB and SA performed computations; AD performed literature survey and APD, AD, AHB and SA wrote the paper.

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.