Editing Compartmental models (epidemiology) (section)

===The SIR model without birth and death===
[[File:SIR trajectory.png|thumb|400px|right|A single realization of the SIR epidemic as produced with an implementation of the [[Gillespie algorithm]] and the numerical solution of the ordinary differential equation system (dashed)]]
The dynamics of an epidemic, for example, the [[Influenza|flu]], are often much faster than the dynamics of birth and death, therefore, birth and death are often omitted in simple compartmental models.  The SIR system without so-called vital dynamics (birth and death, sometimes called demography) described above can be expressed by the following system of ordinary [[differential equations]]:<ref name="Beckley"/><ref name="Hethcote2000">{{cite journal |author=Hethcote H |title=The Mathematics of Infectious Diseases |journal=SIAM Review |volume=42 |issue= 4|pages=599–653 |year=2000 |doi=10.1137/s0036144500371907|bibcode=2000SIAMR..42..599H |s2cid=10836889 }}</ref>

:<math>
\left\{\begin{aligned}
& \frac{dS}{dt} = - \frac{\beta}{N} I S, \\[6pt]
& \frac{dI}{dt} = \frac{\beta}{N} I S - \gamma I, \\[6pt]
& \frac{dR}{dt} = \gamma I, 
\end{aligned}\right.
</math>
[[File:SIR model cartoon.png|thumb|right|The SIR model]]

where <math>S</math> is the stock of susceptible population in unit number of people, <math>I</math> is the stock of infected in unit number of people, <math>R</math> is the stock of removed population (either by death or recovery) in unit number of people, and <math>N</math> is the sum of these three in unit number of people. <math>\beta</math> is the infection rate constant in the unit number of people infected per day per infected person, and <math>\gamma</math> is the recovery rate constant in the unit fraction of a person recovered per day per infected person, when time is in unit day.

This model was for the first time proposed by [[William Ogilvy Kermack]] and [[Anderson Gray McKendrick]] as a special case of what we now call [[Kermack–McKendrick theory]], and followed work McKendrick had done with [[Ronald Ross]].{{citation needed|date=May 2023}}

This system is [[non-linear]], however it is possible to derive its analytic solution in implicit form.<ref name="Harko"/> Firstly note that from:

:<math> \frac{dS}{dt} + \frac{dI}{dt} + \frac{dR}{dt}  = 0,</math>

it follows that:

:<math> S(t) + I(t) + R(t) = \text{constant} = N,</math>

expressing in mathematical terms the constancy of population <math> N </math>. Note that the above relationship implies that one need only study the equation for two of the three variables.

Secondly, we note that the dynamics of the infectious class depends on the following ratio:

:<math> R_0 = \frac{\beta}{\gamma},</math>

the so-called [[basic reproduction number]] (also called basic reproduction ratio). This ratio is derived as the expected number of new infections (these new infections are sometimes called secondary infections) from a single infection in a population where all subjects are susceptible.<ref name=Bailey1975>{{cite book |author=Bailey, Norman T. J. |title=The mathematical theory of infectious diseases and its applications |publisher=Griffin |location=London |year=1975 |isbn=0-85264-231-8 |edition=2nd}}</ref><ref name=nunn2006>{{cite book |author1=Sonia Altizer |author2=Nunn, Charles |title=Infectious diseases in primates: behavior, ecology and evolution |publisher=Oxford University Press |location=Oxford [Oxfordshire] |year=2006 |isbn=0-19-856585-2 |series=Oxford Series in Ecology and Evolution}}</ref> This idea can probably be more readily seen if we say that the typical time between contacts is <math>T_{c} = \beta^{-1}</math>, and the typical time until removal is <math>T_{r} = \gamma^{-1}</math>. From here it follows that, on average, the number of contacts by an infectious individual with others ''before'' the infectious has been removed is: <math>T_{r}/T_{c}.</math>

By dividing the first differential equation by the third, [[Separation of variables|separating the variables]] and integrating we get

:<math> S(t) = S(0) e^{-R_0(R(t) - R(0))/N}, </math>

where <math>S(0)</math> and <math>R(0)</math> are the initial numbers of, respectively, susceptible and removed subjects. 
Writing <math>s_0 = S(0) / N</math> for the initial proportion of susceptible individuals, and
<math>s_\infty = S(\infty) / N</math> and 
<math>r_\infty = R(\infty) / N</math> for the proportion of susceptible and removed individuals respectively in the limit <math>t \to \infty,</math> one has

:<math>s_\infty = 1 - r_\infty = s_0 e^{-R_0(r_\infty - r_0)}</math>

(note that the infectious compartment empties in this limit).
This [[transcendental equation]] has a solution in terms of the [[Lambert W function|Lambert {{mvar|W}} function]],<ref>{{cite web |author1=Wolfram Research, Inc. |title=Mathematica, Version 12.1 |url=https://www.wolfram.com/mathematica |publisher=Champaign IL, 2020}}</ref> namely

:<math>s_\infty = 1-r_\infty = - R_0^{-1}\, W(-s_0 R_0 e^{-R_0(1-r_0)}).</math>

This shows that at the end of an epidemic that conforms to the simple assumptions of the SIR model, unless <math>s_0=0</math>, not all individuals of the population have been removed, so some must remain susceptible. A driving force leading to the end of an epidemic is a decline in the number of infectious individuals. The epidemic does not typically end because of a complete lack of susceptible individuals.

The role of both the basic reproduction number and the initial susceptibility are extremely important. In fact, upon rewriting the equation for infectious individuals as follows:

:<math> \frac{dI}{dt} = \left(R_0 \frac{S}{N}  - 1\right) \gamma I,</math>

it yields that if:

:<math> R_{0} \cdot S(0) > N,</math>

then:

:<math> \frac{dI}{dt}(0) >0 ,</math>

i.e., there will be a proper epidemic outbreak with an increase of the number of the infectious (which can reach a considerable fraction of the population). On the contrary, if

:<math> R_{0} \cdot S(0) < N,</math>

then

:<math> \frac{dI}{dt}(0) <0 ,</math>

i.e., independently from the initial size of the susceptible population the disease can never cause a proper epidemic outbreak. As a consequence, it is clear that both the basic reproduction number and the initial susceptibility are extremely important.

====The force of infection====

Note that in the above model the function:

:<math> F = \beta I,</math>

models the transition rate from the compartment of susceptible individuals to the compartment of infectious individuals, so that it is called the [[force of infection]]. However, for large classes of communicable diseases it is more realistic to consider a force of infection that does not depend on the absolute number of infectious subjects, but on their fraction (with respect to the total constant population <math>N</math>):

:<math> F = \beta \frac{I}{N} .</math>

Capasso<ref name="Capasso">{{cite book | vauthors = Capasso V |title=Mathematical Structure of Epidemic Systems |location=Berlin |publisher=Springer |year=1993 |isbn=3-540-56526-4 }}</ref> and, afterwards, other authors have proposed nonlinear forces of infection to model more realistically the contagion process.

====Exact analytical solutions to the SIR model====
In 2014, Harko and coauthors derived an exact so-called analytical solution (involving an integral that can only be calculated numerically) to the SIR model.<ref name="Harko" /> In the case  without vital dynamics setup, for <math>\mathcal{S}(u)=S(t)</math>, etc., it corresponds to the following time parametrization

:<math>\mathcal{S}(u)= S(0)u </math>
:<math>\mathcal{I}(u)= N -\mathcal{R}(u)-\mathcal{S}(u) </math>
:<math>\mathcal{R}(u)=R(0) -\rho \ln(u)</math>

for

:<math>t= \frac{N}{\beta}\int_u^1 \frac{du^*}{u^*\mathcal{I}(u^*)} , \quad \rho=\frac{\gamma N}{\beta},</math>

with initial conditions

:<math>(\mathcal{S}(1),\mathcal{I}(1),\mathcal{R}(1))=(S(0),N -R(0)-S(0),R(0)), \quad  u_T<u<1,</math>

where <math>u_T</math> satisfies <math>\mathcal{I}(u_T)=0</math>. By the transcendental equation for <math>R_{\infty}</math> above, it follows that <math>u_T=e^{-(R_{\infty}-R(0))/\rho}(=S_{\infty}/S(0)</math>, if <math>S(0) \neq 0)</math> and <math>I_{\infty}=0</math>.

An equivalent so-called analytical solution (involving an integral that can only be calculated numerically) found by Miller<ref>{{cite journal | vauthors = Miller JC | title = A note on the derivation of epidemic final sizes | journal = Bulletin of Mathematical Biology | volume = 74 | issue = 9 | pages = 2125–2141 | date = September 2012 | pmid = 22829179 | pmc = 3506030 | doi = 10.1007/s11538-012-9749-6 | quote = Section 4.1 }}</ref><ref>{{cite journal | vauthors = Miller JC | title = Mathematical models of SIR disease spread with combined non-sexual and sexual transmission routes | journal = Infectious Disease Modelling | volume = 2 | issue = 1 | pages = 35–55 | date = February 2017 | pmid = 29928728 | pmc = 5963332 | doi = 10.1016/j.idm.2016.12.003 | quote = Section 2.1.3 }}</ref> yields

:<math>
\begin{align}
S(t) & = S(0) e^{-\xi(t)} \\[8pt]
I(t) & = N-S(t)-R(t) \\[8pt]
R(t) & = R(0) + \rho \xi(t) \\[8pt]
\xi(t) & = \frac{\beta}{N}\int_0^t I(t^*) \, dt^*
\end{align}
</math>

Here <math>\xi(t)</math> can be interpreted as the expected number of transmissions an individual has received by time <math>t</math>.  The two solutions are related by <math>e^{-\xi(t)} = u</math>.

Effectively the same result can be found in the original work by Kermack and McKendrick.<ref name="Kermack–McKendrick">{{cite journal |vauthors=Kermack WO, McKendrick AG |date=1927 |title=A Contribution to the Mathematical Theory of Epidemics |journal=Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character |volume=115 |issue=772 |pages=700–721 |bibcode=1927RSPSA.115..700K |doi=10.1098/rspa.1927.0118 |doi-access=free}}</ref>

These solutions may be easily understood by noting that all of the terms on the right-hand sides of the original differential equations are proportional to <math>I</math>. The equations may thus be divided through by <math>I</math>, and the time rescaled so that the differential operator on the left-hand side becomes simply <math>d/d\tau</math>, where <math>d\tau=I dt</math>, i.e. <math>\tau=\int I dt</math>. The differential equations are now all linear, and the third equation, of the form <math>dR/d\tau =</math> const., shows that <math>\tau</math> and <math>R</math> (and <math>\xi</math> above) are simply linearly related.

A highly accurate analytic approximant of the SIR model as well as exact analytic expressions for the final values <math>S_{\infty}</math>, <math>I_{\infty}</math>, and <math>R_{\infty}</math> were provided by [[Martin Kröger|Kröger]] and Schlickeiser,<ref name="KrogerSchlickeiser"/> so that there is no need to perform a numerical integration to solve the SIR model (a simplified example practice on [[COVID-19]] numerical simulation using [[Microsoft Excel]] can be found here <ref name="SIR_COVID19">{{cite journal | vauthors = Hart KD, Thompson C, Burger C, Hardwick D, Michaud AH, Bulushi A, Pridemore C, Ward C, Chen J |title = Remote Learning of COVID-19 Kinetic Analysis in a Physical Chemistry Laboratory Class | journal = ACS Omega | volume = 6 | issue= 43 | pages = 29223–29232 |year=2021 |doi = 10.1021/acsomega.1c04842|pmid = 34723043 |pmc = 8547164 }}</ref>), to obtain its parameters from existing data, or to predict the future dynamics of an epidemics modeled by the SIR model. The approximant involves the [[Lambert W function|Lambert {{mvar|W}} function]] which is part of all basic data visualization software such as Microsoft Excel, [[MATLAB]], and [[Wolfram Mathematica|Mathematica]].

While [[David George Kendall|Kendall]]<ref name="Kendall">{{cite book |title=Contributions to Biology and Problems of Health |vauthors=Kendall DG |date=1956 |publisher=University of California Press |volume=4 |pages=149–165 |chapter=Deterministic and Stochastic Epidemics in Closed Populations |doi=10.1525/9780520350717-011 |mr=0084936 |zbl=0070.15101 |chapter-url=http://projecteuclid.org/euclid.bsmsp/1200502553}}</ref> considered the so-called all-time SIR model where the initial conditions <math>S(0)</math>, <math>I(0)</math>, and <math>R(0)</math> are coupled through the above relations, Kermack and McKendrick<ref name="Kermack–McKendrick"/> proposed to study the more general semi-time case, for which <math>S(0)</math> and <math>I(0)</math> are both arbitrary. This latter version, denoted as semi-time SIR model,<ref name="KrogerSchlickeiser"/> makes predictions only for future times <math>t>0</math>. An analytic approximant and exact expressions for the final values are available for the semi-time SIR model as well.<ref name="KrogerSchlickeiser_partB"/>

====Numerical solutions to the SIR model with approximations====
Numerical solutions to the SIR model can be found in the literature. An example is using the model to analyze [[COVID-19]] spreading data.<ref name="SIR_COVID19"/><ref name="SIRV_COVID19">{{Cite journal| vauthors = Smith DK, Lauro K, Kelly D, Fish J, Lintelman E, McEwen D, Smith C, Stecz M, Ambagaspitiya TD, Chen J | title=Teaching Undergraduate Physical Chemistry Lab with Kinetic Analysis of COVID-19 in the United States |journal= Journal of Chemical Education |year=2022| volume=99 | issue=10 | pages=3471–3477 |doi=10.1021/acs.jchemed.2c00416| pmid=36589277 | pmc=9799982 | bibcode=2022JChEd..99.3471S | s2cid=251484056 }}</ref> Three reproduction numbers can be pulled out from the data analyzed with numerical approximation,
:the [[basic reproduction number]]:

::<math>R_0=\frac{\beta_0}{\gamma_0}</math>

:the real-time reproduction number:

::<math>R_t=\frac{\beta_t}{\gamma_t}</math>

:and the real-time effective reproduction number:

::<math>R_e=\frac{\beta_tS}{\gamma_tN}</math>

<math>R_0</math> represents the speed of reproduction rate at the beginning of the spreading when all populations are assumed susceptible, e.g. if <math>\beta_0 = 0.4  day^{-1}</math> and <math>\gamma_0 = 0.2  day^{-1}</math> meaning one infectious person on average infects 0.4 susceptible people per day and recovers in 1/0.2=5 days. Thus when this person recovered, there are two people still infectious directly got from this person and <math>R_0 = 2</math>, i.e. the number of infectious people doubled in one cycle of 5 days. The data simulated by the model with <math>R_0 = 2</math> or real data fitted will yield a doubling of the number of infectious people faster than 5 days because the two infected people are infecting people. From the SIR model, we can tell that <math>\beta</math> is determined by the nature of the disease and also a function of the interactive frequency between the infectious person <math>I</math> with the susceptible people <math>S</math> and also the intensity/duration of the interaction like how close they interact for how long and whether or not they both wear masks, thus, it changes over time when the average behavior of the carriers and susceptible people changes. The model use <math>SI</math> to represent these factors but it indeed is referenced to the initial stage when no action is taken to prevent the spread and all population is susceptible, thus all changes are absorbed by the change of <math>\beta</math>.

<math>\gamma</math> is usually more stable over time assuming when the infectious person shows symptoms, she/he will seek medical attention or be self-isolated. So if we find <math>R_t</math> changes, most probably the behaviors of people in the community have changed from their normal patterns before the outbreak, or the disease has mutated to a new form. Costive massive detection and isolation of susceptible close contacts have effects on reducing <math>1/\gamma</math> but whose efficiencies are under debate. This debate is largely on the uncertainty of the number of days reduced from after infectious or detectable whichever comes first to before a symptom shows up for an infected susceptible person. If the person is infectious after symptoms show up, or detection only works for a person with symptoms, then these prevention methods are not necessary, and self-isolation and/or medical attention is the best way to cut the <math>1/\gamma</math> values. The typical onset of the [[COVID-19]] infectious period is in the order of one day from the symptoms showing up, making massive detection with typical frequency in a few days useless.

<math>R_t</math> does not tell us whether or not the spreading will speed up or slow down in the latter stages when the fraction of susceptible people in the community has dropped significantly after recovery or vaccination. <math>R_e</math> corrects this dilution effect by multiplying the fraction of the susceptible population over the total population. It corrects the effective/transmissible interaction between an infectious person and the rest of the community when many of the interaction is immune in the middle to late stages of the disease spreading. Thus, when <math>R_e > 1</math>, we will see an exponential-like outbreak; when <math>R_e = 1</math>, a steady state reached and no number of infectious people changes over time; and when <math>R_e < 1</math>, the disease decays and fades away over time.

Using the differential equations of the SIR model and converting them to numerical discrete forms, one can set up the recursive equations and calculate the S, I, and R populations with any given initial conditions but accumulate errors over a long calculation time from the reference point. Sometimes a [[convergence test]] is needed to estimate the errors. Given a set of initial conditions and the disease-spreading data, one can also fit the data with the SIR model and pull out the three reproduction numbers when the errors are usually negligible due to the short time step from the reference point.<ref name="SIR_COVID19"/><ref name="SIRV_COVID19"/> Any point of the time can be used as the initial condition to predict the future after it using this numerical model with assumption of time-evolved parameters such as population, <math>R_t</math>, and <math>\gamma</math>. However, away from this reference point, errors will accumulate over time thus [[convergence test]] is needed to find an optimal time step for more accurate results.

Among these three reproduction numbers, <math>R_0</math> is very useful to judge the control pressure, e.g., a large value meaning the disease will spread very fast and is very difficult to control. <math>R_t</math> is most useful in predicting future trends, for example, if we know the social interactions have reduced 50% frequently from that before the outbreak and the interaction intensities among people are the same, then we can set <math>R_t = 0.5R_0</math>. If social distancing and masks add another 50% cut in infection efficiency, we can set <math>R_t = 0.25R_0</math>. <math>R_e</math> will perfectly correlate with the waves of the spreading and whenever <math>R_e>1</math>, the spreading accelerates, and when <math>R_e<1</math>, the spreading slows down thus useful to set a prediction on the short-term trends. Also, it can be used to directly calculate the threshold population of vaccination/immunization for the [[herd immunity]] stage by setting <math>R_t = R_0</math>, and <math>R_E = 1</math>, i.e. <math>S = N/R_0</math>.