Ziegler and Nichols meet Kermack and McKendrick: Parsimony in dynamic models for epidemiology.

The COVID-19 crisis popularized the importance of mathematical modeling for managing epidemics. A celebrated pertinent model was developed by Kermack and McKendrick about a century ago. A simplified version of that model has long been used and became widely popular recently, even though it has limitations that its originators had clearly articulated and warned against. A basic limitation is that it unrealistically assumes zero time to recovery for most infected individuals, thus underpredicting the peak of infectious individuals in an epidemic by a factor of as much as about 2. One could avoid this limitation by returning to the original comprehensive model, at the cost of higher complexity. To remedy that, we blend Ziegler-Nichols modeling ideas, developed for automatic controller tuning, with Kermack-McKendrick ideas to develop novel model structures that predict infectious peaks accurately yet retain simplicity. We illustrate these model structures with computer simulations on real epidemiological data.

Introduction

The COVID-19 crisis popularized the importance of mathematical modeling for managing infectious disease epidemics (Adam 2020; Giordano et al., 2020; Jewell et al., 2020; Kucharski et al., 2020; Wang et al., 2020). Concepts such as herd immunity, or flattening the curve entered the global vernacular shortly after the viral epidemic started its ominous exponential spread (Tufekci 2020). Underlying these concepts is seminal modeling work originating in the 1920s, most notably by Kermack and McKendrick (1927) (K-M) who contributed a remarkably general model to the mathematical theory of epidemics. That model partitions a fixed-size population into compartments S (susceptible to infection), I (infectious as a result of infection), and R (the remaining fraction of the population, comprising individuals removed from the I compartment by recovery or death or already being immune to the disease from the outset1 ) and keeps track of the corresponding population fractions, over time by balancing loading and discharge (removal) rates for each compartment (Fig. 1 ). The result is a set of three general nonlinear integrodifferential equations (IDEs, APPENDIX A).

Fig. 1

The SIR model structure. Note that

In developing their general IDEs, K-M monitor at time the fraction of individuals in compartment I who have been infectious for a time period since their infection at time , and express their removal rate from I as They further simplify their IDEs by assuming that , which they combine with constant infectivity, , to get the celebrated SIR model (Anderson and May 1979; Anderson et al., 1992; Diekmann et al., 1995; Murray 2002; Keeling and Rohani 2008; Brauer and Castillo-Chavez 2012; Brauer 2017), comprising the equations where are constants (counterparts of the K-M parameters ) related to infectivity and removal rates, and the ordinary differential equation (ODE) may trivially be used in place of Eq. (3). The SIR Eqs. (1)-(3), as well as the general IDEs from which they were derived, lead K-M to state two important principles about (a) how epidemics spread (“No epidemic can occur if the [susceptible] population density is below this threshold value.”) and (b) how epidemics end (“… from a particular relation between the population density, and the infectivity, recovery, and death rates.”). In fact, K-M showed that the corresponding values of both the threshold value for epidemic spread and the long-term values of the three fractions, , of a population that just went through an epidemic are robust, namely not affected by the particular form of the kernels (infectious duration distributions) involved in the IDEs referred to above. The simplicity, intuitive appeal, and descriptive power of the SIR model made it enormously successful over the years, either as a stand-alone or as a module in numerous combinations and permutations of basic compartments and associated dynamics, as discussed by Hethcote (1994) who figuratively refers to about an order of magnitude below a myriad of epidemic models.

Recent use of the simple SIR model to guide the management of COVID-19 has brought to the forefront the importance of the above two K-M principles through two basic insights, with profound implications for designing interventions against this epidemic (Nikolaou 2021, and references therein): (a) to contain the spread of COVID-19 without any intervention, the susceptible fraction of the population would have to drop below 40%, and (b) to achieve this level of immunity in the population (above 60%) by contacting the disease (in the absence of a vaccine) would have meant that about 90% of the population would have to be infected throughout the epidemic. With mortality rate of the infected at about 0.5%, the immunity building option just described would have meant fatalities for about 0.5% of the population, an inordinately large number. This vital insight was provided by a simple model that is strictly speaking “wrong” yet utterly useful (Box 1979).

