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 7, 2020

#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

#np.random.random(size=(p,days))
#np.log(4)
j=11
over = 10
#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
    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)

    b=beta
    #b2=beta[0]
    b2=np.array([],dtype=float)
    b1=b

    # signal diminshing beta
    over = 0 

    # signal end of 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)

    # 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 
        #StoE_index = df.loc[(df.State == 'S') & (rand[:,j] < b[j]*len(np.where(df.State=='I')[0])/p)].index
        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
            StoS_index = df.loc[(df.State == 'S') & (df.Age < 49)].index

        StoRem_index = df.loc[(df.State == 'S') & (df.Age == 49)].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
        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 > 49)].index
            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
        #Epd = np.append(Epd,len(StoE_index))
        #Ipd = np.append(Ipd,len(EtoI_index))
        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.iloc[ItoRem_index] = ['S', j, 15]
        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.iloc[StoRem_index] = ['S', j, 15]
        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.iloc[RtoRem_index] = ['S', j, 15]
        df.loc[RtoR_index, 'Age'] = df.loc[RtoR_index, 'Age'] + 1
        
        
        # 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.015):
            art1 = 1
            if over == 0:
              over = j
          
          if art1 == 1:
            if j > over + 15:    
            #if Ipd[-2] > Ipd[-1]:
              art2 = 1

          if over != 0:
            #b = beta2+(b1/np.exp((j+3-over)/15))
            b = beta2+(b1/np.exp((j+1-over)/10))

        if control == 2:
          if (I[-1]/p > 0.3):
            art1 = 1
            if over == 0:
              over = j
              #print(over)
          
          if art1 == 1:
            if j > over + 15:    
            #if Ipd[-2] > Ipd[-1]:
              art2 = 1

          if over != 0:
            #b = beta2+(b1/np.exp((j+3-over)/15))
            b = beta2+(b1/np.exp((j+1-over)/10))


        xxbeta2 = ((S[j-1]+I[j-1])/I[j-1])*Ipd[j]/S[j-1]
        #xxbeta2 = 0.5
        #print(xxbeta2)
        b2 = np.append(b2, xxbeta2)
                
    #Epd[0]+=num_E
    Ipd[0]+=num_I
    Rpd[0]+=num_R

    #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
days = 200
p = 10000
num_E = 0
num_I = 1
num_R = 0
beta_stoch = 0.3*np.ones(days)
beta_stoch2 = 0.05

# Run 3 stochastic simulations
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, 1)
results_stoch3 = seir_model_stoch(beta_stoch, beta_stoch2, p, num_I, num_R, years, 'gamma', 0, 2)
results_stoch4 = seir_model_stoch(beta_stoch, beta_stoch2, p, num_I, num_R, years, 'gamma', 0, 2)
results_stoch1[8]
#collapse_hide
fig = go.Figure(data=[    
    go.Scatter(name='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.Bar(name='Ip_stoch1', x=np.arange(len(results_stoch1[0])), y=results_stoch1[4]*10/p, legendgroup="stoch1"),
    go.Scatter(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.Bar(name='Ip_stoch2', x=np.arange(len(results_stoch2[0])), y=results_stoch2[4]*10/p, legendgroup="stoch2"),
    go.Scatter(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.Bar(name='Ip_stoch3', x=np.arange(len(results_stoch3[0])), y=results_stoch3[4]*10/p, legendgroup="stoch3"),
    go.Scatter(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.Bar(name='Ip_stoch4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[4]*10/p, legendgroup="stoch4"),
    go.Scatter(name='R_stoch4', x=np.arange(len(results_stoch4[0])), y=results_stoch4[2]/p, line={'dash':'dot', 'color':'green'}, legendgroup="stoch4")
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
    title={
        'text':r'$\text{Effect of stochasticity on Deterministic SEIR model}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()