The SIR dynamic model of infectious disease transmission and its analogy with chemical kinetics

The reactions

Susceptible folks can contract the disease if they come into contact with an infectious individual. Once infected, they move into the infectious compartment, assuming zero delay between infection and the ability to transmit the disease. This is analogous to an irreversible, autocatalytic chemical reaction (Schuster, 2019; Scott, 1983) between a reactant, S, and catalyst, I: {1} $$\mathrm{S}+\mathrm{I}\to 2\mathrm{I}$$

Infectious individuals eventually recover or die from the disease, entering the removed compartment. Folks in the removed compartment do not participate in disease transmission. That is, they cannot transmit the disease, nor can they contract it again, assuming that recovery from the disease confers immunity to reinfection. This is analogous to a reaction where the catalyst, I, irreversibly degrades or converts to a deactivated form, R: {2} $$\mathrm{I}\to \mathrm{R}$$

We assume permanent protective immunity is conferred upon recovery from the disease, thus neglecting the possibility of an R → S reaction.

So, the SIR model of an epidemic is analogous to an autocatalytic reaction (rxn. {1}) with catalyst deactivation (rxn. {2}). An infectious individual (the catalyst, I), (i) upon contacting (colliding with) a susceptible individual (the reactant, S), can convert them into another infectious individual (another catalyst particle) and (ii) recovers or dies (deactivates) with time. Figure 1 depicts the flow of individuals between compartments under the SIR model, induced by rxns. {1} and {2} .

The dynamic mathematical model

Mathematically, the SIR model (Kermack & McKendrick, 1927; Edelstein-Keshet, 1988; Murray, 1993; Martcheva, 2015a) is equivalent to a dynamic model of a well-mixed, isothermal batch reactor carrying out the two homogeneous, elementary rxns. {1} and {2} .

As in a (closed) batch reactor, we neglect immigration and emigration (hence, the absence of flow to/from external populations in Fig. 1). Moreover, we take births and deaths not caused by the disease to be negligible over the (short) time scale of the epidemic.

Let [S](t), [I](t) and [R](t) be the fraction of the population that is susceptible, infectious, and removed, respectively, at time t. Considering a large population, we treat [S], [I], and [R] as continuous variables. So, [S], [I], [R] ∈ [0, 1].

The replacement and basic reproduction numbers

Two important numbers aid our characterization and understanding of SIR model dynamics: the replacement number, r, and the basic reproduction number, ${\mathcal{R}}_0$.

The replacement number, r = r(t), is the expected number of folks (directly) infected by a typical infectious individual, mixing in the population, over the course of their infectiousness (Hethcote, 2000). Because the concentration of susceptible folks [S] = [S](t) influences the frequency that a typical infectious individual contacts a susceptible individual, r changes over the course of an epidemic. In the SIR model, a typical infectious individual is expected to produce β [S](t) new infections per unit time (incidence rate per infectious individual) for an infectious duration of γ⁻¹. The replacement number is therefore: (4) $$r=r(t)={\displaystyle \frac{\mathrm{\beta }}{\mathrm{\gamma }}[\mathrm{S}](t).}$$

The basic reproduction number, ${\mathcal{R}}_0$, is defined as the initial replacement number when one infectious individual is introduced into an all-susceptible population (Edelstein-Keshet, 1988; Hethcote, 2000). In this context, ${\mathcal{R}}_0$ in the SIR model is the replacement number in Eq. (4) when [S] ≈ 1: (5) $${\mathcal{R}}_0={\displaystyle \frac{\mathrm{\beta }}{\mathrm{\gamma }}.}$$

That is, ${\mathcal{R}}_0$ is the expected number of infections directly caused by a single infectious individual introduced into an entirely susceptible population (Hethcote, 2000).

The dimensionless numbers r and ${\mathcal{R}}_0$ are properties of both the disease and the population (Delamater et al., 2019). While r = r(t) changes with time, ${\mathcal{R}}_0$ is constant and defined only at the initial stage of a particular context: when one infectious individual is introduced to an all-susceptible population. Notably, the two numbers are related via $r(t)={\mathcal{R}}_0[\mathrm{S}](t)$.

