A quick overview of control measures in times of pandemic/epidemic
modeling
SEIR
epidemiology
stochastic
COVID-19
real-world
Author
Jeffrey Post
Published
April 2, 2020
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:
Reaching herd immunity all while controlling peak infectious individuals to within hospital and sanitary resources so as to limit morbidity and mortality
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.
#collapse_hide!pip install plotly==4.8.2import pandas as pdimport numpy as npimport plotly.graph_objects as goimport plotly.express as pxfrom scipy.stats import exponfrom scipy.stats import gammafrom scipy.stats import weibull_minfrom numpy.random import default_rngrng = default_rng()
#collapse_hide# 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
#collapse_hide# 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 3if isolation isTrue: 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 runfor j inrange(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 valuesif lockdown isTrue:if ((I[-1] >100) & (Ipd[-1] >39)):if lockdown_date ==0: lockdown_date = j+1 b = beta2elif ((I[-1] >100) & (Ipd[-1] <40)): b = beta3 Epd[0]+=num_E Ipd[0]+=num_I Rpd[0]+=num_Rreturn 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.
#collapse_hide# Define parameters for stochastic modeldays =300p =10000num_E =1num_I =0# Run 2 simulations, one above HIT, and one below:num_R1 =8000num_R2 =9200beta_stoch =0.5*np.ones(days)# Run 4 stochastic simulations, 2 with exponential sigma, 2 with gamma sigmalation, iso_n, contact_tracing, lockdownresults_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_n1, contact_tracing=False, lockdown=False)
NameError: ignored
#collapse_hidefig = 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