Epidemic modeling - Part 6

Motivation for write-up

This is the 6th part of a multi-part series blog post on modeling in epidemiology.

The COVID-19 pandemic has brought a lot of attention to the study of epidemiology and more specifically to the various mathematical models that are used to inform public health policies. Everyone has been trying to understand the growth or slowing of new cases and trying to predict the necessary sanitary resources. This blog post attempts to explain the foundations for some of the most used models and enlighten the reader on two key points.

Following the first 5 parts of this blog series, we are left wondering what possible measures can be put in place to control the epidemic.

This 6th installment focuses on this and attempts to elucidate the subject.

Why do we need control measures during an epidemic?

In the previous sections we have seen that without control measures in place i.e. when injecting an exposed person into the simulation, depending on the values of \(\beta\) and \(\gamma\), the virus will spread until it has infected everyone.

There are many reasons why this is bad with a virus as virulent as SARS-CoV-2:

• Even with a low Case Fatality Rate (CFR), the total death toll will be unacceptable
• The strain on sanitary resources of peak sick people will lead to an increase in all-cause deaths
• The strain on sanitary resources will lead to increased long-term morbidity
• Many economic implications
• Many many others…

Control measures: the basics

While the previous sections might have left us wondering what we can possibly do to control the pread of the epidemic, the evidence left from those studies in fact leave us with a lot of clues as to where to begin.

Only one way to stop an epidemic

We have seen there is only one way to stop an epidemic, and that is by having:

\[R \leq 1\] \[\leftrightarrow R_0~s(t) \leq 1\] \[\leftrightarrow \frac{\beta}{\gamma}~s(t) \leq 1\]

Two methods to save lives

From the equation above, we see an epidemic will either continue until herd immunity is reached:

\[s(t) \leq \frac{1}{R_0}\] \[\leftrightarrow Immune(t) \geq 1-\frac{1}{R_0}\]

Or until measures are in place so that:

\[R_0 \leq \frac{1}{s(t)}\]

In the worst case scenario where \(s(t) = 1\) (completely susceptible population), we need:

\[R_0 \leq 1\] \[\leftrightarrow \beta \leq \gamma\]

From this results two main ideas of control:

1. Reaching herd immunity all while controlling peak infectious individuals to within hospital and sanitary resources so as to limit morbidity and mortality

2. Reducing \(\beta\) to smaller than \(\gamma\) to stop the epidemic even before herd immunity

While the difference between the two is noted here, in practice they are effectively the same.

Option 1 is the same as option 2, with the difference that in 1 we don’t quite have \[R\leq 1\] and so the difference between the two results from different levels of imlementation.

Control measures in practice

We have seen three things influence the total number of people that end up infected with the virus:

• Proportion of Susceptible in population \(s(t)\)
• Value \(\beta\)
• Value of \(\gamma\) (\(T_{Infectious}\))

We have also seen that the peak of infectious individuals could be affected by:

• Value of \(\sigma\) (\(T_{Latent}\))

And indeed, reducing the peak of infectious at any point comes down to either:

• Reducing the number of S:
• vaccination
• prophylactic treatment when potentially exposed
• Reducing \(\beta = r * \rho\) by:

1. reducing \(r\) i.e. reducing the average number of contacts a person has per day:

• lockdown measures
• work from home
• closing places where people gather (restaurants, bars, places of worship, etc)

2. reducing \(\rho\) - reducing the pobability of transmitting infection from an infectious to a susceptible via:

• physical distancing
• hygiene measures
• wearing personal protective equipment (PPE i.e. masks, gloves, etc)
• Reducing \(\gamma\) by:
• Isolation of the sick (mass testing of symptomatic and immediate isolation of positive cases)
• Contact-tracing: tracing and quarantining all contacts from infectious people as to quarantine potential exposed before they become infectious
• Chemotherapy: treatment to shorten duration of sickness and infectious period
• Increasing \(\sigma\) by:
• prophylactic treatment when potentially exposed

All these are the same requirements to reduce \(R \leq 1\).

Let’s try to quantify the impact of each measure below.

Code

import pandas as pd
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
from scipy.stats import expon
from scipy.stats import gamma
from scipy.stats import weibull_min
from numpy.random import default_rng
rng = default_rng()
import plotly
import plotly.io as pio
from IPython.display import display, HTML

## Tomas Mazak's workaround
plotly.offline.init_notebook_mode()
display(HTML(
    #'<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_SVG"></script>'
    '<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_SVG"></script>'
    '<script src="https://cdn.plot.ly/plotly-3.0.1.js" charset="utf-8"></script>'
))

