Model Description

We use a compartmental epidemiological model, based on the classic SEIR model, to describe the spread and clinical progression of COVID-19. A nice primer to this sort of model is available on Wikipedia. It is important to track the different clinical outcomes of infection, since they require different level of healthcare resources to care for and may be tested and isolated at different rates. Susceptible (\(S\)) individuals who become infected start out in an exposed class \(E\), where they are asymptomatic and do not transmit infection. The rate of progressing from the exposed stage to the infected stage \(I\), where the individual is symptomatic and infectious, occurs at rate \(a\). The clinical descriptions of the different stages of infection are given below. Infected individuals begin with mild infection (\(I_1\)), from which they either recover, at rate \(\gamma_1\), or progress to severe infection (\(I_2\)), at rate \(p_1\). Severe infection resolves at rate \(\gamma_2\) or progresses to a critical stage (\(I_3\)) at rate \(p_2\). Individuals with critical infection recover at rate \(\gamma_3\) and die at rate \(\mu\). Recovered individuals are tracked by class \(R\) and are assumed to be protected from re-infection for life. Individuals may transmit the infection at any stage, though with different rates. The transmission rate in stage \(i\) is described by \(\beta_i\) .

Equations

\begin{equation} \dot{S} = -(\beta_1 I_1 -\beta_2 I_2 - \beta_3 I_3) S \end{equation}

\begin{equation} \dot{E} =(\beta_1 I_1 +\beta_2 I_2 + \beta_3 I_3) S - a E \ \end{equation}

\begin{equation} \dot{I_1} = a E - (\gamma_1 + p_1) I_1 \ \end{equation}

\begin{equation} \dot{I_2} = p_1 I_1 -(\gamma_2 + p_2) I_2 \ \end{equation}

\begin{equation} \dot{I_3} = p_2 I_2 -(\gamma_3 + \mu) I_3 \ \end{equation}

\begin{equation} \dot{R} = \gamma_1 I_1 + \gamma_2 I_2 + \gamma_3 I_3 \ \end{equation}

\begin{equation} \dot{D} = \mu I_3 \end{equation}

Variables

\(S\): Susceptible individuals

\(E\): Exposed individuals - infected but not yet infectious or symptomatic

\(I_i\): Infected individuals in severity class \(i\). Severity increases with \(i\) and we assume individuals must pass through all previous classes \(I_1\): Mild infection \(I_2\): Severe infection \(I_3\): Critical infection

\(R\): individuals who have recovered from disease and are now immune

\(D\): Dead individuals

\(N=S+E+I_1+I_2+I_3+R+D\) Total population size (constant)

Parameters

\(\beta_i\) rate at which infected individuals in class \(I_i\) contact susceptibles and infect them

\(a\) rate of progression from the exposed to infected class

\(\gamma_i\) rate at which infected individuals in class \(I_i\) recover from disease and become immune

\(p_i\) rate at which infected individuals in class \(I_i\) progress to class \(I_{i+1}\)

\(\mu\) death rate for individuals in the most severe stage of disease

All rates are per day.

Clinical stages

Mild infection - These individuals have symptoms like fever and cough and may have mild pneumonia. Hospitalization is not required (though in many countries such individuals are also hospitalized)

Severe infection - These individuals have more severe pneumonia that leads to dyspnea, respiratory frequency <30/min, blood oxygen saturation <93%, partial pressure of arterial oxygen to fraction of inspired oxygen ratio <300, and/or lung infiltrates >50% within 24 to 48 hours. Hospitalization and supplemental oxygen are generally required.

Critical infection - These individuals experience respiratory failure, septic shock, and/or multiple organ dysfunction or failure. Treatment in an ICU, often with mechanical ventilation, is required.

Relating clinical observations to model parameters

To determine the model parameters consistent with current clinical data, we collect the following values from the slider values chosen by the user, and then use the formulas below to relate them to the rate parameters in the model. Note that the slider inputs for time intervals are average durations.

IncubPeriod: Average incubation period, days

DurMildInf: Average duration of mild infections, days

FracMild: Average fraction of (symptomatic) infections that are mild

FracSevere: Average fraction of (symptomatic) infections that are severe

FracCritical: Average fraction of (symptomatic) infections that are critical

