#collapse_hide
# Need this new function for model below:
def make_df(p, num_I, num_R):
= pd.DataFrame(np.full((p,1), 'S').T[0], columns=['State'])
df 'Year'] = 0
df['Age'] = (np.random.random(p)*35+15).astype(int)
df[=df.loc[rng.choice(p, size=num_I+num_R, replace=False),'State'].index
tochange0:num_I],'State'] = 'I'
df.loc[tochange[+num_R],'State'] = 'R'
df.loc[tochange[num_I:num_Ireturn df
Building the model
#np.random.random(size=(p,days))
#np.log(4)
=11
j= 10
over #10/np.cumsum(np.ones(100))
0.05 + (0.3/np.exp((j+1-over)/10))
#collapse_hide
def seir_model_stoch(beta, beta2, p, num_I, num_R, years, T_Infectious, ART, control):
# Initialize population dataframe with data given by user
= make_df(p, num_I, num_R)
df
# This variable is used to track daily value of beta if it varies over time
=np.array([],dtype=float)
xxbeta
# Initialize the arrays to return
# Below are numbers of S, I, R total
=np.array([],dtype=int)
S=np.array([],dtype=int)
I=np.array([],dtype=int)
R# Below are the daily additions in S, I, R
=np.array([],dtype=int)
Spd=np.array([],dtype=int)
Ipd=np.array([],dtype=int)
Rpd
=beta
b#b2=beta[0]
=np.array([],dtype=float)
b2=b
b1
# signal diminshing beta
= 0
over
# signal end of deaths due to ART
= 0
art1 = 0
art2
# Stochastic model so use random values to decide on progression
= np.random.random(size=(p,years))
rand
# Depending if you want exponential, gamma, or Weibull distribution for T_Infectious
# Uses distributions found on blog part 3
if T_Infectious == 'expon':
= expon.rvs(loc=0,scale=10,size=p)
ItoR elif T_Infectious == 'gamma':
= gamma.rvs(4,loc=3,scale=2,size=p)
ItoR else:
= weibull_min.rvs(2.3, loc=2, scale=20.11, size=p)
ItoR
# Iterate over every day the simulation is run
for j in range(0,years-1):
# Record daily beta values
=np.append(xxbeta, b[j])
xxbeta
# 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
if ART < 2:
= 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
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 elif ART == 2:
if art2 == 0:
= 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
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 elif art2 == 1:
= df.loc[(df.State == 'S') & (df.Age > 55)].index
StoI_index = df.loc[(df.State == 'S') & (df.Age < 49)].index
StoS_index
= df.loc[(df.State == 'S') & (df.Age == 49)].index
StoRem_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
= df.loc[(df.State == 'I') & (df.Age == 49)].index
ItoRem_index if ART == 0: #don't use ART
= df.loc[(df.State == 'I') & (j-df.Year >= ItoR) & (df.Age < 49)].index
ItoR_index = df.loc[(df.State == 'I') & (j-df.Year < ItoR) & (df.Age < 49)].index
ItoI_index elif ART > 0:
if art2 == 0:
= df.loc[(df.State == 'I') & (j-df.Year >= ItoR) & (df.Age < 49)].index
ItoR_index = df.loc[(df.State == 'I') & (j-df.Year < ItoR) & (df.Age < 49)].index
ItoI_index elif art2 ==1:
= df.loc[(df.State == 'I') & (df.Age > 49)].index
ItoR_index = df.loc[(df.State == 'I') & (df.Age < 49)].index
ItoI_index
= df.loc[(df.State == 'R') & (df.Age == 49)].index
RtoRem_index = df.loc[(df.State == 'R') & (df.Age < 49)].index
RtoR_index
# Use indexes collected above to populate per day values
#Epd = np.append(Epd,len(StoE_index))
#Ipd = np.append(Ipd,len(EtoI_index))
= np.append(Ipd,len(StoI_index))
Ipd = np.append(Rpd,len(ItoR_index))
Rpd
# Now we use the indexes collected above randomly to change the actual population dataframe to the new states
= ['S', j, 15]
df.iloc[ItoRem_index] 'State','Year']] = ['S', j]
df.loc[ItoR_index, ['Age'] = df.loc[ItoR_index, 'Age'] + 1
df.loc[ItoR_index, 'Age'] = df.loc[ItoI_index, 'Age'] + 1
df.loc[ItoI_index, = ['S', j, 15]
df.iloc[StoRem_index] 'State','Year']] = ['I', j]
df.loc[StoI_index, ['Age'] = df.loc[StoI_index, 'Age'] + 1
df.loc[StoI_index, 'Age'] = df.loc[StoS_index, 'Age'] + 1
df.loc[StoS_index,
= ['S', j, 15]
df.iloc[RtoRem_index] 'Age'] = df.loc[RtoR_index, 'Age'] + 1
df.loc[RtoR_index,
# Append the S, I, and R arrays
=np.append(S,len(np.where(df.State=='S')[0]))
S=np.append(I,len(np.where(df.State=='I')[0]))
I=np.append(R,len(np.where(df.State=='R')[0]))
R
# Code below for control measures to reduce beta values
if control == 1:
if (I[-1]/p > 0.015):
= 1
art1 if over == 0:
= j
over
if art1 == 1:
if j > over + 15:
#if Ipd[-2] > Ipd[-1]:
= 1
art2
if over != 0:
#b = beta2+(b1/np.exp((j+3-over)/15))
= beta2+(b1/np.exp((j+1-over)/10))
b
if control == 2:
if (I[-1]/p > 0.3):
= 1
art1 if over == 0:
= j
over #print(over)
if art1 == 1:
if j > over + 15:
#if Ipd[-2] > Ipd[-1]:
= 1
art2
if over != 0:
#b = beta2+(b1/np.exp((j+3-over)/15))
= beta2+(b1/np.exp((j+1-over)/10))
b
= ((S[j-1]+I[j-1])/I[j-1])*Ipd[j]/S[j-1]
xxbeta2 #xxbeta2 = 0.5
#print(xxbeta2)
= np.append(b2, xxbeta2)
b2
#Epd[0]+=num_E
0]+=num_I
Ipd[0]+=num_R
Rpd[
#return S,E,I,R, Epd, Ipd, Rpd, xxbeta
return S, I, R, Spd, Ipd, Rpd, xxbeta, b2, over
Testing the model
#collapse_hide
# Define parameters for stochastic model
= 200
days = 10000
p = 0
num_E = 1
num_I = 0
num_R = 0.3*np.ones(days)
beta_stoch = 0.05
beta_stoch2
# Run 3 stochastic simulations
= seir_model_stoch(beta_stoch,beta_stoch2, p, num_I, num_R, years, 'gamma', 0, 1)
results_stoch1 = seir_model_stoch(beta_stoch, beta_stoch2, p, num_I, num_R, years, 'gamma', 0, 1)
results_stoch2 = seir_model_stoch(beta_stoch, beta_stoch2, p, num_I, num_R, years, 'gamma', 0, 2)
results_stoch3 = seir_model_stoch(beta_stoch, beta_stoch2, p, num_I, num_R, years, 'gamma', 0, 2) results_stoch4
8] results_stoch1[
#collapse_hide
= go.Figure(data=[
fig ='Beta_stoch1', x=np.arange(len(results_stoch1[0])), y=results_stoch1[6], line={'dash':'dot','color':'yellow'}, legendgroup="stoch1"),
go.Scatter(name='Beta_meas1', x=np.arange(len(results_stoch1[0])), y=results_stoch1[7], line={'dash':'dot','color':'yellow'}, legendgroup="stoch1"),
go.Scatter(name='I_stoch1', x=np.arange(len(results_stoch1[0])), y=results_stoch1[1]/p, line={'dash':'dot', 'color':'red'}, legendgroup="stoch1"),
go.Scatter(name='Ip_stoch1', x=np.arange(len(results_stoch1[0])), y=results_stoch1[4]*10/p, legendgroup="stoch1"),
go.Bar(name='R_stoch1', x=np.arange(len(results_stoch1[0])), y=results_stoch1[2]/p, line={'dash':'dot', 'color':'green'}, legendgroup="stoch1"),
go.Scatter(name='Beta_stoch2', x=np.arange(len(results_stoch2[0])), y=results_stoch2[6], line={'dash':'dot','color':'yellow'}, legendgroup="stoch2"),
go.Scatter(name='Beta_meas2', x=np.arange(len(results_stoch2[0])), y=results_stoch2[7], line={'dash':'dot','color':'yellow'}, legendgroup="stoch2"),
go.Scatter(name='I_stoch2', x=np.arange(len(results_stoch2[0])), y=results_stoch2[1]/p, line={'dash':'dot', 'color':'red'}, legendgroup="stoch2"),
go.Scatter(name='Ip_stoch2', x=np.arange(len(results_stoch2[0])), y=results_stoch2[4]*10/p, legendgroup="stoch2"),
go.Bar(name='R_stoch2', x=np.arange(len(results_stoch2[0])), y=results_stoch2[2]/p, line={'dash':'dot', 'color':'green'}, legendgroup="stoch2"),
go.Scatter(name='Beta_stoch3', x=np.arange(len(results_stoch3[0])), y=results_stoch3[6], line={'dash':'dot', 'color':'yellow'}, legendgroup="stoch3"),
go.Scatter(name='Beta_meas3', x=np.arange(len(results_stoch3[0])), y=results_stoch3[7], line={'dash':'dot','color':'yellow'}, legendgroup="stoch3"),
go.Scatter(name='I_stoch3', x=np.arange(len(results_stoch3[0])), y=results_stoch3[1]/p, line={'dash':'dot', 'color':'red'}, legendgroup="stoch3"),
go.Scatter(name='Ip_stoch3', x=np.arange(len(results_stoch3[0])), y=results_stoch3[4]*10/p, legendgroup="stoch3"),
go.Bar(name='R_stoch3', x=np.arange(len(results_stoch3[0])), y=results_stoch3[2]/p, line={'dash':'dot', 'color':'green'}, legendgroup="stoch3"),
go.Scatter(name='Beta_stoch4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[6], line={'dash':'dot', 'color':'yellow'}, legendgroup="stoch4"),
go.Scatter(name='Beta_meas4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[7], line={'dash':'dot','color':'yellow'}, legendgroup="stoch4"),
go.Scatter(name='I_stoch4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[1]/p, line={'dash':'dot', 'color':'red'}, legendgroup="stoch4"),
go.Scatter(name='Ip_stoch4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[4]*10/p, legendgroup="stoch4"),
go.Bar(name='R_stoch4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[2]/p, line={'dash':'dot', 'color':'green'}, legendgroup="stoch4")
go.Scatter(name
])
fig.update_layout(= 'Day',
xaxis_title = 'Proportion of population',
yaxis_title ={
title'text':r'$\text{Effect of stochasticity on Deterministic SEIR model}$',
'x':0.5,
'xanchor':'center'
}
)
fig.show()