CHARACTERIZING AND MANAGING AN EPIDEMIC: A FIRST PRINCIPLES MODEL AND A CLOSED FORM SOLUTION TO THE KERMACK AND MCKENDRICK EQUATIONS

Supplement 1. The complete SIR model and its closed-form solution

Equations are numbered sequentially within each supplement. Equation numbers from the main text appear in bold.

Supplement 2. Verification of the CSIR model

Equations S1-38 and S1-39 can be combined purposefully to define an important relationship: where and t₀ is the time between the start of the epidemic and the time at which the first data point in the epidemic is measured.

The RCO expression is convenient in that it transforms Equation S1-38 into an equation linear in K₁. We call the expression the rate of change operator because it is the rate of change of new cases scaled by the current total of cases. It is Equation S2-1 that we fit to the nine data points after the imposition of containment measures in the country data in Figure 1. We also used this equation to find the parameters in Table 1.

We developed the solutions to Equations S1-1 through S1-5 assuming that K_T is a constant dependent only on the disease and is the same for all countries. We also assumed that P_c(t) is dependent on the population’s behavior. We used these assumptions when generating the correlations of the case data for various countries in Figures 2 and 3. Both of these assumptions can be tested independently of the country correlations.

To demonstrate that K_T is indeed a constant, we need to first further define P_c(t). As previously defined, P_c(t) is the number of specific infectious contacts a member of subpopulation N(t) has across the entire population. This is a function of the population’s behavior. Initially, we assume this is a function of population density and further, that people’s mobility extends over a constant average effective area per unit of time. We define this area as the effective area rate, A_1r(t). Using these definitions, we can write an expression for P_cr(t): where N_p = the entire population of the region with the infection, A = the area of the region, and = the population density. From Equation S2-2, we can see that P_cr(t) is proportional to both the population’s behavior, A_1r(t), and the population density, .

Similar to the way we defined P_c(t) using P_cr(t), we now define a quantity, A₁(t), in terms of A_1r(t): where A₁(t)is the effective specific area traversed by an individual. In this case, “specific” has the same meaning as it has for P(t)_c; that is, each person traverses only and exactly the same area for the duration of the time under consideration. We also call A₁(t) the “effective area” because the population is typically only dispersed within ∼1% of the land within a given region of a country (12). If we take this into account, then .

From the preceding discussion, we can now write an expression for K₁(t): The physical meanings of P_c(t) and A₁(t) may not be intuitive. As defined, neither are rates, but they can both vary in time and their values depend on the population’s behavior. They are, respectively, the number of specific infectious people who have been contacted and the specific area traversed by any index person within a given time interval. P_c(t) and A₁(t) are constant when the number of specific infectious people or the traversed area remain constant. However, if different people are contacted within a given time interval, the rates they depend on change, and therefore, P_c(t) and A₁(t) may change even if the total number of people contacted or area covered did not change during that time interval.

We can now check the assumption that K_T is a constant by substituting Equation S2-4 into Equation S1-39 and solving for K_Tt. Doing this, we find the following expression: If we define , then we can also write this expression as If K_T is a constant, then Equation S2-6 predicts that F(N(t)) is a linear function of time. Excepting A₁, all the quantities on the left-hand side of Equation S2-5 can be found for each country in the time before containment measures are enacted and are listed in Table 2. Therefore, we can estimate whether K_T is a constant by first using Equation S2-5 (or S2-6) to find if there is a constant value of A₁ that creates a straight line. An implicit assumption in this process is that the behavior of the population, P_c, is constant and therefore A₁ is constant, during the initial phase of the epidemic before containment measures were put in place.

Applying these assumptions, through a process of iteration, a value was found for A₁ (= 0.48 km²) which created a straight line with a correlation coefficient of 0.96 (Figure 4). Using Equation S2-5 (or S2-6) we can then determine that the slope of the line in Figure 4 indicates that the value of K_T is 0.26. This analysis strongly supports the assumption that is a constant.