If the basic reproduction number ${\mathcal{R}}_0$ in Eq. (5) is large (small), the infected are infectious for a long (short) period of time, the disease is (not) easily transmitted, and/or the average frequency of contacts in the population is high (low). Under the analogy with chemical kinetics, since the activity and longevity of the catalyst, I, are embedded in β and γ, respectively: ${\mathcal{R}}_0$ is large (small) if the catalyst has a high (low) activity and/or remains active for a long (short) time. These remarks also hold for the replacement number, $r={\mathcal{R}}_0[\mathrm{S}]$. However, r decreases as the concentration of the reactant, [S], decreases, owing to the reduced frequency that any given catalyst particle encounters a reactant particle to catalyze its conversion into another catalyst particle by rxn. {1} .

The replacement number r(t) determines the sign of [I]′(t)

The replacement number r(t) in Eq. (4) is key to understanding SIR model dynamics. By inspection of Eq. (2), [I](t) is increasing at time t if the replacement number r(t) > 1 and decreasing if r(t) < 1. This is intuitive: if the typical infectious person mixing in the population is expected to infect less than one susceptible person before they recover, they are not expected to replace themselves with a new infectious individual to propagate the disease, and, consequently, [I](t) is decreasing. Note $r(t)<1\iff [\mathrm{S}](t)<{\mathcal{R}}_0^{-1}$.

Does an epidemic ensue?

We first address a qualitative question: given the initial conditions in Eqs. (6)–(8), does the disease invade the population? The outcome depends on the initial replacement number $r_0:=r(0)={\mathcal{R}}_0[\mathrm{S}{]}_0$ in a threshold manner. If r₀ > 1, [I](t) initially increases. That is, the disease spreads; an epidemic ensues. If r₀ < 1, [I](t) initially decreases and, since [S](t) is a decreasing function of time (see Eq. (1)), decays to zero. That is, the disease dies out; an epidemic does not ensue.

Under the chemical reaction analogy, if r₀ < 1 (r₀ > 1), the injected catalyst particles deactivate via rxn. {2} faster (slower) than they catalyze rxn. {1} to convert the reactant, S, into more catalyst, I, to propagate autocatalytic rxn. {1} ; the reaction dies out (proceeds).

For our remaining analysis, we take r₀ < 1 and further analyze the dynamics of an SIR epidemic.

Initial exponential growth

Early in the epidemic, the number of infectious folks grows approximately exponentially with growth rate (r₀ − 1)γ: (9) $$[\mathrm{I}](t)\approx [\mathrm{I}{]}_0{e}^{(r_0-1)\mathrm{\gamma }t}.$$

This follows from Eq. (2) if we neglect the depletion of the susceptible pool by approximating [S](t) ≈ [S]₀, valid only in the initial stage of the epidemic; as the disease spreads, [S] decreases and diminishes the replacement number. Equation (9) is thus an overestimate.

Since Eq. (9) is also a valid initial approximation for r₀ < 1, it reinforces that an epidemic will not ensue if r₀ < 1, since [I](t) would then decay approximately exponentially.

A simulation of an SIR model epidemic

Figure 2 shows the numerical solution to Eqs. (1)–(3) with initial conditions in Eqs. (6)–(8) for ${\mathcal{R}}_0=2$ and [I]₀ = 10⁻⁵. Initially, the concentration of infectious folks, [I], grows approximately exponentially (see Fig. S1 for a comparison of [I](t) with Eq. (9)). The concentration of susceptible folks, [S], decreases monotonically as the disease invades the population. As a result, the frequency with which any given infectious individual comes into contact with a susceptible individual decreases. In conjunction with the infectious folks recovering, this eventually causes the (net) growth rate of [I] to diminish and, when the replacement number r in Eq. (4) drops below one, causes [I] to decay. The epidemic self-extinguishes, that is, ${lim}_{t\to \mathrm{\infty }}[\mathrm{I}](t)=0$ (Bjørnstad et al., 2020a). The R category accumulates the folks that have been infected by and have recovered or died from the disease, $[R](t)=\mathrm{\gamma }{\int }_0^t[\mathrm{I}](\mathrm{\tau })d\mathrm{\tau }$ (see Eq. (3) with initial condition 8). Notably, the disease does not infect the entire population, even after an infinite amount of time (${lim}_{t\to \mathrm{\infty }}[\mathrm{S}](t)\ne 0$). That is, the epidemic self-extinguishes not because the population is depleted of susceptible folks, but rather because it is depleted of infectious folks (Murray, 1993; Weiss, 2013). An alternative visualization of SIR model dynamics in Fig. S2 emphasizes [S](t) + [I](t) + [R](t) = 1, ∀t ≥ 0.