As important as the standard SIR model is, there is also a subtle shortcoming in Eqs. (1)-(3) that had not been fully appreciated before (Kemper 1980) but emerges naturally when one considers “Flattening the Curve” (Qualls et al., 2017; Ferguson et al., 2020)): As the infectious fraction, goes through a peak value to subsequently decline asymptotically to zero, is not accurately characterized by the SIR model. The values of produced by the SIR and by the full IDEs may differ by as much as a factor of about 2, as pointed out by Nikolaou (2021), the essence of whose arguments are summarized below (Section 2.2). In the same fashion, the rise time to suggested by the SIR model can also be commensurately inaccurate (Nikolaou 2021). These shortcomings are of conceptual and practical importance when designing epidemic management strategies to lower the infectious fraction peak, , and thus lessen the burden on resources for treatment of infected patients in need of intensive care. More importantly, the issue extends to the numerous compartment-based model variants spawned by the SIR model for management of various kinds of infectious disease epidemics (Hethcote 1994) and even concerns efforts to capture the dynamics of epidemics by sophisticated model structures that may embed SIR modules (Horrocks and Bauch 2020).

It turns out that these shortcomings emerge from uncritically following the assumption of constant (in the original K-M notation), which intuitively corresponds to an exponential distribution of infectious duration (Fig. 2 ). An exponential distribution suggests that most infectious individuals are removed from I in zero time after infection. This is a realistically untenable proposition, which is not merely inaccurate but has important consequences, as we discuss next.

Fig. 2

Sample cumulative distribution functions (top) and corresponding probability distribution functions (bottom) for discharge time from the I compartment of a population, ( in dimensionless form). Curves follow the formulas and (see APPENDIX B). The exponential distribution with (top) and (bottom) corresponds to , whereas the impulse distribution with (unit step at , bottom) and (unit impulse at , top) correspond to

It can be easily demonstrated (Nikolaou 2021) that as the infectious duration probability density function (PDF) peaks away from (corresponding to the cumulative density function (CDF) approaching a step, Fig. 2) the response of the original system of IDEs quickly approaches that of a system of delay differential equations (DDEs) as shown in Fig. 3 . For the latter, Nikolaou (2021) developed a simple counterpart of the standard SIR model using Padé approximations for replacement of Eq. (2) of the SIR model by an equally simple ODE.

Fig. 3

Response of the infectious fraction, , according to the model of Eqs. (1), (4), and (3) for distributions shown in Fig. 2. Note that the distribution for is closer to the exponential distribution than to the impulse distribution in Fig. 2, yet the response of for is a lot closer to the response for rather than to that for

Nevertheless, situations arise (Wallinga and Lipsitch 2007) where the infectious duration CDF, albeit far from exponential, is still not close enough to a step and, as a result, DDEs (and, accordingly, their Padé approximations by ODEs mentioned above) may not be a reasonably accurate approximation of the general IDEs. The purpose of this publication is to develop a simple approximation for such a case. Specifically, we borrow from the field of automatic control, where a high-order system may be fruitfully approximated by a first-order-plus-time-delay (FOPTD) model. Ziegler and Nichols (Z-N) showed that the FOPTD model is convenient both to obtain from data and to use in PID controller tuning following the celebrated tuning rules developed by these authors (Ziegler and Nichols 1993). Applying the Z-N idea to the K-M IDEs, we develop a continuum of simple models which can accurately approximate the original general IDEs by simple ODEs that retain the simple SIR model structure and feature only one additional parameter easily obtainable from epidemiological data on infectious duration CDFs.

In the rest of the paper, we first provide a summary of relevant background facts and an overview of subsequent developments. We then show our main results. We use a case study based on relevant data easily available in literature to illustrate these results. Finally, we present a broader context, conclusions, and recommendations for further study.

Background and overview

The basic model for general infectious duration distributions