pio.renderers.default = "plotly_mimetype+notebook_connected"

pio.templates.default = "plotly_dark"

Code

# Need this new function for model below:
def make_df(p,num_E, num_I, num_R):
  df = pd.DataFrame(np.full((p,1), 'S').T[0], columns=['State'])
  df['Day'] = 0
  tochange=df.loc[rng.choice(p, size=num_E+num_I+num_R, replace=False),'State'].index
  df.loc[tochange[0:num_E],'State'] = 'E'
  df.loc[tochange[num_E:num_I+num_E],'State'] = 'I'
  df.loc[tochange[num_E+num_I:num_E+num_I+num_R],'State'] = 'R'
  return df

Code

# regular Stoachastic SEIR model below:
def seir_model_stoch_ctrl(beta, p, num_E, num_I, num_R, days, isolation, iso_n, contact_tracing, lockdown):

    # Initialize population dataframe with data given by user
    df = make_df(p,num_E, num_I, num_R)
    
    # This variable is used to track daily value of beta
    xxbeta=np.array([],dtype=float)

    # Initialize the arrays to return
    # Below are numbers of S, E, I, R total
    S=np.array([],dtype=int)
    E=np.array([],dtype=int)
    I=np.array([],dtype=int)
    R=np.array([],dtype=int)
    # Below are the daily additions in S, E, I, R
    Spd=np.array([],dtype=int)
    Epd=np.array([],dtype=int)
    Ipd=np.array([],dtype=int)
    Rpd=np.array([],dtype=int)

    b=beta
    beta2=b/10
    beta3=b/10

    lockdown_date=0
    
    # Stochastic model so use random values to decide on progression
    rand = np.random.random(size=(p,days))

    # Depending if you want exponential or gamma distribution for sigma
    EtoI = gamma.rvs(1.8,loc=0.9,scale=(5.2-1.8)/0.9,size=p)

    # Depending if you want exponential or gamma distribution for gamma and if you have isolation or not
    # Uses distributiosn found on blog part 3
    if isolation is True:
      ItoR = iso_n*np.ones(p)
    else:
      ItoR = weibull_min.rvs(2.3, loc=2, scale=20.11, size=p)   


    # Iterate over every day the simulation is run
    for j in range(0,days):

        # Record daily beta values
        xxbeta=np.append(xxbeta, b[j])

        # First we get the index of the individuals that will change state today:

        # Random number tells you which 'S' have been exposed on this day 
        StoE_index = df.loc[(df.State == 'S') & (rand[:,j] < b[j]*len(np.where(df.State=='I')[0])/p)].index

        # For each row, if a person has been a certain number of days in E, they will go to I
        # This follows EtoI variable which is either exponential or gamma distributed according to above
        EtoI_index = df.loc[(df.State == 'E') & (j-df.Day >= EtoI)].index
        
        # Similaraly as above
        # For each row, if a person has been a certain number of days in I, they will go to R
        # This follows EtoI variable which is either exponential or gamma distributed according to above
        ItoR_index = df.loc[(df.State == 'I') & (j-df.Day >= ItoR)].index

        # Use indexes collected above to populate per day values
        Epd = np.append(Epd,len(StoE_index))
        Ipd = np.append(Ipd,len(EtoI_index))
        Rpd = np.append(Rpd,len(ItoR_index))

        # Append the S, E, I, and R arrays
        S=np.append(S,len(np.where(df.State=='S')[0]))
        E=np.append(E,len(np.where(df.State=='E')[0]))
        I=np.append(I,len(np.where(df.State=='I')[0]))
        R=np.append(R,len(np.where(df.State=='R')[0]))

        # Now we use the indexes collected above randomly to change the actual population dataframe to the new states
        df.iloc[ItoR_index] = ['R', j]
        df.iloc[EtoI_index] = ['I', j]
        df.iloc[StoE_index] = ['E', j]

        # Code below for control measures to reduce beta values
        if lockdown is True:
          if ((I[-1] > 100) & (Ipd[-1] > 39)):
            if lockdown_date == 0:
              lockdown_date = j+1
            b = beta2
          elif ((I[-1] > 100) & (Ipd[-1] < 40)): 
            b = beta3
        
    Epd[0]+=num_E
    Ipd[0]+=num_I
    Rpd[0]+=num_R

    return S,E,I,R, Epd, Ipd, Rpd, xxbeta, lockdown_date

