Building a new stochastic SEIR model to deal with probability distributions
probability distributions
modeling
SEIR
epidemiology
stochastic
HIV
AIDS
Author
Jeffrey Post
Published
May 5, 2020
Building a stochastic model
This model is going to focus on individual level disease dynamics to model the disease propagation.
It models DHS dataset which contains a homogeneous population between 15 and 49 years old.
The basic idea of this model is to have a dataframe with the number of rows equal to the population size (each individual is a row) and two columns: * State column to describe the state of each individual (S, I, or D) * Year column to save the day of transition of the individual into that state * Age column to know the age of the individuals
However, the population-level rates of transmission still apply here i.e. a person goes from S → I following two points: 1. the effective contact rate \(\beta\), which is itself given by: - the number of contacts the person has per unit time (given by \(r\)) - the chance of an S contracting the disease from a contact with an I (given by \(\rho\)) 2. the chance a given contact is with an I - infectious individual (the higher the number of I, the higher the chance)
This is done stochastically.
Once a person becomes I, their progression is unique to them. This progression is calculated in advance for computational reason, but it allows to use the time ditributions we want.
#collapse_hide!pip install plotly==4.14.3import pandas as pdimport numpy as npimport mathimport 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()import tqdmimport time
Below is a function to create the initial population dataframe: * \(p\) is the population number * \(num_I\) is the number of infectious on day 0 * \(num_R\) is the number of people recovered on day 0
#collapse_hide# Need this new function for model below:def make_df(p, num_I, num_R): df = pd.DataFrame(np.full((p,1), 'S').T[0], columns=['State']) df['Year'] =0 df['Age'] = (np.random.random(p)*35+15).astype(int) tochange=df.loc[rng.choice(p, size=num_I+num_R, replace=False),'State'].index df.loc[tochange[0:num_I],'State'] ='I' df.loc[tochange[num_I:num_I+num_R],'State'] ='R'return df
Building the model
# Modelling the decrease of beta over time#np.random.random(size=(p,days))#np.log(4)j=12over =10#10/np.cumsum(np.ones(100))b1 =0.25# original beta = beta value before epidemicb2 =0.05# end beta = beta at the end of epidemicb2 + (b1/np.exp((j+(b1*2.9)-over+1)/(b1*27)))
0.19397059336956546
8#collapse_hidedef seir_model_stoch(beta, beta2, p, num_I, num_R, years, T_Infectious, ART, control):#################################### Explanation of inputs ######################################## As seen in SSA, beta has a starting value, but after a certain threshold (as soon as incidence or prevalence reaches a certain threshold) behaviours change and beta decreases# beta is initial value of beta at start of epidemic (usually 0.3, but can range from 0.2 to 0.5 as seen in SSA)# beta2 is final value (usually around 0.05)# p is total number of individuals in population# num_I is initial number of PLWHA in population (for simulations start with something between 1 and 10 depending on size of p)# num_R is initial number of people deceased from HIV/AIDS# years is number of years you want to run simulation for# T_infectious is distribution of progression of HIV in an individual (use 'gamma' for HIV)# ART is to emulate ART usage:# ART == 0 means no ART# ART == 1 means ART stops evolution of I to R but does not stop spread from I to S# ART == 2 means ART stops both I to R, and S to I# control sets the threshold at which beta above will decrease# control == 0 means no control i.e. beta never decreases# control == 1 means beta decreases once incidence is 15 per 1 thousand population# control == 2 means beta decreases once incidence is 30 per 1 thousand population##################################### Set up the dataframe ###################################### Initialize population dataframe with data given by user df = make_df(p, num_I, num_R)# This variable is used to track daily value of beta if it varies over time xxbeta=np.array([],dtype=float)# Initialize the arrays to return# Below are numbers of S, I, R total S=np.array([],dtype=int) I=np.array([],dtype=int) R=np.array([],dtype=int)# Below are the daily additions in S, I, R Spd=np.array([],dtype=int) Ipd=np.array([],dtype=int) Rpd=np.array([],dtype=int)# Beta values to track spread b=beta b2=np.array([],dtype=float) b1=b# Signal to initiate decrease of beta over =0# signal to end transmission and deaths due to ART art1 =0 art2 =0# Stochastic model so use random values to decide on progression rand = np.random.random(size=(p,years))# Depending if you want exponential, gamma, or Weibull distribution for T_Infectious# Uses distributions found on blog part 3if T_Infectious =='expon': ItoR = expon.rvs(loc=0,scale=10,size=p)elif T_Infectious =='gamma': ItoR = gamma.rvs(4,loc=3,scale=2,size=p) else: ItoR = weibull_min.rvs(2.3, loc=2, scale=20.11, size=p)####################################### Simulation code ######################################### Iterate over every day the simulation is runfor j inrange(0,years-1):# 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 if ART <2: StoI_index = df.loc[(df.State =='S') & (df.Age <49) & (rand[:,j] < b[j]*len(np.where(df.State=='I')[0])/(len(np.where(df.State=='I')[0])+len(np.where(df.State=='S')[0])))].index StoS_index = df.loc[(df.State =='S') & (df.Age <49) & (rand[:,j] >= b[j]*len(np.where(df.State=='I')[0])/(len(np.where(df.State=='I')[0])+len(np.where(df.State=='S')[0])))].indexelif ART ==2: if art2 ==0: StoI_index = df.loc[(df.State =='S') & (df.Age <49) & (rand[:,j] < b[j]*len(np.where(df.State=='I')[0])/(len(np.where(df.State=='I')[0])+len(np.where(df.State=='S')[0])))].index StoS_index = df.loc[(df.State =='S') & (df.Age <49) & (rand[:,j] >= b[j]*len(np.where(df.State=='I')[0])/(len(np.where(df.State=='I')[0])+len(np.where(df.State=='S')[0])))].indexelif art2 ==1: StoI_index = df.loc[(df.State =='S') & (df.Age >55)].index # cannot happen so put an impossible condition like df.Age > 55 to emulate StoS_index = df.loc[(df.State =='S') & (df.Age <49)].index # anyone S under 49 will stay S StoRem_index = df.loc[(df.State =='S') & (df.Age >=49)].index# For each row, if a person has been a certain number of years in I, they will go to R (progression to AIDS and death)# This follows ItoR variable which is either exponential or gamma distributed according to above ItoRem_index = df.loc[(df.State =='I') & (df.Age >=49)].indexif ART ==0: #don't use ART ItoR_index = df.loc[(df.State =='I') & (j-df.Year >= ItoR) & (df.Age <49)].index ItoI_index = df.loc[(df.State =='I') & (j-df.Year < ItoR) & (df.Age <49)].indexelif ART >0:if art2 ==0: ItoR_index = df.loc[(df.State =='I') & (j-df.Year >= ItoR) & (df.Age <49)].index ItoI_index = df.loc[(df.State =='I') & (j-df.Year < ItoR) & (df.Age <49)].indexelif art2 ==1: ItoR_index = df.loc[(df.State =='I') & (df.Age >55)].index # cannot happen so impossible condition ItoI_index = df.loc[(df.State =='I') & (df.Age <49)].index RtoRem_index = df.loc[(df.State =='R') & (df.Age >=49)].index RtoR_index = df.loc[(df.State =='R') & (df.Age <49)].index# Use indexes collected above to populate per day values Ipd = np.append(Ipd,len(StoI_index)) Rpd = np.append(Rpd,len(ItoR_index))# Now we use the indexes collected above randomly to change the actual population dataframe to the new states df.loc[ItoR_index, ['State','Year']] = ['S', j] df.loc[ItoR_index, 'Age'] = df.loc[ItoR_index, 'Age'] +1 df.loc[ItoI_index, 'Age'] = df.loc[ItoI_index, 'Age'] +1 df.loc[StoI_index, ['State','Year']] = ['I', j] df.loc[StoI_index, 'Age'] = df.loc[StoI_index, 'Age'] +1 df.loc[StoS_index, 'Age'] = df.loc[StoS_index, 'Age'] +1 df.loc[RtoR_index, 'Age'] = df.loc[RtoR_index, 'Age'] +1 df.iloc[ItoRem_index] = ['S', j, 15] df.iloc[StoRem_index] = ['S', j, 15] df.iloc[RtoRem_index] = ['S', j, 15]# Append the S, I, and R arrays S=np.append(S,len(np.where(df.State=='S')[0])) I=np.append(I,len(np.where(df.State=='I')[0])) R=np.append(R,len(np.where(df.State=='R')[0]))# Code below for control measures to reduce beta valuesif control ==1:if (I[-1]/p >0.006): art1 =1if over ==0: over = jif art1 ==1:if j > over +15: art2 =1if over !=0: b = beta2 + (b1/np.exp((j+(b1*2.9)-over+1)/(b1*27)))if control ==2:if (I[-1]/p >0.01): art1 =1if over ==0: over = jif art1 ==1:if j > over +15: art2 =1if over !=0: b = beta2 + (b1/np.exp((j+(b1*2.9)-over+1)/(b1*27))) xxbeta2 = ((S[j-1]+I[j-1])/I[j-1])*Ipd[j]/S[j-1] b2 = np.append(b2, xxbeta2) Ipd[0]+=num_I Rpd[0]+=num_Rreturn S, I, R, Spd, Ipd, Rpd, xxbeta, b2, over
Testing the model
#collapse_hide# Define parameters for stochastic modelyears =50p =100000num_E =0num_I =50num_R =0beta_stoch = [0.17,0.17,0.26,0.26,0.36,0.36]#beta_stoch = np.linspace(0.2,0.5,num=10)#beta_stoch = [0.1,0.1,0.1,0.1,0.1,0.1]beta_stoch2 = [0.05,0.05,0.05,0.05,0.05,0.05]#beta_stoch = beta_stoch2control= [1,1,1,1,2,2]n =len(beta_stoch)#results_stoch = []# Run n stochastic simulationsfor i in tqdm.tqdm(range(n)): res = seir_model_stoch(beta_stoch[i]*np.ones(years),beta_stoch2[i], p, num_I, num_R, years, 'gamma', 0, control[i]) results_stoch.append(res)