For a fixed-size population with initially susceptible, infectious, and immune fractions respectively, one can stack the values of and monitor their changes over discrete time by performing corresponding mass balances, as shown in Fig. 4 Nikolaou 2021). For infectious duration CDF, , remaining invariant throughout an epidemic, the model resulting from Fig. 4 comprises Eqs. (1), ((3), and the IDE(Cushing 1977; Hethcote 2000).

Fig. 4

Schematic of time-varying susceptible ( green), infectious ( orange), and removed (blue) fractions of a fixed-size population after an initial infection, at discretized time (Nikolaou 2021). Each new part of the infectious fraction (thick-framed orange rectangles) moves to the removed fraction, (thick-framed blue rectangles) piecewise in a number of time steps, following a certain distribution of infectious duration (age). The population eventually reaches a steady state at and (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Models for specific infectious duration distributions

For in Eq. (4), combination with Eqs. (1), (3) immediately yields the SIR model Eqs. (1)-(3) withrepresenting the average infectious duration and or, equivalently, representing dimensionless time. Similarly, for CDF being the delayed Heaviside (step) function as shown in Fig. 2, one can immediately get the algebraic delay equation (ADE) (Hethcote 2000)

which, along with Eqs. (1) and (3), constitutes a full dynamic model for

It was recently shown, using first- or second-order Padé approximations in the Laplace domain Nikolaou 2021), that the d–SIR model, comprising the above Eqs. (6), ((1), and (3), can be well approximated by ODEs, namely Eqs. (1), (3), and either of the following two ODEs:

orIn addition to the exponential or step CDF discussed above, a variety of physically meaningful CDFs between these two extremes have been considered for (Byrne et al., 2020). While these intermediate CDFs have a few minor differences, they all follow simple profiles, as shown in Fig. 2. It is for CDF profiles qualitatively following such patterns that we show in the next section how to create simple approximations using Z-N ideas.

In fact, we show that the entire span of IDE models resulting from CDFs ranging from exponential to step (e.g. as shown in Fig. 2) can be approximated by a family of models similar, in form, to the SIR model Eqs. (1)-((3)) with the simple addition of a corresponding time delay to Eq. (2). The delay renders the counterpart of Eq. (2) a DDE.

We further show that approximation can remarkably simplify that DDE to the simple ODEwhere the parameter takes values in the interval and approximately produces the range of responses for shown in Fig. 3. More specifically, at the lower end of that interval Eq. (9) trivially yields Eq. (2) of the standard SIR model (corresponding to exponential CDF) whereas at the upper end of that interval Eq. (9) trivially yields the first-order Padé-SIR model comprising Eqs. (1), (7), and (3). Behavior between these two extremes can be captured by knowledge of this single parameter.

It should be stressed, as mentioned in Introduction, that the threshold value for epidemic spread and long-term values suggested by all models discussed coincide. It is in transients where significant differences emerge.

The basic reproductive ratio

Making time dimensionless in the above models, asor, equivalently, as trivially results in retention of only one parameter in the corresponding equations, namely . That parameter, implicitly suggested in the K-M paper and explicitly introduced as such later MacDonald 1957; Heesterbeek 2002; Brauer 2017) is called “the basic reproductive ratio” and is widely considered “one of the most critical epidemiological parameters” (Dietz 1993; Keeling and Rohani 2008). In fact, making time dimensionless in Eqs. (1) and ((2) immediately yields This suggests that the dynamics of the SIR model is essentially governed by a single parameter, (hence this parameter's paramount importance, as already mentioned). This observation corroborates the use of dimensionless time in Fig. 2 and Fig. 3, a convention that will also be followed in the sequel.

The value of depends on many factors, both natural and resulting from human interventions, and is typically estimated from epidemiological data (Dietz 1993). Nikolaou (2021) presented a detailed analysis of the systematic errors that emerge in predictions by the SIR model if is estimated from data collected in the early exponential growth of an epidemic and is combined with estimates of the average infectious duration, . It is also noted that numerous refinements of may be considered in more elaborate models (the proverbial “alphabet soup” of models such as SEIR, SIR, and many others) that account for factors such as (a) compartments and corresponding interactions different from the standard S, I, and R (Hethcote 1994), (b) distinct subpopulations based on age, social contact structure, or other clustering factors (Anderson and May 1991; Anderson et al., 1992; Hethcote 1996; Ferguson et al., 2006; Keeling and Rohani 2008)), (c) variation in space (Keeling and Rohani 2008)), or (d) combinations thereof, which collectively lead to diverse stratification patterns (Kestenbaum 2019). For all such models, the term capturing discharge from a corresponding compartment typically corresponds to an exponential CDF of time from loading to discharge (even in advanced modeling exercises purporting to develop sophisticated tools such as automated algorithmic discovery (Horrocks and Bauch 2020)) with all limitations already elaborated on.

In the ensuing discussion we will consider that takes values above 1, because the “Threshold Theorem” of K-M, referred to in Introduction, asserts that an epidemic spreads iff

Main results

The main idea in the subsequent developments, previewed in summary in Section 2.2, is that the family of CDFs shown in Fig. 2 can be thought of, in approximation, as the unit-step response of a FOPTD model – the celebrated Z-N idea developed for control systems. We show next how this idea can be exploited to yield a dynamic model for that approximates Eqs. (1), (4), and (3), and how further simplification by subsequent approximations may arise to span a broad spectrum of dynamics entailed by the K-M IDEs.

Z-N for K-M

Taking Laplace transforms in Eq. (4) withand using the convolution theorem yieldswhere is the Laplace domain variable in place of the most commonly used (reserved here for the susceptible fraction of the population).

Because and can be thought of as a unit step response, as mentioned above, one can write the simple approximation(Fig. 5 ) which, combined with Eq. (14), yieldswhere .

Fig. 5

Approximation of by the CDF , the response of a system comprising a first-order element (of steady-state gain 1 and time constant ) followed by a time delay () element to a unit-step-change on the system's input at time . Note that the average corresponding to the above CDF is the area between the CDF and 1 for , which can be shown to be (APPENDIX C).

Therefore, Eqs. (13) and (16) imply that the dynamics of is captured by Eqs. (1), (3), and the simple DDE

The FOPTD-SIR model

Substituting by in Eq. (17) and using Eq. (1) immediately yieldsThe above Eq. (18) generalizes Eq. (2) of the standard SIR model in a direct and simple manner, through replacement of the term (corresponding only to an exponential distribution of infectious duration) by which can approximately account for a variety of infectious duration distributions. The resulting FOPTD-SIR model comprises Eqs. (1), (18), and (3). Note that Eq. (2) is trivially recovered from Eq. (18) for Similar to Eqs. (11) and (12), which render the standard SIR model dimensionless, thus linking its dynamics to alone, the FOPTD-SIR model can be made dimensionless as well. This can be done by considering the natural counterpart of the standard SIR model parameter (inverse of the average of infectious duration when the latter follows an exponential distribution). This counterpart is the average, of infectious duration when the latter follows the CDF shown in Fig. 5 can be easily shown (APPENDIX C) to bewhere and are the dead time and time constant of the first-order lag, respectively. Making time dimensionless in Eqs. (1) and (18), and defining asimmediately yields whereandThe above suggests that the dynamics of the FOPTD-SIR model is essentially governed by a single additional parameter, , compared to of standard SIR dynamics. This is a remarkable simplification of dynamics representation by approximation, compared, for example, to representations based on full account of infectious duration distributions (Wallinga and Lipsitch 2007).

Properties of the FOPTD-SIR model

It is worth exploring the form of the FOPTD-SIR model for different values of , given infectious duration distribution profiles such as shown in Fig. 2.

An exponential distribution of infectious duration ( in Fig. 2) corresponds to and , yielding . This is the lowest possible value of , corresponding to the standard SIR model, as immediately evident from Eqs. (21) and (22), which revert to Eqs. (11) and (12), respectively.

An impulse distribution of infectious duration (with unit step CDF, in Fig. 2) corresponds to and , yielding which by simple manipulation turns Eqs. (22) to the algebraic delay equationthus resulting in a d-SIR model comprising Eqs. (11), (25), and (3). Incidentally, the above Eq. (25) can also be directly deduced from the graph in Fig. 4, with the infectious duration CDF being a step rather than S-shaped. In fact, we will make the case in section 3.6 that the behaviors of FOPTD-SIR models for virtually coincide, thus leaving the as the range of relevant values for .

The properties of the FOPTD-SIR model can be studied using methods for DDEs well studied in biology (Gopalsamy 1992; Kuang 1993)). Even though DDEs – possessing infinite spectra and derivative discontinuities, among other features – may be considered more difficult to analyze than ODEs (Cushing 1977, p. 5) the corresponding theory “does not present substantial additional difficulties” (Bellen and Zennaro 2003, p. 6). Therefore, we will only present next two important properties of the FOPTD-SIR model, to confirm that this model is in agreement with the standard SIR model: The corresponding value for herd immunity suggested by the Threshold Theorem, in section 3.4; and the long-term values of after an epidemic, which capture the total number of infections throughout an epidemic, hence the number of potential fatalities. To establish a difference between FOPTD-SIR and standard SIR, an analytical expression for characterization of the maximum of the infectious fraction along with the exponential growth rate leading to it will be developed through introduction of a further simplified SIR-like model, in section 3.6.