Independent evidence that P_c is a measure of the population behavior can be developed by first using Equations S2-1 and S3-8 (Supplement 3) to predict that the RCO measure will be proportional to . We reasoned that if an independent measure of people’s mobility during the epidemic could be found and was linearly related to the RCO, we could have additional confidence in the veracity of the CSIR solution.

Google has compiled different measures, derived from mobile phone data, of people’s mobility (11). One of these measures is termed the Residential Mobility Measure (RMM). The RMM is a measure of the percentage change in the degree to which people stayed in their residence during the pandemic relative to a baseline measured over 5 weeks starting on January 3, 2020. Since and the RMM are both inversely proportional to the population’s mobility, we hypothesized that the RMM would be a good proxy for the value of . To test this, we plotted the integral over time of the daily RMM for the six countries whose data we analyzed, against the daily RCO. These plots appear in Figure 5, which clearly validates the hypothesized linear relationship.

Supplement 3. An analysis of the ASIR model

In this supplement, we address the questions: 1) What are the implied assumptions within the conventional approximation to Kermack and McKendrick’s original equations (1), known as the SIR model? and 2) How do β and γ need to vary with time for the ASIR model equations to mimic the CSIR model?

The SIR model, with non-time-varying parameters β and γ, is described by the following equations using the same definitions of S, I, and R as in Supplement 1: where N_p = total number of people in the population, β = rate of contact and transmission, t_r = time of infectiousness, and γ = rate of recoveries = 1/t_r.

This formulation from Kermack and McKendrick (1) is widely promoted as illustrative of the qualitative behavior of the more complex model with time-varying parameters presented in the paper. Because the assumption that β and γ are constant strongly limits the capability to describe epidemics, we call this form of the model the Approximate SIR (ASIR) model. While not intended to model epidemic dynamics perfectly, Equations S3-1 to S3-4 have always been assumed to reflect the general trends of an epidemic (1).

To understand the implications of the ASIR approximation, we must first understand the implications of the assumption that φ(t) and ψ(t) are constant over time. The implications are easily illustrated using Equations S1-49 and S1-50. With these equations and the prior definitions that and , we can find expressions for the time varying K_T(t) and P_c(t) when φ(t) and ψ(t) are assumed to be the constants φ and ψ: These two equations answer the first question posed at the beginning of this section. The implied assumption in the ASIR model formulation is that K_T and P_c vary with time in the manner described by Equations S3-5 and S3-6.

To apply the conditions in S3-5 and S3-6 to the CSIR solution, we must first rederive the CSIR solution under the assumption that both K_T(t) and P_c(t) do vary with time. We do this following the procedure used in Supplement 1, which results in the following set of relationships for N(t), , I(t) and R(t): The results of the simulation of both N(t) and are shown in Figure 7. We created these plots by imposing the functional relationships in Equations S3-5 and S3-6, as well as using Equations S3-7 and S3-8. In Figure 8, we then used Equations S3-5 and S3-6 to plot the time variation of both K_T and P_c resulting from the assumption that φ and ψ are the constants β and γ, respectively.

We now explore the opposite construct, the time varying nature of φ(t) and ψ(t) (and therefore β and γ) required when K_T and P_c are constants in Equations S1-49 and S1-50. Similar to the preceding analysis, we used simulations of the CSIR solution, Equations S1-49 and S1-50, and the ASIR model. Since the CSIR solution predicts the country data well, we also compared the simulation results to the data from two countries (Italy and New Zealand) and demonstrate that only if φ(t) and ψ(t) vary as prescribed in Equations S1-49 and S1-50 can the ASIR approximation model the country data with any degree of accuracy.

We first used the values of K₁ and K₂ in Table 1 to derive values of K_T and K₁ for Italy and New Zealand and used the CSIR solution to project the results. We then simulated the ASIR model with the assumption that the values of φ and ψ (and therefore β and γ) were constant and directly related to the values of K_T and K₁ used in the CSIR solution. The results of the CSIR simulation and the ASIR simulation are plotted in Figure 9, along with the country data. The CSIR model accurately models the country data and the ASIR model does not.