CFR: Case fatality rate (fraction of infections that eventually result in death)

DurHosp: Average duration of hospitalization (time to recovery) for individuals with severe infection, days

TimeICUDeath: Average duration of ICU admission (until death or recovery), days

(Note g=\(\gamma\))

a=1/IncubPeriod g1=(1/DurMildInf)*FracMild p1=(1/DurMildInf)-g1 p2=(1/DurHosp)*(FracCritical/(FracSevere+FracCritical)) g2=(1/DurHosp)-p2 u=(1/TimeICUDeath)*(CFR/FracCritical) g3=(1/TimeICUDeath)-u

Basic reproductive ratio

Idea: \(R_0\) is the sum of

the average number of secondary infections generated from an individual in stage \(I_1\) the probability that an infected individual progresses to \(I_2\) multiplied by the average number of secondary infections generated from an individual in stage \(I_2\) the probability that an infected individual progresses to \(I_3\) multiplied by the average number of secondary infections generated from an individual in stage \(I_3\)

\begin{equation} R_0 = N\frac{\beta_1}{p_1+\gamma_1} + \frac{p_1}{p_1 + \gamma_1} \left( \frac{N \beta_2}{p_2+\gamma_2} + \frac{p_2}{p_2 + \gamma_2} \frac{N \beta_3}{\mu+\gamma_3}\right) \end{equation}

\begin{equation} = N\frac{1}{p_1+\gamma_1} \left(\beta_1 + \frac{p_1}{p_2 + \gamma_2} \left( \beta_2 + \beta_3 \frac{p_2}{\mu + \gamma_3} \right) \right) \end{equation}

Calculations using the next generation matrix give the same results.

Early epidemic growth rate

Early in the epidemic, before susceptibles are depleted, the epidemic grows at an exponential rate \(r\), which can also be described with doubling time T\(_2\)=ln(2)\(/r\). During this phase all infected classes grow at the same rate as each other and as the deaths and recovered individuals. The cumulative number of infections that have happened since the outbreak started also grows at the same rate. This rate can be calculated from the dominant eigenvalue of the linearized system of equations in the limit that \(S=N\).

During this early exponential growth phase, there will be a fixed ratio of individuals between any pair of compartments. This expected ratio could be used to estimate the amount of underreporting in data. For example, we might think that all deaths are reported, but that some mild infections might not be reported, since these patients might not seek healthcare or might not be prioritized for testing. These ratios have expected values under the model for a fixed set of parameters. They can be calculated by finding the eigenvector corresponding to the dominant eigenvalue (\(r\)) for the linearized system described above. Ratios that deviate from these values suggest either a) underreporting of cases relative to deaths, or b) local differences in the clinical parameters of disease progression. The expected ratios are

\begin{equation} \frac{I_3}{D} = \frac{r}{\mu} \end{equation}

\begin{equation} \frac{I_2}{D} = \frac{(\mu+\gamma_3+r)}{p_2}\frac{r}{\mu} \end{equation}

\begin{equation} \frac{I_1}{D} = \frac{(p_2+\gamma_2+r)}{p_1}\frac{(\mu+\gamma_3+r)}{p_2}\frac{r}{\mu} \end{equation}

\begin{equation} \frac{\textrm{Total symptomatic}}{D} =_ \sum I_i = \frac{r}{\mu} \left[1 + \frac{(\mu+\gamma_3+r)}{p_2} \left(1+\frac{(p_2+\gamma_2+r)}{p_1} \right) \right] \end{equation}

\begin{equation} \frac{E}{D} = \frac{(p_1+\gamma_1+r)}{a}\frac{(p_2+\gamma_2+r)}{p_1}\frac{(\mu+\gamma_3+r)}{p_2}\frac{r}{\mu} \end{equation}

Assumptions

This model is formulated as a system of differential equations and the output therefore represents the expected values of each quantity. It does not take into account stochastic events, and so the epidemic cannot go extinct even when it gets to very low values (except when an intervention is stopped, at which case the number of individuals in each state is rounded to the nearest integar). The model does not report the expected variance in the variables, which can sometimes be large.

Individuals must pass through a mild stage before reaching a severe or critical stage

Individuals must pass through a severe stage before reaching a critical stage

Only individuals in a critical stage die