Under the chemical kinetics analogy, the autocatalytic rxn. {1} begins to slow down after the concentration of the reactant, [S], decreases sufficiently. Once the replacement number r(t) drops below one (i.e., once [S] drops below ${\mathcal{R}}_0^{-1}$), the catalyst, I, is deactivating via rxn. {2} faster than it is converting the reactant, S, into more catalyst to replenish itself via rxn. {1}. Consequently, the catalyst concentration, [I], begins to drop and decays to zero. Owing to catalyst deactivation (rxn. {2}), not all reactant, S, is consumed, even after an infinite amount of time.

Solution trajectory in the ([S], [I]) phase space

We can analytically find the trajectory of the solution to Eqs. (1) and (2) in the ([S], [I]) phase plane. Dividing Eq. (2) by Eq. (1) takes us into the phase plane by giving a differential equation with a ${\displaystyle \frac{d[\mathrm{I}]}{d[\mathrm{S}]}}$ derivative, with the view of [I] as a function of [S]. Separating, integrating, and applying the initial conditions in Eqs. (6) and (7), we arrive at the solution path (Brauer, 2008): (10) $$[\mathrm{I}](t)=1-[\mathrm{S}](t)+{\displaystyle \frac{1}{{\mathcal{R}}_0}\log \left({\displaystyle \frac{[\mathrm{S}](t)}{{[\mathrm{S}]}_0}}\right).}$$

Figure 3 shows the trajectory of the solution given by Eq. (10) (${\mathcal{R}}_0=2$, [I]₀ = 10⁻⁵), traveling from the initial condition in the bottom right corner to the final condition on the bottom left. The trajectory reinforces that [S](t) decreases monotonically with time; that [I](t) increases, peaks, then diminishes to zero; and that a fraction of the population remains susceptible after the epidemic dies out.

Final size of the epidemic

Given the epidemic (autocatalytic reaction) dies out before all of the susceptible folks (reactant) have been infected, what fraction of the population (reactant) remains susceptible (unreacted) after the epidemic ends?

We find an implicit equation for $[\mathrm{S}{]}_{\mathrm{\infty }}:={lim}_{t\to \mathrm{\infty }}[\mathrm{S}](t)$ by taking the limit $t\to \mathrm{\infty }$ in Eq. (10): (11) $$0=1-[\mathrm{S}{]}_{\mathrm{\infty }}+{\displaystyle \frac{1}{{\mathcal{R}}_0}\log \left({\displaystyle \frac{{[\mathrm{S}]}_{\mathrm{\infty }}}{{[\mathrm{S}]}_0}}\right).}$$

We used the fact that the epidemic eventually dies out, ${lim}_{t\to \mathrm{\infty }}[\mathrm{I}](t)=0$. [S]_∞ is the unique root of Eq. (11) in $(0,{\mathcal{R}}_0^{-1})$ (Hethcote, 2000). The fraction of the population infected over the course of the epidemic is ${lim}_{t\to \mathrm{\infty }}[R](t)=:[R{]}_{\mathrm{\infty }}=1-[\mathrm{S}{]}_{\mathrm{\infty }}$ since the R category accumulates those that have recovered or died from the disease. Figure 4A shows that, as ${\mathcal{R}}_0$ increases from one, more of the population will be infected over the course of the epidemic.

Peak prevalence of infectious folks

The peak prevalence of infectious folks, max_t [I](t), is important because it determines the maximum strain on the healthcare system. Both Eqs. (10) and (2) show that the maximum in [I](t) (look for ${\displaystyle \frac{d[\mathrm{I}]}{d[\mathrm{S}]}=0}$ or ${\displaystyle \frac{d[\mathrm{I}]}{dt}=0}$, respectively) occurs when the replacement number r(t) = 1. Before the peak is reached, the replacement number r(t) > 1 and [I](t) is increasing; after the peak, r(t) < 1 and [I](t) is decreasing to zero. Substituting $[\mathrm{S}]={\mathcal{R}}_0^{-1}$ into Eq. (10), we find: (12) $$\underset{t}{max}[\mathrm{I}](t)=1-{\displaystyle \frac{1}{{\mathcal{R}}_0}\left[1+\log \left({\mathcal{R}}_0{[\mathrm{S}]}_0\right)\right].}$$