In a second step, we recast Equations S1-49 and S1-50 in terms of β(t) and γ(t): Using Equations S3-11 and S3-12, we then calculated the time series of β and γ necessary to generate the curves in Figure 9 and plotted them in Figure 10. Next, we used the values of β and γ plotted in Figure 10 in the ASIR model to generate the curves in Figure 11. Clearly, when β and γ are allowed vary in accord with the time variations of φ(t) and ψ(t), exposited above, the ASIR equations model country data quite well.

Equations S3-11 and S3-12 answer the second question posed at the beginning of this supplement. If β and γ are forced to vary according to these equations, then the ASIR model will model the epidemic. In a sense, Equations S3-11 and S3-12 are corrections to the ASIR model.

Supplement 4. Important relationships derived from the CSIR solution

Taking the limit of Equation S1-39 as t → ∞, we obtain the expression for the total number of individuals who will be infected throughout the entire epidemic: Mathematically, Equation S4-1 demonstrates that the total number of infections expected, N_∞, and the behavior of the population, P_c, are interrelated in the sense that changes in P_c have an exponential effect on the final number of infections. The exponential nature of the relationship underscores that small changes in population behavior dramatically affect the epidemic’s outcome. However, the reciprocity also means that the eventual number of cases produced by the epidemic is not foreordained, but rather, a strong function of interventions introduced.

We gain insight into the meaning of the CSIR model solution by looking at Equation S1-43 in more detail. Equations S1-38 and S1-40 can be used to rewrite Equation S1-43 as The left-hand side of Equation S4-2 is the rate of change in the number of new infections per person currently infected. The first two terms on the furthest right-hand side of Equation S4-2, , describe the net rate of successful infection. Since K_T is the rate at which an infected person causes infections per infectable contact, the terms must represent the rate of recovery per infected person.

We gain an additional intuitive insight about the solution from the following relationship, derived from Equations S1-38, S1-39, and S4-1: In words, the form of Equation S4-3 is or Equations S4-4 and S4-5 illustrate the logic of the CSIR solution in terms of probabilities.

Finally, using the prior definition, , and if t₀ = 0, we can use Equation S4-1 to write this simple expression for the CSIR solution for total cases:

Supplement 5. Controlling epidemics early

The quantitative mathematical relationships derived from the CSIR solution in Supplements 2 and 4 characterize the dynamics of an epidemic and illustrate that strong and early intervention is critical. Equation S4-1 quantifies that the ultimate number of individuals infected in an epidemic, N_∞, will be exponentially dependent on the number of people with whom each person interacts.

The real-world country data provide vivid examples. Both South Korea and New Zealand enacted strong and early interventions compared to other countries (3, 4), as reflected by their K₁ values (Table 1). These strong interventions led to earlier peaks in new cases and to far fewer total cases than in other countries (Figures 2 and 3): the peak number of new cases in both South Korea and New Zealand was 90–99% lower than in other countries, a compelling validation of the explicit statement in the CSIR solution that strong intervention leads to exponentially more favorable outcomes.

In the USA, interventions that began on March 16 started to have an effect around March 23, 2020 (Figure 3B); the number of active cases on March 23, 2020 (6) was 46,136 (Table 1). Using the values of K₁ and K₂ from Table 1, Equation S4-1 predicts that the ultimate number of cases would have been approximately 1.22 million. If the same intervention had been implemented and sustained starting on March 10, when there were 59 times fewer (782) cases (6), the model predicts that the ultimate number of cases would also have been 59 times lower, or 20,725. Thus, earlier action could have reduced the ultimate number of projected cases by more than 98%. Of course, the projected estimate of approximately 1.22 million total USA cases would only have occurred if the effectiveness of the interventions that were launched on March 16 had been sustained. Unfortunately, a marked reduction in effective interventions occurred in many parts of the USA in mid-April, well before the official reopening of the economy (13). This caused a second surge in new cases in late April and is why the observed data and the model prediction diverge in Figure 3B.

