 ## 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.3
import pandas as pd
import numpy as np
import math
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 tqdm
import time

Collecting plotly==4.14.3
|████████████████████████████████| 13.2MB 254kB/s
Requirement already satisfied: retrying>=1.3.3 in /usr/local/lib/python3.7/dist-packages (from plotly==4.14.3) (1.3.3)
Requirement already satisfied: six in /usr/local/lib/python3.7/dist-packages (from plotly==4.14.3) (1.15.0)
Installing collected packages: plotly
Found existing installation: plotly 4.4.1
Uninstalling plotly-4.4.1:
Successfully uninstalled plotly-4.4.1
Successfully installed plotly-4.14.3


### Creating the initial population dataframe

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, 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=12
over = 10

#10/np.cumsum(np.ones(100))

b1 = 0.25  # original beta = beta value before epidemic
b2 = 0.05 # end beta = beta at the end of epidemic
b2 + (b1/np.exp((j+(b1*2.9)-over+1)/(b1*27)))

0.19397059336956546
8#collapse_hide
def 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 3
if 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 run
for j in range(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'))/(len(np.where(df.State=='I'))+len(np.where(df.State=='S'))))].index
StoS_index = df.loc[(df.State == 'S') & (df.Age < 49) & (rand[:,j] >= b[j]*len(np.where(df.State=='I'))/(len(np.where(df.State=='I'))+len(np.where(df.State=='S'))))].index
elif 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'))/(len(np.where(df.State=='I'))+len(np.where(df.State=='S'))))].index
StoS_index = df.loc[(df.State == 'S') & (df.Age < 49) & (rand[:,j] >= b[j]*len(np.where(df.State=='I'))/(len(np.where(df.State=='I'))+len(np.where(df.State=='S'))))].index
elif 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)].index
if 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)].index
elif 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)].index
elif 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')))
I=np.append(I,len(np.where(df.State=='I')))
R=np.append(R,len(np.where(df.State=='R')))

# Code below for control measures to reduce beta values
if control == 1:
if (I[-1]/p > 0.006):
art1 = 1
if over == 0:
over = j

if art1 == 1:
if j > over + 15:
art2 = 1

if 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 = 1
if over == 0:
over = j

if art1 == 1:
if j > over + 15:
art2 = 1

if 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+=num_I
Rpd+=num_R

return S, I, R, Spd, Ipd, Rpd, xxbeta, b2, over


## Testing the model

#collapse_hide
# Define parameters for stochastic model
years = 50
p = 100000
num_E = 0
num_I = 50
num_R = 0
beta_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_stoch2

control= [1,1,1,1,2,2]

n = len(beta_stoch)

#results_stoch = []

# Run n stochastic simulations
for 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)

100%|██████████| 6/6 [00:41<00:00,  6.95s/it]


#collapse_hide
fig = go.Figure()

for i in range(len(results_stoch)):
#fig.add_trace(go.Scatter(name='Beta_stoch1', x=np.arange(len(results_stoch1)), y=results_stoch1, line={'dash':'dot','color':'yellow'}, legendgroup="Sim_"+str(i)))
fig.add_trace(go.Scatter(name='Beta_meas'+str(i), x=np.arange(len(results_stoch[i])), y=results_stoch[i], line={'dash':'dot','color':'yellow'}, legendgroup="Sim_"+str(i)))
fig.add_trace(go.Scatter(name='I_stoch'+str(i), x=np.arange(len(results_stoch[i])), y=results_stoch[i]/p, line={'dash':'dot', 'color':'red'}, legendgroup="Sim_"+str(i)))
fig.add_trace(go.Bar(name='Ip_stoch'+str(i), x=np.arange(len(results_stoch[i])), y=results_stoch[i]*10/p, legendgroup="Sim_"+str(i)))
fig.add_trace(go.Scatter(name='R_stoch'+str(i), x=np.arange(len(results_stoch[i])), y=results_stoch[i]/p, line={'dash':'dot', 'color':'green'}, legendgroup="Sim_"+str(i)))

fig.update_layout(
xaxis_title = 'Years',
yaxis_title = 'Proportion of population',
title={
'text':r'$\text{Stochastic HIV SIR model}$',
'x':0.5,
'xanchor':'center'
}
)

fig.show()