The threshold theorem and herd immunity for the FOPTD-SIR model

Applying standard stability theory for the FOPTD-SIR model linearized around the stationary point using Taylor series yields (APPENDIX D) that the threshold value of , above which the epidemic is guaranteed to spread (i.e. the population does not possess herd immunity), isas anticipated by general K-M results.

Long-term values of for the FOPTD-SIR model

It can be shown (APPENDIX E) that at the end of an epidemic that started at and the total fraction of infected throughout the epidemic isaccording to the FOPTD model, where is the Lambert function (Kesisoglou et al., 2021). These values are the same as produced by the standard SIR model (Keeling and Rohani 2008).

The model

While Eq. (18) is a useful generalization of Eq. (2) (in that it accounts for a variety of infectious duration distributions other than exponential), it is in the form of a DDE, which, as we discussed, might be less appealing than an ODE for analytical computations. To get the analytical computation benefits of a simple ODE from Eq. (18), assuming no significant compromise of accuracy, one can show (APPENDIX G) that the following approximation of Eq. (18) can be used, which immediately turns out to be Eq. (9) presented in dimensionless form:According to its definition in Eq. (23), the parameter in the above Eq. (28) takes values in , as the CDFs shown in Fig. 2 traverse the range from exponential to unit step . Fig. 3 makes the case that the resulting dynamics for a contiguous family of these CDFs, as the index increases above a certain value towards , is tightly clustered very close to the dynamics resulting from the CDF corresponding to , i.e. to the dynamics of the d-SIR model presented above. The d-SIR model has been analyzed in detail by Nikolaou (2021), who also introduced an approximate ODE version of the d-SIR model using first- and second-order Padé approximations, eventually leading to Eqs. (7) or (8), respectively. Simple inspection of the above Eqs. (28) and (7) immediately indicates that Eq. (28) for is exactly Eq. (7) (in dimensionless time). Therefore, the relevant range of values for the parameter , over which appreciable differences can be seen in the dynamics of respective models, isThe corresponding CDFs for FOPTD-SIR models in dimensionless time are shown in Fig. 6 . This figure indicates that all dynamics resulting from dimensionless dead time ( by Eq. (24)) remain essentially indistinguishable, thus leading to Eq. (29); dimensionless dead times result in a range of distinguishable dynamics.