As shown in Supplement 4, the CSIR solution provides an estimate of the time to the peak of new cases, t_max. Using Equation S3-10 and the values of K₁ and K₂ from Table 1, the predicted peak in new cases in the USA would have occurred near March 24 if the intervention had begun on March 10. Instead, a 6-day delay in effective intervention shifted the initial peak to April 11, 16 days later, as projected, and that peak was much higher (Figure 3B).

As shown, too, in Supplement 4, epidemic acceleration, the instantaneous potential to change the pace of the epidemic, can be determined at any point in the epidemic and depends on the social containment actions in effect at that time (Equation S4-11). What is perhaps less apparent, but predicted by the model, is that two countries with identical numbers of cases on a given day can, in fact, have different accelerations on the same day, and exhibit different dynamics immediately after that day.

South Korea and New Zealand (Figure 2A and F) had nearly identical case counts when each imposed strong containment measures (204 cases in South Korea on February 21, and 205 in New Zealand on March 25). Their models suggest that their interventions were about equally effective (K₁ = 0.24 in South Korea and 0.17 in New Zealand; see Table 1). However, since South Korea has a much higher population density than New Zealand ((7), data in Table 2), it had a much higher number of interactions when the interventions were imposed and, therefore, a higher rate of acceleration, as evidenced by its higher RCO at the time of intervention. Indeed, the rate of change of new cases was higher in South Korea than in New Zealand, and the later number of cases in South Korea was higher than in New Zealand (Figure 2A and 2F).

Equation S4-8 clearly illustrates these lessons. As social distancing is strengthened (lower P_c and therefore, higher K₁), the Effective Replication Number decreases, and the epidemic slows. Early and strong interventions, especially in countries with indigenously high levels of social interaction, are necessary to stop an epidemic in the initial stages. Re-openings, enacted too early, can reignite the epidemic, dramatically increasing the number of cases. The astonishing magnitude of the effects, driven by only a few days of delay, derives from the doubly exponential nature of the underlying relationships.

Supplement 6. Ending an ongoing epidemic

We can use the CSIR solution to design measures to end an epidemic in an advanced stage. The management plan is built by first using Equation S4-13 to predict the number of days a given level of intervention, K₁, is needed to reduce the new daily cases by a target fraction: where t_target is the time to the desired reduction and D_tf is the target as a fraction of the current level of new cases per day.

For example, a country targeting a 90% reduction of new cases per day (e.g., from 50,000 to 5,000 cases per day, D_tf = 0.1), can attain its target in about 12 days by imposing a containment level of K₁ = 0.2. The South Korea and New Zealand data demonstrate that Equation S4-15 is valid and that K₁ = 0.2 is achievable for this duration. Both countries achieved a value of K₁ close to 0.2 for the time necessary to produce a 90% reduction. It took 13 days in South Korea (March 3–16) and 15 days in New Zealand (April 2–15), New Zealand (6).

Since , K₁ = 0.2 characterizes a lockdown in which people in a country can each have only one plausibly infectious contact with a little over one specific person for the containment duration. This does not mean they cannot contact anyone other than the one person; but they must use care, masks and proper distancing, to ensure there is no plausibly infectious contact with anyone other than the one person.

Returning to the planning example, after achieving the initial 90% reduction, a reasonable next step might be to relax social containment to a level that allows the economy to remain viable, while preventing the epidemic from erupting again. We can again find the level of K₁ necessary to achieve a chosen target, using Equation S4-13. If an additional 90% reduction in new cases per day is desired, and a period of 90 days is tolerable for that reduction, then a new level of approximately K₁ = 0.025 is needed. This equates to a 90-day period during which each person can be in contact with seven specific people, in an infectable way. Note that this is three times less stringent than the original USA shutdown level in April 2020 as shown by the level of K₁ calculated for the United States in that period (Table1). Thus, with a well-planned approach, a country can reduce its new daily cases by 99% in approximately 100 days, enabling the country to control, and essentially end the epidemic, while simultaneously maintaining economic viability.