Figure 4B shows how the peak prevalence of infectious folks increases as ${\mathcal{R}}_0$ increases from one. In addition, Fig. 4C shows that the (non-dimensional) time to reach peak prevalence, argmax_t α[I](t), decreases as ${\mathcal{R}}_0$ increases from one.

Herd immunity

A population achieves herd immunity when a sufficient fraction are immune to the disease so as to confer indirect, population-level protection from an invasion of the disease upon the introduction of an infectious individual. Notably, the immunity could be acquired by either previous infection or by vaccination (Hethcote, 2000).

The administration of a vaccine that confers perfect and permanent immunity to a susceptible individual is modeled by introducing an S → R reaction to the SIR model, that is, by allowing flow in Fig. 1 from the susceptible category directly to the removed category (Yusuf & Benyah, 2012; Martcheva, 2015a; Shulgin, Stone & Agur, 1998).

What fraction v of a population must be immune to the disease in order to achieve herd immunity? [I](t) will decrease upon introducing an infectious individual if the replacement number r(t) < 1, that is, if $[\mathrm{S}]<{\mathcal{R}}_0^{-1}$. Thus, a fraction $v>1-{\mathcal{R}}_0^{-1}$ of the population must be immune to achieve herd immunity. Figure 5 shows the region $v>1-{\mathcal{R}}_0^{-1}$ (gray) and illustrates that, as ${\mathcal{R}}_0$ increases from one, more of the population must be immune to prevent an epidemic via herd immunity.

In the chemical kinetics analogy, herd immunity results from reducing the concentration of the reactant, [S], so that a catalyst particle, I, fed to the reactor is expected to deactivate before it encounters an S particle and autocatalyzes rxn. {1} to replace itself. That is, to achieve herd immunity, [S] must be reduced sufficiently to make the replacement number r less than one.

We now make a comparison between herd immunity achieved exclusively (i) by vaccination (S → R) and (ii) by infection and removal (recovery or death) (S → I → R). Consider an all-susceptible population. Under strategy (i), vaccinating a fraction $v=1-{\mathcal{R}}_0^{-1}$ of the population will suffice to confer herd immunity. Under strategy (ii), by introducing infectious individuals and allowing the epidemic to run its course, a fraction [R]_∞ of the population will be infected over the course of the epidemic. Though, herd immunity will be achieved when $[\mathrm{S}]={\mathcal{R}}_0^{-1}$. At this point, [I](t) is maximal (see Eq. (10) and Fig. 3). These infectious individuals will infect more susceptible folks before they all recover/die, even though, on average, each is expected to recover/die before they infect a susceptible individual (r < 1 after this point). The comparison of $v=1-{\mathcal{R}}_0^{-1}$ and [R]_∞ in Fig. 5 confirms that a significantly greater fraction of the population will be infected by the disease if we let it run its course than would have needed vaccination to achieve herd immunity. However, the most obvious benefit of the path S → R over S → I → R is avoidance of disease-induced suffering and, possibly, death.

Conclusion: ${\mathcal{R}}_0$ is a useful tool

The dimensionless basic reproduction number ${\mathcal{R}}_0$ is a property of an infectious disease within a population. The average frequency of contacts of an individual in the population, the transmissibility of the disease, and the average duration of infectiousness are all embedded in ${\mathcal{R}}_0$. In the SIR model, ${\mathcal{R}}_0$ influences whether or not the disease invades the population, the initial exponential growth rate of infectious folks, how many are infected over the course of the epidemic, the peak prevalence of infectious folks, and how many must be vaccinated to achieve herd immunity.

Births, deaths and loss of immunity