All individuals have equal transmission rates and equal susceptiblity to infection

Updates

Mar 21 2020

The model now includes the possibility for asymptomatic infection. After leaving the \(E\) class, a fraction \(f\) of individuals develop asymptomatic infection (enter \(I_0\) class), whereas the remaining fraction \(1-f\) develop symptomatic infection (enter \(I_1\) class). Asymptomatic infection never progresses to more severe stages. The rate of recovery from asymptomatic infection is \(\gamma_0\). Asymptomatically-infected individuals may transmit to others at rate \(\beta_0\). The original sliders that control the fractions of infections that are mild vs severe vs critical now have the interpretation as being the fraction of symptomatic infections that enter each of these stages.

The model now also includes the possibility that exposed individuals who have not yet developed symptoms may still be able to transmit the virus (“presymptomatic transmission”). To model this, we divide the \(E\) class into two separate classes, \(E_0\) (no symptoms or transmission) and \(E_1\) (no symptoms but can transmit). The rate of exit from \(E_0\) is \(a_0\) and from \(E_1\) is \(a_1\).

We now include the option for seasonality in the tranmission rates. All transmission rates are modified by a factor \(\sigma(t) = 1 + \epsilon \cos(2 \pi (t-\phi))\) where \(\epsilon \in [0,1]\) is the relative amplitude of the seasonal oscillations and and \(\phi \in [-\infty, \infty]\) is the phase, and determines the time (in years) of the peak in transmisison relative to the time the simulation starts. The values the user inputs for the transmisison rates are interpreted as the rates at time zero of the simulation. This input will be equal to the peak transmission if \(\phi = 0\), as the minimum transmission of if \(\phi=365/4 \sim 90\), and as the time-average transmission if \(\phi=365/2 \sim 180\), for example.

The updated model equations are

\begin{equation} \dot{S} = -(\beta_e E_1 + \beta_0 I_0 + \beta_1 I_1 +\beta_2 I_2 + \beta_3 I_3) S \sigma(t) \end{equation}

\begin{equation} \dot{E_0} =(\beta_e E_1 + \beta_0 I_0 +\beta_1 I_1 +\beta_2 I_2 + \beta_3 I_3) S \sigma(t) - a_0 E_0 \ \end{equation}

\begin{equation} \dot{E_1} = a_0 E_0 - a_1 E \ \end{equation}

\begin{equation} \dot{I_0} = f a_1 E_1 - \gamma_0 I_0 \ \end{equation}

\begin{equation} \dot{I_1} = (1-f) a_1 E_1 - (\gamma_1 + p_1) I_1 \ \end{equation}

\begin{equation} \dot{I_2} = p_1 I_1 -(\gamma_2 + p_2) I_2 \ \end{equation}

\begin{equation} \dot{I_3} = p_2 I_2 -(\gamma_3 + \mu) I_3 \ \end{equation}

\begin{equation} \dot{R} = \gamma_0 I_0 + \gamma_1 I_1 + \gamma_2 I_2 + \gamma_3 I_3 \ \end{equation}

\begin{equation} \dot{D} = \mu I_3 \end{equation}

The extra slider inputs are

FracAsym: Fraction of all infections that are asymptomatic

PresymPeriod: Length of infectious phase of incubation period

DurAsym: Duration of asympatomatic infection

And the formula for extracting the rate constants from these inputs are

a1=1/PresymPeriod a0=(IncubPeriod-PresymPeriod)^(-1) f=FracAsym g0=1/DurAsym

The basic reproductive ratio becomes

\begin{equation} R_0 = N \left[ \frac{\beta_e}{a_1} + f \frac{\beta_0}{\gamma_0} + (1-f) \frac{1}{p_1+\gamma_1} \left(\beta_1 + \frac{p_1}{p_2 + \gamma_2} \left( \beta_2 + \beta_3 \frac{p_2}{\mu + \gamma_3} \right) \right) \right] \end{equation}

Output

Rate parameters of dynamic model

These parameters can be changed using the sliders in the other tabs. The values in this table represent the current values chosen via the sliders. Note that the transmission rates chosen by the sliders are always scaled by \(N\), so that \(\beta*N\) is constant as \(N\) changes.

Ratios of cases during early growth phase

These values are calculated based on the current model parameters