If even K₁ = 0.025 is too restrictive, we can choose a still lower K₁, but it must be large enough to avoid a new outbreak. A bound for the new value of K₁, low enough to prevent an outbreak and yet continue decreasing the new cases per day, can be found using Equation S4-12.

We can easily monitor the progress of interventions using the RCO, as the curve for South Korea illustrates (Figure 1A). Had this country maintained the implemented level of distancing measures, the data would have followed the initial slope. However, the actual data departed from the slope, heralding failures in (or relaxation of) social distancing, which were later documented to have occurred during the indicated time frame (3) (circled data, Figure 1A). Because it summarizes epidemic dynamics, we can use the RCO to continuously determine the effectiveness of implemented measures and whether they need adjustment.

Supplement 7: Understanding the assumption of contiguousness

To enhance the understanding of the CSIR solution and explain what the assumption of a contiguous epidemic means, we applied a perspective derived from the expression for N_∞ in Equation S4-1 to Equations S1-1 to S1-5. Rather than look at the epidemic as involving the total population, N_p, this perspective focuses on only that portion of the population that will eventually become infected in the epidemic. This is the subpopulation N_∞. We also introduce a subpopulation of N_∞ called N_S(t), which we define as the sub-population at time t that is in contact with the epidemic and under threat of infection. This change in perspective views the epidemic mathematically from inside the bounds of the already infected population, which is ever expanding, rather than mathematically describing what is happening within a fixed total population.

The reinterpretation begins by first recognizing that, in Equations S1-1 to S1-3, the initial number of infections introduced to the population—which we designate I_i—is merely the starting value of the epidemic and can be any value at all. The second step is to imagine that at every time, t, the epidemic starts again and the initial infections introduced into the population, I_i(t), are equal to the then-current infections. The third step is to recognize that the values of the integrals in Equations S1-3 to S1-5 are equal to zero and B(θ) =1 whenever the epidemic starts. With these perspectives in mind, if Equations S1-1 to S1-3 and S1-5 are first all multiplied by A_P, they can be rewritten as where σ(t) = φ(t) N_S(t) and ω(t) = ψ(t),

Su(t) is the remaining portion of the population that will become infected, and since there are only susceptible people in N_S(t) when the new infections are introduced, Of course, .

If we also now define a quantity, R_i(t), as the number of recovered persons introduced to the population at the same time as I_i(t), we can write an equation for N_∞ and define a new quantity, N(t), as the number of people either currently infected or recovered: We have not assumed that any of the quantities are constant; therefore, we can now explicate the time-varying nature of these quantities. We begin by considering what happens to Equations S7-1 to S7-3 during the time interval Δt from a time t to t +Δ t, and by rewriting these equations as difference equations: Taking the limit as Δt → 0, we obtain the following differential equations: From the preceding, it can be easily seen that I_i(t) = I(t) and R_i(t) = R(t). Note that Su(∞) = 0 unlike S(∞).

The perspective in the immediately preceding part of the analysis is that the epidemic can be considered to start over again at each instant in time. In this perspective, the susceptible population Su(t) is not fixed by initial conditions, but rather is the population that will eventually become infected. Embedded in this concept is the assumption that during each Δt, the susceptible population is always in contact with those people who have been previously infected or who will become infected. That is, the epidemic remains contiguous. This is not restrictive when considering the initial stage of the epidemic; however, as explained in Supplement 6, it plays a critical role in recognizing an outbreak that starts outside the population N_∞.

We have recast the equations by shifting the population and area under consideration to N_∞ and A_∞, respectively, rather than N_p and A_p. Further, we have defined the susceptible portion of the population during the epidemic as those who will eventually become infected under the conditions in place at each instance in time. This shift in perspective retains the mathematical equivalence to the equations derived by Kermack and McKendrick (1) because the portion of N_p that is not a part of N_∞ never becomes infected, and therefore never affects the values of I(t) or R(t) in the solutions to Equations S1-1 to S1-5 or S1-8 to S1-11.