To model infectious disease transmission over longer time scales, we can modify the SIR model to account for births and deaths not caused by the disease by allowing flow into the S compartment and flow out of all compartments, respectively. See Fig. 6A (top). While Eqs. (1)–(3) resemble a dynamic model of a closed batch reactor carrying out rxns. {1} and {2}, the modified SIR model with demographics resembles a dynamic model of a continuously stirred tank reactor (Fogler, 2010): births are represented by a feed stream of pure S flowing into the reactor, while deaths are represented by an effluent stream. In addition, as opposed to assuming recovery from infection confers permanent immunity, we can model temporary immunity by introducing an R → S reaction to represent loss of immunity. Because births and the loss of immunity continuously add to the susceptible pool, [I](t) in the SIR model with births, deaths, and loss of immunity can exhibit damped oscillations that settle on an endemic equilibrium, where the disease remains in the population indefinitely (i.e., [I] is maintained at a non-zero value). See Fig. 6A (bottom). At an endemic equilibrium, the replacement number is maintained at one so that [I](t) is neither increasing nor decreasing (Hethcote, 2000; Bjørnstad et al., 2020b; Earn, 2008).

Time-varying parameters

A time-varying transmission rate constant β = β(t) can model (voluntary or policy-induced) changes in behavior during an epidemic that affect (i) the average frequency of contacts of an individual in the population or (ii) the transmissibility of the disease. Examples include staying home, social distancing, and taking hygiene measures (Fenichel et al., 2011; Reluga, 2010; Bjørnstad et al., 2020a; Bootsma & Ferguson, 2007). A periodic β(t) can model seasonality of an infectious disease (Fisman, 2007; Aron & Schwartz, 1984). A time-varying recovery rate constant γ = γ(t) can model changes in the average time period of infectiousness, for example, a reduction of γ⁻¹ by administering a drug to infected patients (Bjørnstad et al., 2020a). Evolution of the pathogen could also cause β and γ to change with time (Lion & Metz, 2018). Reducing the transmission rate constant β and/or the average time period of infectiousness γ⁻¹ reduces the replacement number r (see Eq. (4)); if r is reduced below one, the prevalence of infectious folks will decrease. Analogously, in chemical kinetics, time-varying reaction rate constants arise when the temperature in the reactor changes with time (Fogler, 2010).

Stochasticity

As opposed to the deterministic differential Eqs. (1)–(3), we can introduce randomness into the SIR model to account for the stochastic and uncertain nature of human interaction and disease transmission (Allen, 2008, 2017; Britton, 2010). Stochastic epidemic models aim to describe the probabilistic distribution of outcomes, for example, the distributions of [R]_∞ and max_t [I](t) (Allen, 2008; Britton, 2010). Stochasticity can be particularly important to model for small populations and in the early stage of an epidemic, when there are small numbers of infectious individuals (Britton, 2010). Analogously, when modeling chemical reaction dynamics, if the reactants are not abundant, such as in a biological cell, discreteness and stochasticity can be important to model (Erban, Chapman & Maini, 2007; Gillespie, Hellander & Petzold, 2013). See Fig. 6B for an example of a stochastic simulation of the SIR model.

More realistic probabilistic distributions of the infectious time period

The probabilistic implication of first-order decay in Eq. (3) for a single infectious individual is that their time period of infectiousness is an exponentially distributed random variable with mean γ⁻¹ (Brauer, 2008). More realistic probabilistic distributions of the time period of infectiousness can be built into the SIR model (Hethcote & Tudor, 1980; Lloyd, 2001). In the realm of biochemical kinetics, such non-exponential waiting times could arise, for example, between gene transcription events, owing to several steps involving transcription factors, cofactors, chromatin remodeling, etc. (Pedraza & Paulsson, 2008; Haberle & Stark, 2018).

Additional compartments

We can introduce additional compartments and further divide the population to (i) account for heterogeneity in the population, (ii) model public health interventions, and (iii) include other characteristics of the disease. Members of different compartments can then have different contact patterns, transmissibilities, and/or recovery times. For example, the SEQIJR model (Gumel et al., 2004) introduces three additional compartments: Exposed (E), Quarantined (Q), and Isolated (J). Usually, the E compartment is included to model the latent period of a disease; exposed individuals have been exposed to and infected by the disease but are not yet infectious (Brauer, 2008; Martcheva, 2015b). The Q and J compartments are included to model the control measures of quarantining exposed individuals and isolating infectious individuals, respectively. Members of the Q and J compartments contact other members of the population with reduced frequencies, compared to those in the E and I compartments, respectively (Brauer, 2008). See Fig. 6C. As another example, we can partition the infectious compartment into two distinct categories of infectious folks: symptomatic and asymptomatic (Vivas-Barber, Castillo-Chavez & Barany, 2014), which can have different frequencies of contacts, transmissibilities, and recovery rates. Finally, to account for different mixing patterns among different age groups (Brauer, 2008), age-structured compartmental models partition the S and I compartments into age cohorts. Notably, introducing additional compartments to the SIR model is much like introducing additional reactive chemical species into a chemical kinetics model.