Fig. 6

Family of infectious duration CDFs of the form

for or (gray area)

The dynamics of FOPTD-SIR or -SIR models resulting from CDFs with (hatched area) are similar.

Eq. (28), combined with Eqs. (1)/(11), and (3), offers an exceptionally simple parametrization of an SIR-like model, which we call -SIR, for brevity. That model approximates a wide range of epidemic spread dynamics, without necessarily involving the specific CDF of infectious duration, resorting instead to its FOPTD approximation.

Eq. (28) trivially suggests that the threshold value for resulting from the -SIR model is exactly the same as produced by the standard SIR model, for the same average of the infectious duration distribution. It is also easy to show (APPENDIX E) that the long-term values of resulting from the -SIR model are exactly the same as those produced by the standard SIR model. These findings confirm that -SIR model provides the same insight as the standard SIR model regarding (a) herd immunity, and (b) projections of total number of infected throughout an epidemic and consequent number of deaths. However the -SIR model (therefore, by close approximation, the FOPTD model as well) provide better insight on how to “Flatten the Curve”, by providing more relevant representation of the infectious fraction peak, , and the initial exponential rate of rise of towards . This is discussed next.

Peak of infectious fraction

It can be easily shown (Kermack and McKendrick 1927) that the peak value of the infectious fraction in the standard SIR model is . A similar expression cannot be derived for the peak value suggested by the FOPTD-SIR model. However, the -SIR model can be used to derive the approximate expression(APPENDIX F) . The above approximation suggests that as increases, the infectious fraction peak, of the FOPTD-SIR model can be as high as about twice the infectious fraction peak, of the SIR model (Nikolaou 2021) as demonstrated in Fig. 3.