Reducing the proportion of Susceptible in the population

As we have seen, reducing the proportion of susceptible in the population helps reduce the impact of the epidemic.

In the first blog post of the series, we derived the threshold for herd immunity as being:

\[1-\frac{1}{R_0}\]

In our simulations we have:

\[R_0 = \frac{\beta}{\gamma} = \frac{0.5}{\frac{1}{20.62}} = 10.31\]

And so, the herd immunity threshold in our simulation should be:

\[HIT = \frac{9.31}{10.31}~100\% = 90.3\%\]

Let’s run the model with the initial condition that 92% are in the R state already.

Code

# Define parameters for stochastic model
days = 300
p = 10000
num_E = 1
num_I = 0

# Run 2 simulations, one above HIT, and one below:
num_R1 = 8000
num_R2 = 9200
beta_stoch = 0.5*np.ones(days)
iso_n = 0

# Run 4 stochastic simulations, 2 with exponential sigma, 2 with gamma sigmalation, iso_n, contact_tracing, lockdown
results_stoch0 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R1, days, isolation=False, iso_n=iso_n, contact_tracing=False, lockdown=False)
results_stoch1 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R2, days, isolation=False, iso_n=iso_n, contact_tracing=False, lockdown=False)

Code

fig = go.Figure(data=[    
    go.Scatter(name='I_below_HIT', x=np.arange(len(results_stoch0[0])), y=results_stoch0[2]/p),
    go.Scatter(name='I_above_HIT', x=np.arange(len(results_stoch1[0])), y=results_stoch1[2]/p),
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
    title={
        'text':r'$\text{Effect of herd immunity on SEIR model}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

Reducing \(\beta\)

See blog post 2 for effect of \(\beta\) on the SEIR model.

Reducing \(\gamma\)

Isolating positive tests

Say you isolate positive tests and are able to test everyone.

What if tests are positive \(n\) days after people are infectious ?

You can isolate people after \(n\) days and this effectively reduces \(T_{Infectious}\) to \(n\) days.

Let’s plot the total infectious individuals with:

• n = 10 days
• n = 7 days
• n = 5 days
• n = 2 days

Code

# Define parameters for stochastic model
days = 200
p = 10000
num_E = 1
num_I = 0
num_R = 0
beta_stoch = 0.5*np.ones(days)

# Isolate after iso_n days
iso_n1 = 10 # Test and isolate all infectious after 10 days
iso_n2 = 7
iso_n3 = 5
iso_n4 = 2

# Run 4 stochastic simulations, 2 with exponential sigma, 2 with gamma sigmalation, iso_n, contact_tracing, lockdown
results_stoch0 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R, days, isolation=False, iso_n=iso_n, contact_tracing=False, lockdown=False)
results_stoch1 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R, days, isolation=True, iso_n=iso_n1, contact_tracing=False, lockdown=False)
results_stoch2 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R, days, isolation=True, iso_n=iso_n2, contact_tracing=False, lockdown=False)
results_stoch3 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R, days, isolation=True, iso_n=iso_n3, contact_tracing=False, lockdown=False)
results_stoch4 = seir_model_stoch_ctrl(beta_stoch, p, num_E, num_I, num_R, days, isolation=True, iso_n=iso_n4, contact_tracing=False, lockdown=False)

Code

fig = go.Figure(data=[    
    go.Scatter(name='I_no_isolation', x=np.arange(len(results_stoch0[0])), y=results_stoch0[2]/p),
    go.Scatter(name='I_iso10', x=np.arange(len(results_stoch1[0])), y=results_stoch1[2]/p),
    go.Scatter(name='I_iso7', x=np.arange(len(results_stoch2[0])), y=results_stoch2[2]/p),
    go.Scatter(name='I_iso5', x=np.arange(len(results_stoch3[0])), y=results_stoch3[2]/p),
    go.Scatter(name='I_iso2', x=np.arange(len(results_stoch4[0])), y=results_stoch4[2]/p),
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
    title={
        'text':'Effect of early testing and isolation on SEIR model',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

Clearly, early testing and isolation of cases reduces transmission and this the peak of infectious individuals. The faster we do this the better.

If we catch infectious individuals early enough, we can completely control the epidemic (see I_iso_2 above).

Discussion

This was a quick introduction to the basics of control measures and their impact on the SEIR model.

While there are many different possible approaches, the accumulation of the different solutions is the way to achieve control over an epidemic.

Maybe I will make an interactive app to see the effect of each effect.