HIV modeling

• Building a stochastic model
  • Creating the initial population dataframe
  • Building the model
• Testing the model

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
  Downloading https://files.pythonhosted.org/packages/1f/f6/bd3c17c8003b6641df1228e80e1acac97ed8402635e46c2571f8e1ef63af/plotly-4.14.3-py2.py3-none-any.whl (13.2MB)
     |████████████████████████████████| 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[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=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')[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])))].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')[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])))].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')[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 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[0]+=num_I
    Rpd[0]+=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[0])), y=results_stoch1[6], 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][0])), y=results_stoch[i][7], 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][0])), y=results_stoch[i][1]/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][0])), y=results_stoch[i][4]*10/p, legendgroup="Sim_"+str(i)))
  fig.add_trace(go.Scatter(name='R_stoch'+str(i), x=np.arange(len(results_stoch[i][0])), y=results_stoch[i][2]/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()