Exponential growth rate of infectious fraction

All models discussed in this manuscript produce approximately exponential growth of the infectious fraction at early stages of a spreading epidemic. Expressions for these initial growth rates can be easily written following Taylor-series linearization of each model around a steady state (equilibrium point) comprising and . Taking the -SIR model with to represent the entire range of models considered here (from standard SIR for to d-SIR for ), linearizing Eq. (28) around the equilibrium point immediately yieldsfor and . The above Eq. (31) implies that, for a certain , the standard SIR model, with a growth rate , may underpredict the initial growth rate of the infectious fraction of the FOPTD-SIR model by as much as a factor of 2.

It is also of practical significance to examine the values of predicted by the entire range of -SIR models when values of are estimated from available data on early exponential growth of the infectious fraction. (The average infectious duration can be estimated relatively precisely from epidemiological data (Keeling and Rohani 2008, section 2.1.1), thus time can be easily made dimensionless.) Eq. (31)indicates that the corresponding dimensionless-time growth rate isSolving Eq. (32) for and substituting into Eq. (30) for yieldsThe ratio of over of the standard SIR model is shown in Fig. 7 . That figure indicates that as moves from 1 towards 2 discrepancies between and increase, except on a single line where the surface intersects the plane at 1.

Fig. 7

Ratio of the maximum infectious fraction produced by the -SIR model, , over the maximum infectious fraction produced by the standard SIR model, , each tuning its corresponding value of based on the same rate of exponential growth of the infectious fraction.

Applicability range of approximations

While all approximations considered above invariably include approximation errors, they must also conform to the fundamental requirements , and . We examine below the implications of these inequalities for the -SIR model.

The inequality is trivially satisfied by Eq. (1)/(11), as long as .

Satisfying the inequality is more interesting, as Eqs. (28), (1)/(11), and (3), imply , which requireswith . For the above inequality is guaranteed for which covers a range of infection spreads of practical significance, particularly if containment measures are taken. Nevertheless, in many practical situations as well. This would lead to for initial times , which is clearly infeasible. That situation can be handled by the simple -SIR model adjustmentand combination of the above Eq. (35) with Eqs. (1)/(11), and (3). The resulting dynamics will quickly revert to the original -SIR model dynamics as decreases and makes in Eq. (35).

The inequalities can be handled in a similar manner. For example, Eq. (30) suggests that .

A full investigation of the limits of the above investigations (e.g. for very large ) is beyond the scope of the current analysis, and will be pursued elsewhere. Nevertheless, the analysis presented above already covers important aspects of the proposed approximations. Additionally, the case study presented next indicates that the -SIR model works remarkably well for a wide range of values of

Case study

We use literature data analyzed by Anderson et al. (1986) to illustrate the concepts developed in the previous sections. The analysis presented therein is based on use of the standard SIR model for predictions of corresponding compartment sizes, even though data on a corresponding CDF is also shown. The data refer to the spread of HIV infection and development of AIDS. Specifically, the data shown in Fig. 8 refer to the infectious duration CDF governing the incubation period, namely the time from infection of previously susceptible individuals (entering the I compartment) to appearance of full-blown AIDS symptoms (individuals leaving the I and entering the R compartment). Note that entering the R compartment for this case refers to mere removal from the I compartment rather than recovery or death. The latter issue, as gravely important as it is, is not the focus here; if of interest, that issue can be modelled by augmenting this SIR model (Keeling and Rohani 2008, section 2.2); for a potential structure in series with SIR see Fig. 5 in Lipsitch et al. (2003). The data set of this case study is deliberately selected to correspond to an infectious duration CDF that is not close to a step, so that the approach developed here can be meaningfully applied.

