key: cord-102966-7vdz661d authors: Nikolaou, M. title: A Fundamental Inconsistency in the SIR Model Structure and Proposed Remedies date: 2020-05-01 journal: nan DOI: 10.1101/2020.04.26.20080960 sha: doc_id: 102966 cord_uid: 7vdz661d The susceptible-infectious-removed (SIR) compartmental model structure and its variants are a fundamental modeling tool in epidemiology. As typically used, however, this tool may introduce an inconsistency by assuming that the rate of depletion of a compartment is proportional to the content of that compartment. As mentioned in the seminal SIR work of Kermack and McKendrick, this is an assumption of mathematical convenience rather than realism. As such, it leads to underprediction of the infectious compartment peaks by a factor of about two, a problem of particular importance when dealing with availability of resources during an epidemic. To remedy this problem, we develop the dSIR model structure, comprising a single delay differential equation and associated delay algebraic equations. We show that SIR and dSIR fully agree in assessing stability and long-term values of a population through an epidemic, but differ considerably in the exponential rates of ascent and descent as well as peak values during the epidemic. The novel Pade-SIR structure is also introduced as a approximation of dSIR by ordinary differential equations. We rigorously analyze the properties of these models and present a number of illustrative simulations, particularly in view of the recent coronavirus epidemic. Suggestions for further study are made. In their landmark 1927 publication Contribution to the Mathematical Theory of Epidemics, 1, 2 Kermack and McKendrick developed a general, if elaborate model structure to capture the dynamics of a fixed-size population comprising compartments of individuals susceptible (S) to a spreading infection, infectious (I), and removed (R) from the preceding two compartments by recovery or death. Propagation of an individual from S to I to R underlies the basic context of the exercise. In a modeling tour-de-force, the authors eventually present, in equations (11) through (15) of their paper (ibid.), the general structure of the elaborate mathematical model they derive. They proceed to examine the implications of their model for special cases, and finally present a very special case resulting in a set of three relatively simple ordinary differential equations (ODEs, equations (29) ibid.), which were destined to form the basis for a genre of mathematical models in epidemiology, the celebrated SIR model and its many variants. 3 The three ODEs are where prime denotes time derivative. 2, 4, 5 Consistent with the importance of the SIR ODEs, the basic reproductive ratio, 0 ≝ / , is widely considered "one of the most critical epidemiological parameters" 2 and has even become a household name in the recent coronavirus epidemic. 6 In the sentence right before they present their SIR model in equation (29), (ibid., cf. the above eqns. The assumption about constant is plausible, as it refers to the rate of spread of the epidemic (cf. eqn. (1) ). While that parameter might change over time as a result of measures However, the rate of removal depends more on the duration over which individuals have remained infected and less on the size of that group. In the conceptually simple case where all individuals recover or die at a single number of days, , after their infection, the rate of removal is 1 or 0, depending on whether is greater than or not ( Figure 1 ). Of course, in reality will likely follow a distribution ( Figure 2 ) rather than being a single number. However, in that case as well, the removal rate will depend on comparison between and the distribution of , rather than on the size of the group remaining infectious for time . Starting with the assumption that individuals leave the infectious group at time after infection, we develop in this paper a corresponding mathematical model structure, named delay SIR (dSIR), in the form of a single delay differential equation (DDE) for , and two associated delay algebraic equations, for and in terms of . In the rest of the paper we first introduce the dSIR model structure in section 2, and provide an intuitive exposition of its basic properties in section 3, where we also introduce the Padé SIR model structure as an ODE approximation of dSIR. Rigorous analysis follows in section 4. In that section we explain that certain SIR and dSIR properties are exactly similar (e.g. herd immunity, total number of infected), while others are quite different (e.g. exponential rates, peak of the infectious group). Extension to models with additional compartments beyond S, I, and R is discussed in section 5, with presentation of the dSPIR model. Finally, the significance of this work and future extensions are discussed. To explain the derivation of the dSIR model structure, we will rely on the detailed version of Figure 2 shown in Figure 3 . The schematic shows the evolution of , , over discrete time steps = 0, 1, 2, … of length each. The thick-green bordered rectangle in Figure 3 suggests by visual inspection that where ( , ) is the rate at which the infection spreads, for example in proportion to the product , as shown above. Assuming that each new part of the infectious fraction, , gained at a given time step, , moves to the removed fraction, , in time steps, the thick-red bordered rectangle in Figure 3 suggests that Eqns. (4) and (5) yield where ′ ≝ / and − ≝ ( − ). Eqn. (7) is a single nonlinear DDE that involves only and is decoupled from equations for and . As such, it captures the entire dynamics of the dSIR system. medRχiv.org Michael Nikolaou 3 The remaining two fractions of the population, and , can simply be inferred by algebraic equations, as eqn. (5) implies and + + = 1 implies To put the schematic in Figure 3 in context, observe that summation of eqns. (1)-(3) and discretization yields This suggests that, according to the SIR model, each thickblack bordered green rectangle in Figure 3 would have to be a constant fraction of the immediately previous orange column. However, with each new part of the infectious fraction moving to the removed fraction after steps (thickblack bordered orange rectangles moving to thick-black bordered blue rectangles) this cannot be true as an emerging fact, by the simple observation of the shapes of and + , which follow a convex to concave pattern after an inflection point. Therefore, the assumption of constant , equivalent to constant in eqn. (10) , is not compatible with the assumption of time to transition from infectious to removed being independent of the infectious fraction size. This will be further illustrated with simulations in the next section. Incidentally, eqns. (1)-(3) of the SIR model can also be decoupled to a single if non-intuitive nonlinear ODE in terms of , solved in a form implicit in time, . 7 Before any theoretical properties of the dSIR model structure are analyzed, a simple visual comparison between SIR and dSIR is presented. Unless specifically stated otherwise, all numerical simulations have been conducted with 8,9 Figure 4 is the standard { , , } plot, for SIR and dSIR. Note the faster dynamics of dSIR compared to SIR, the approximately twice as high peak of i for dSIR compared to SIR, as well as the asymmetric profile of for SIR, compared to the symmetric dSIR profile of . (Details in section 4.) To further illustrate the dSIR/SIR relationship, Figure 5 and Figure 6 show the continuous-time counterparts of Figure 3 for the SIR and dSIR models, respectively. In this stacked representation of { , , } over time, it is clear that the SIR model corresponds to a time-varying infectious period for each individual. In fact, the value of calculated according to eqn. (3) as with ′ and produced by the dSIR model, is time varying, as shown in Figure 7 . This discrepancy suggests that the standard interpretation of 1/ as the average infectious period … estimated relatively precisely from epidemiological data 2 is increasingly inaccurate as 0 increases above 1. What is estimated from epidemiological data is the delay , rather than , and if 1/ is used as an estimate of , as is typically done, the SIR model response will be too slow, with a peak value for lower than its dSIR counterpart by about half. Depending on . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 1, 2020. (9) . Note that remains constant. Figure 5 is also superimposed, for comparison. The correspondence = 1/ was used. (9) . Note that all lines start at = , as ′ ( < ) = 0. Note also small spurious deviations from constant shortly after = , due to numerical approximation of ( ≤ 0), by the continuous function As pointed out in the preceding sections, interpreting the value of in the SIR model as the average infectious period is not accurate and produces misleading results. It turns out (Appendix A) that the following simple remedy can be used to retain the ODE structure of the standard SIR model, while better approximating the DDE dynamics of the more realistic dSIR model structure: The SIR equations for { ′ , ′ }, eqns. (2) and (3), can be replaced by the equally simple ODEs or by the more accurate set of ODEs with 2 (0) = 0. Figure 8 confirms that the trajectories for the dSIR and the modified SIR model structures -eqns. (7)-(9) and (1)-(3) with eqn. (2) replaced by (13) or (14), respectively -are close to one another, both in terms of the time to peak and the value of the peak. As already mentioned, these properties are important when such models are used to anticipate hospitalization needs for the infected during an epidemic. Figure 8 . Comparison of the profiles for the first-order Padé SIR, second-order Padé SIR, dSIR, and SIR models. Note the improve approximation of dSIR by the second-order Padé SIR, compared to first-order Padé SIR, as anticipated by eqns. (13) or (14), respectively. Note (Appendix A) that the essence of eqns. (13) and (14) is in approximating the pulse profile ( ) − ( − ) in Figure 2 by a transfer function approximation based on firstor second-order Padé-approximants (in the Laplace domain) (rather than by the decaying exponential of Figure 2 ) as shown in Figure 9 . Padé-approximation has long been a popular . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. . https://doi.org/10.1101/2020.04.26.20080960 doi: medRxiv preprint medRχiv.org Michael Nikolaou approach for approximating transcendental transfer functions by polynomial rational fractions in process control. 11 Figure 9 . Reduction of the infectious fraction from (0) to 0 (a) in a single front in terms of the Heaviside step function with = 1/ , (b) following a first-order Padé approximation, and (c) following a second-order Padé approximation. (cf. Figure 2) . Because of the critical role of Padé approximation in deriving eqn. (14) for the ODEs of the modified SIR model to approximate the dSIR model, we will use the term Padé SIR to denote the modified SIR model structure. Finally, it should be noted that eqn. (13) may seem counter-intuitive, as it appears to suggest that the generation and depletion rates of are 2 and −2 / , respectively. However, eqn. (13) rather suggests that while the susceptible depletion rate remains − , the infectious depletion rate appears as ( − 2/ ) (which is -dependent) rather than − / (which is not) due simply to replacement of the delay term − in the Laplace domain by its Padé approximant. Similarly, the rate of increase for in eqn. (13) is -dependent, rather than not. Standard theory of DDEs (e.g. equation (9.42) in Ch. 4 and equation (1.2) in Ch. 5 of Kuang, 12 or Gopalsamy 13 ) can be applied to establish rigorous properties for the dSIR model, such as global stability, convergence to a final steady state, and others. A complete analysis is beyond the scope of this paper. However, some important theoretical properties of the dSIR model structure of practical interest are discussed next, particularly in comparison to their SIR counterparts. The analysis of the Padé SIR model follows standard ODE analysis and is presented more briefly, except when it has important implications for either theoretical or practical issues. Eqn. (7) can be used to show (Appendix B) that an equilibrium point ̅ is stable and an epidemic outbreak does not occur iff Note that the stability upper bound 1/( ) for the dSIR model in the above eqn. (17) coincides with the well known bound / = 1/ 0 dictated by the SIR model under the widely used correspondence The same stability bound can be derived for the Padé SIR model structure using standard ODE analysis based on linearization around an equilibrium point. Accepting for now without proof the global stability of eqn. where 0 ≝ and is the Lambert function 14,15 of order 0. The standard plot for ∞ , equal to the total fraction of infected by the end of an epidemic, 2 is shown in Figure 10 for completeness. Note that the plot is valid for both SIR and dSIR models, with 0 = / and 0 = 1/ , respectively. Interestingly, while use of the Lambert function to solve problems such as the above was pointed out as early as 1996, 14 it may have escaped the attention of most literature in this field. 2 The Lambert function in its various forms will turn up in a number of results below. It should also be noted that eqns. (19) and (20) are the same for the SIR and Padé SIR model (Appendix C) under the correspondence between and in eqn. (18) . For the initial part of a spreading epidemic, starting from a perturbation of the steady state ( , , ) = ( ̅ , 0,1 − ̅ ) as ( < 0) = 0, (0) = , ̅ ≝ ̅ 0 > 1 it can be shown (Appendix D) that the infectious fraction initially grows approximately at an exponential rate as . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. where Δ ≝ − ̅ and the constants , are in terms of 0 ≝ (Appendix D) with as shown in Figure 11 . Similarly, the initial exponential rate of the Padé SIR model is given by where 0 = / , and for the standard SIR model, the exponential rate is where 0 = / , shown in Figure 11 as well. This has immediate implications for the early rate of rise of the infectious fraction to its peak, * , as illustrated in Figure 12 . Note that while the initial SIR rate is half of the Padé SIR rate and about half of the dSIR rate, all three models eventually reach the exact same steady-state values, as captured by eqns. (19) and (20). While it is not obvious to the author whether the peak * can be easily obtained for the dSIR model, a good approximation can be obtained (Appendix E) through the Padé SIR model, as * = 1 , * = 2 � −1 ln( ) − 1 + 1� (25) Note that the above * , exact for the first-order Padé SIR model and approximate for the second-order Padé SIR and dSIR models, is double the * of the standard SIR model, as confirmed in Figure 8 . Once more, there are obvious practical implications from this discrepancy. Note also that in case an upper bound is placed on * , to avoid overwhelming hospitalization facilities during an epidemic, eqn. (25) has an explicit solution for the corresponding maximum 0 ≝ as For comparison, the standard SIR model yields The values of 0 indicated by eqns. (26) and (27), with corresponding definitions, are shown in Figure 13 . It is evident that the Padé SIR model places twice as tight a restriction on 0 as the standard SIR model, if is not to exceed the * value. The dSIR model structure developed in Figure 3 can be easily extended to include additional compartments. In fact, practically all population models of infections developed to date using the concept of exchange between compartments 3 can be immediately translated (a) from ODEs to DDEs through replacement of compartment drain rates proportional to the drained quantity by drainage of amounts that have resided for a certain time, , in that compartment, or (b) from ODEs to Padé approximations that maintain the ODE structure but are more realistic. We illustrate these ideas on . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. . https://doi.org/10.1101/2020.04.26.20080960 doi: medRxiv preprint medRχiv.org Michael Nikolaou the SPIR model structure, in view of its importance for the recent coronavirus epidemic. [16] [17] [18] A distinct feature of the coronavirus causing COVID-19 is that it enables infection transmission at the pre-symptomatic stage. 19 Therefore, the SPIR model structure (Figure 14 ) comprises the usual population fractions { , , } along with the pre-symptomatic infectious fraction, , with symptomatic infectious being . SPIR differs from the standard SEIR structure 2,4,5,20 by the way the four compartments interact. While it is formidable to practically monitor infections on presymptomatic infectious individuals, monitoring symptomatic infectious is more reasonable, as symptoms are clear and can be confirmed by testing. Therefore, tighter restrictions can be placed on the I group, in addition to the P group typically following general restrictions placed on the general population to curb the spread of the epidemic, as captured by the spread factors and , respectively, in Figure 14 . The dynamics of the dSPIR structure is shown in Figure 15 , which follows the pattern of Figure 3 , with the addition of the pre-symptomatic infectious fraction . Following the same technique as in section 2, we can immediately write the following equations for the dSPIR model structure by visual inspection of Figure 15 : Combining the above equations and letting → ∞ yields the nonlinear DDE and associated delay algebraic equations where and − are the durations of an individual's staying in the P and I compartment, respectively, and typically Note that the form of the SPIR model following the standard SIR pattern, eqns. (1)- (3), is the coupled ODEs 21 In the form of these ODEs is the following more realistic Padé SPIR model, which approximates the dSPIR DDEs, eqns. (31)-(35), significantly better than the SIR ODEs (Appendix F): . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. The illustration and analysis presented in sections 3 and 4 can be easily repeated for the dSPIR and Padé SPIR structures. To maintain the scope of this publication, only a few properties will be explored below. The rest will be explored in more detail in forthcoming publications. We present here only a few simulations comparing the { , , , } profiles resulting from numerical solution of the SPIR, dSPIR, and Padé SPIR models. The values are used in all simulations. 8,9 Figure 16 illustrates the differences in the infectious peaks, * and * , similar to these in Figure 8 . Following the same approach as in section 4.1, eqn. (31) can be used to show (Appendix G) that an epidemic outbreak does not occur around an equilibrium point ̅ iff The proof in Appendix G is by approximation. It is conjectured that the bound in the above equation is exact. This conjecture will be examined in subsequent studies. It can be shown (Appendix H) that the dSPIR counterpart of eqn. (19) is where > and The counterpart of the SIR and dSIR plot ( Figure 10 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. A subtle inconsistency in the standard SIR model structure was pointed out. This inconsistency arises from misinterpretation of an assumption explicitly articulated in the original publication of Kermack and McKendrick: 1 that the depletion of the infectious compartment is proportional to the content of that compartment. The depletion rate constant, , is usually interpreted as equal to the inverse of the residence time in the infectious compartment, namely the duration of the infection, , for each individual (eqn. (18) ). To the extent that this duration is about constant, the analysis presented here suggests that the preceding interpretation is incorrect, leading to certain erroneous conclusions. A corresponding model structure, termed dSIR, was developed, to account for each individual leaving the infectious compartment after a certain duration. The dSIR model structure comprises a single DDE for the susceptible fraction, s, and associated algebraic equations capturing the dependence of the remaining population fractions on . While both SIR and dSIR produce the same results for assessment of stability and final values, the SIR model produces a maximum of the infectious fraction, * , about half of its dSIR counterpart. This has profound consequences if the SIR model is used to predict * during an epidemic. It is also noted that even if the SIR model parameters and are estimated based on experimental data fit -albeit under the wrong interpretationmodel predictions are still going to be inaccurate. This is because the standard SIR model, comprising the three ODEs in eqns. (1)-(3), is structurally different from the DDE form of the dSIR model, eqns. (7)-(9), or even from the ODE form of the Padé-based approximation of the dSIR model DDEs, eqns. (1) and (13) or (1) and (14) . The dSIR structure can be easily extented to other compartment-based population models. 3 Such an extension to the dSPIR model structure was presented and briefly illustrated and analyzed. This model is important for infections transmitted by both pre-symptomatic and symptomatic infected individuals. Similarities and differences between SPIR and dSPIR models are of the same nature as between SIR and dSIR models. Numerous additional issues related to this work can be considered, including the following: Rigorous analysis of DDE models; DDE modeling and analysis for a distribution rather than uniform delay (cf. Figure 2) ; resolution of conjectures presented in the text; and implications for different forms of infection kinetics. Such issues will be addressed in forthcoming publications. All computations were done in Mathematica, available at the University of Houston. Sharing of teaching material about the SIR model on Github by Prof. Jeff Kantor of Notre Dame is also gratefully acknowledged. Research reported in this publication was partially supported by the Institute of Allergy and Infectious Diseases of the National Institutes of Health under award number R01AI140287, financed with Federal money. The content is solely the responsibility of the author and does not necessarily represent the official views of the National Institutes of Health. The roots of the transcendental characteristic equation is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. For stability, all must be in the left-half of the complex plane. We show first that for any positive value of ̅ , no complex root = + can cross the imaginary axis to move from stability to instability. Because, if = were a root of eqn. (58) for ̅ > 0, it would be Therefore, only real roots should be considered in eqn. (59) for stability analysis. Furthermore, since the relevant values of are 0, −1 in eqn. (59) for ∈ ℝ, and < 0 implies which is satisfied for < 1 (see figure below) leading immediately to eqn. (17) . Appendix C. Proof of eqn. (19) Dividing eqn. (7) by yields Continuing on the last equation we get which is eqn. (19) . The Padé SIR model also reaches the same result: Dividing eqn. (1) by eqn. (13) and rearranging yields Taking the limit as → ∞ with 1/ = , 0 ≈ 1, and 0 ≈ 0 yields eqn. (19) . Proof that the standard SIR model also reaches the same result follows the same pattern and is omitted for brevity. Therefore, the terms rapidly decay for ≥ 1, and the summation in eqn. (65) quickly becomes approximately equal to . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 1, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 1, 2020. . https://doi.org/10.1101/2020.04.26.20080960 doi: medRxiv preprint A Contribution to the Mathematical Theory of Epidemics Modeling Infectious Diseases in Humans and Animals Population biology of infectious diseases: Part I Mathematical Biology: I. An Introduction Infectious diseases of humans: dynamics and control Evolving Epidemiology and Impact of Non-pharmaceutical Interventions on the Outbreak of Coronavirus Disease Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application Risk Assessment of Novel Coronavirus COVID-19 Outbreaks Outside China An interactive web-based dashboard to track COVID-19 in real time Chemical Process Control: An Introduction to Theory and Practice Delay Differential Equations with Applications in Population Dynamics Stability and Oscillations in Delay Differential Equations of Population Dynamics On the LambertW function The Lambert function belongs in the Function Hall of Fame. Chemical Engineering Education In preparation Report 9 -Impact of nonpharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand Institute for Health Metrics and Evaluation (IHME)