Additional compartments complicate the derivation of the basic reproduction number ${\mathcal{R}}_0$ from the model (Van den Driessche, 2017). Still, the peak prevalence of infectious folks and final size of the epidemic increase with ${\mathcal{R}}_0$ in most models (Van den Driessche & Watmough, 2008). Notably, a more precise definition of ${\mathcal{R}}_0$ is needed if compartments are introduced to account for heterogeneity in the population (Diekmann, Heesterbeek & Metz, 1990).

Spatial heterogeneity

We can model spatial heterogeneity of an epidemic, that is, spatially-dependent [S], [I] and [R], by treating space as discrete (Van den Driessche, 2008) or continuous (Martcheva, 2015a). Modeling the spatial movement of susceptible and infectious individuals in a continuous space as a diffusive process results in reaction-diffusion equations (Martcheva, 2015a), also commonly used in the chemical sciences. Compartmental, epidemic models of metapopulations allow travel between spatially segregated regions and resemble models of multiple reactors connected with pipes that permit flow between them (Van den Driessche, 2008; Arino et al., 2005; Balcan et al., 2010).

More detail/structure in the contact patterns

Network epidemiological models allow us to study how structure and heterogeneity in the contact patterns within a population impact the dynamics of an epidemic (Kiss, Miller & Simon, 2017; Bansal, Grenfell & Meyers, 2007; Keeling & Eames, 2005; Newman, 2002). In a typical contact network model, each individual in a population is represented by a node, and each interaction/contact between a pair of individuals is represented by an edge. The state of an SIR network at a given time is specified by the disease status (S, I, or R) of each node. The disease can be transmitted along an edge if one of the nodes is infectious and the other is susceptible (Kiss, Miller & Simon, 2017) (See Fig. 6D).

The dynamics of disease transmission on a contact network depend significantly on the structure of the network, for example, on its degree distribution, node clustering, and correlations between the degrees of connected nodes (Bansal, Grenfell & Meyers, 2007). Intuitively, the dynamics of an SIR model on large, static, random k-regular contact networks (a random network where each node has k neighbors) closely approximate those of Eqs. (1)–(3), which assume homogeneous mixing (Bansal, Grenfell & Meyers, 2007). In contrast, the dynamics of SIR models on highly heterogeneous contact networks, such as those with heavy-tailed degree distributions that capture “superspreaders” (Lloyd-Smith et al., 2005; Meyers et al., 2005), deviate significantly from Eqs. (1)–(3).

Interestingly, directed networks can model asymmetric contacts that can only transmit the disease one-way, such as donated blood transfusions (Keeling & Eames, 2005; Allard et al., 2020).

A contact network model such as in Fig. 6D only implicitly accounts for spatial heterogeneity. To more explicitly account for spatial heterogeneity, contact networks can be encoded as a dynamic, bipartite graph with two classes of nodes: individuals and locations. Edges then connect individuals with locations, and disease transmission can occur between people that are in the same location. (Eubank et al., 2004) Other agent-based models are fully spatially explicit and track the location of individuals as they move between households, schools, and workplaces (Ferguson et al., 2005). Notably, to construct such contact network models, very detailed data is required (Ajelli et al., 2010).

Vectors that transmit the pathogen from host to host

Some infectious diseases are primarily transmitted from one host to another host by living vectors that can acquire and carry the infectious agent. For example, mosquitoes can acquire an infectious agent (e.g., the virus that causes dengue fever or the parasite that causes malaria) from feeding on the blood of an infected human, then transmit the infectious agent to another, susceptible human when feeding on their blood. SIR-like models of diseases transmitted by vectors include an incidence term β [S][I_v], where [I_v] is the concentration of infectious vectors and β includes the frequency that the vector bites hosts and the probability of transmission conditioned upon a bite. (Martcheva, 2015c, Smith et al., 2012).