Fig. 8

Fitting of HIV epidemiological data (Anderson et al., 1986) by the CDFs

with (left),

with (right). 95% confidence bands for mean predictions are shown. All calculations and plotting performed on Mathematica (Wolfram Research Inc. 2021).

Excellent data fit by the Gamma and FOPTD CDFs is shown in Fig. 8. It is noted that for the Gamma CDF Fig. 2), suggesting that the neither the SIR model nor the Padé SIR model is the most appropriate for this case, according to Fig. 3. Rather, the full IDE model (Eqs. (1), ((4), (3) featuring the data-fitted Gamma distribution in Fig. 8) or the corresponding FOPTD-SIR and the -SIR counterpart are more appropriate. We illustrate these claims in the following simulation results for different values of :

The SIR model, Eqs. (1)-(3), with .

The Γ–SIR model, Eqs. (1), (4), (3), with .

The FOPTD–SIR model, Eqs. (1), (18), (3), with .

The –SIR model, Eqs. (1), (28)/(9), (3), with

The values of selected correspond to those used in Anderson et al. (1986) for rapid spread of the disease or somewhat tamed Because Eq. (35) is used in all simulations.

A number of interesting patterns are evident in Fig. 9:

Fig. 9

Simulation results for cases (a)-(d) corresponding to different values of in Case Study.

There is obvious discrepancy between the profiles produced by the standard SIR model and the rigorous Γ–SIR model. This discrepancy is most important for the peaks of the infectious fraction, indicated by each of the two respective models, with the peaks rigorously produced by the Γ–SIR model being about higher than their SIR counterparts. This is in agreement with Eq. (30) which suggests would be about higher than for . The significance of this discrepancy for managing epidemics was already emphasized in Introduction.

There is remarkable agreement between the infectious fraction profiles produced by the rigorous Γ–SIR model and the FOPTD-SIR model. This agreement underscores the value of the FOPTD-SIR model, which accurately captures model dynamics by mere addition of one easily interpretable parameter to the standard SIR model.

Although somewhat less accurately matching the profiles of rigorously produced by the FOPTD–SIR model, the –SIR model is still clearly superior to the standard SIR model, while retaining simplicity. This advantage is afforded by the –SIR model with only one additional parameter, , whose value in the interval qualitatively communicates the behavior of the –SIR model between two extremes: The standard SIR model and the d–SIR model .

Conclusions and discussion

Combining ideas from two seminal contributions in their respective fields of epidemiology (K-M) and automatic control (Z-N), we developed simple model structures that capture the basic dynamics of epidemics described by general, rigorous integrodifferential equations. The merit of the proposed modeling approach was illustrated via computer simulations on actual epidemiological data from the literature.

While the results presented were developed for a three-compartment, fixed-size population, they bear relevance to the broader class of compartment-based models that can capture the dynamics of epidemics for any infectious agent. Extensions, study of their properties, and testing on actual epidemiological data would certainly be of interest for future investigations.

The models proposed here are parsimonious, in that they retain a simple structure with the addition of only one parameter, which accepts intuitive interpretation. As such, the proposed models are easy to include in related software for practitioners or researchers (Obadia et al., 2012).

Of course, the proposed models are, by design, empirical, aiming at usefulness rather than ultimate accuracy of representation. Indeed, the message about the utility of models of such nature was emphasized in both publications from which we drew here: After they use their SIR model to fit data from a plague outbreak, Kermack and McKendrick (1927) state that

deductions as to the actual values of the various constants should not be drawn. It may be said, however, that the calculated curve, …, conforms roughly to the observed figures.

Similarly, in developing a model useful for controller design, Ziegler and Nichols (1993) observe that

The difficulty of dealing mathematically with processes involving a series of lags … is very great indeed. An approximate description of the characteristics of a process is given by values of … two quantities. … True, these two are only a rough measure of the entire reaction curve, … but they give enough of the story to allow a prediction … .

There is a similar message here, which resonates with the well known dictum permeating all science and engineering: “All models are wrong, but some are useful” (Box 1979). We hope that the models proposed here will be useful.

Author statement

None.

Declaration of Competing Interest

None.