HIV modeling

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\))
  1. 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.

Code 
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 plotly
import plotly.io as pio
from IPython.display import display, HTML

## Tomas Mazak's workaround
plotly.offline.init_notebook_mode()
display(HTML(
    #'<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_SVG"></script>'
    '<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js?config=TeX-MML-AM_SVG"></script>'
    '<script src="https://cdn.plot.ly/plotly-3.0.1.js" charset="utf-8"></script>'
))

pio.renderers.default = "plotly_mimetype+notebook_connected"

pio.templates.default = "plotly_dark"
import tqdm
import time

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

Code 
# 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

Code 
# 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
Code 
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

Code 
# 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:16<00:00,  2.75s/it]
Code 
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':'Stochastic HIV SIR model',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()