BARS Transmission Model: Chicago

– Aditya, Nick, Jonathan, and Nikki

Today, we will see results from the Chicago model, for a 25,000 day simulation (approx. 68 years).

The model includes:

Currently, we are modeling a testing-and-diagnosis piece. The implementation for the starting network is complete. Temporal simulation of testing and diagnosis is in progress. More details are here.

Following this step, we will include PrEP.

Other points for discussion: Lag time between diagnosis and ART initiation; ART/PrEP cessation; mixing parameters; matching simulated and empirical age structure.

Below, is an “interrogation” of the model, to demonstrate its various features.

These datasets include the sexual network at the last time step, biomarker data with detailed trajectories, infection, death and partnership events, and counts of various quantities at each time step.

Results

Demography

We first compute the annual growth rate.

The final population size is 5301, corresponding to a growth rate of 0.21%. (The rate of entry of new individuals is a free parameter, and will be adjusted to reflect data form the three cities).

The age distribution at the end of the simulation is below.

age <- net%v%"age"
plot(
  qplot(age, geom="histogram", binwidth=5, 
       col=I("red"), 
       alpha=I(.2))
  )

plot of chunk unnamed-chunk-4

Prevalence

The overall prevalence at the end of the simulation was about 13.1%. HIV prevalence as a function of time is shown below.

plot of chunk unnamed-chunk-6

Incidence

The annual incidence rates over the course of the simulation are shown below.

plot of chunk unnamed-chunk-7

The mean annual incidence, averaged over the length of the simulation, is 1.33%.

Sexual Networks

Main Partnerships

Mean degree.

The target main mean degree is 0.46 and a plot of the simulated mean degree, over the course of the simulation, is below.

plot of chunk unnamed-chunk-8

Degree distribution.

The simulated and target momentary degree distributions are below. plot of chunk unnamed-chunk-9

Partnership duration.

## Loading required package: data.table
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, last
## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt

plot of chunk unnamed-chunk-10

The target mean duration is 911. In the simulated data above, the interquartile range is 149, 714.

Casual Partnerships

Mean degree.

The target casual mean degree is 0.59 and a plot of the simulated mean degree, over the course of the simulation, is below.

plot of chunk unnamed-chunk-11

Degree distribution.

The simulated and target momentary degree distributions are below.

## degree0 degree1 degree2 degree3 degree4 degree5 
##    2951    1839     401      90      17       3

plot of chunk unnamed-chunk-12

Partnership duration.

pt_data_casual <- pt_data[network_type == 1]

# group by p1, p2
pt_data_casual[,.(counts = .N),by=.(p1,p2),][,unique(counts)] # 1 2 4 3
## [1] 1 2 3 4
pt_data_casual[,.(counts = .N),by=.(p1,p2),][,.N,by=counts]
##    counts    N
## 1:      1 3257
## 2:      2 9504
## 3:      3    4
## 4:      4    1
# so focus on just 2 is justified
    ## (old DONOTUSE: res <- pt_data_casual[,.(tick, counts = .N),by=.(p1,p2),][counts == 2,abs(diff(tick))])

res <- pt_data_casual[,.(tick, counts = .N),by=.(p1,p2),][counts == 2,.(dur = abs(diff(tick))), by=.(p1,p2)][,dur]

summary(res)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0    81.0   192.0   274.5   383.0  2240.0
res_dt <- as.data.frame(res)
  qplot(res, geom="histogram", binwidth=100, 
       col=I("red"), 
       alpha=I(.2))

plot of chunk unnamed-chunk-13

sum.res <- summary(res)

The target mean duration is 335. In the simulated data above, the interquartile range is 81, 383.

Overlap between main and casual partnerships.

plot of chunk unnamed-chunk-14

Total mean degree: main and casual partnerships combined.

plot of chunk unnamed-chunk-15

Testing and Diagnosis

p_inf_status <- person_data$infection.status
p_time_of_infection <- person_data$time.of.infection
p_art_status <- person_data$art.status
p_time_of_art_initiation <- person_data$time.of.art.initiation

table(p_inf_status)
## p_inf_status
##    0    1 
## 5078  915
table(p_art_status)
## p_art_status
##    0    1 
## 5162  831
xtabs(~factor(p_inf_status, exclude=NULL) + 
        factor(p_art_status, exclude=NULL))
##                                     factor(p_art_status, exclude = NULL)
## factor(p_inf_status, exclude = NULL)    0    1
##                                    0 5078    0
##                                    1   84  831
summary(p_time_of_art_initiation - p_time_of_infection)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2403.0     0.0     0.0   159.6     0.0  3533.0
neg_diff_art.init_inf <- neg.diff <- which(p_time_of_art_initiation - p_time_of_infection < 0)

p_inf_status[neg_diff_art.init_inf]
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1
p_time_of_infection[neg_diff_art.init_inf]
##  [1]    5  337 1022  393 1123  532 1007 1726  701 2280 1544 1501 2354 1853
## [15] 1265 1717 1950 2267 1508 1346  310 2305  500 1629 2288 1124 2204 2263
## [29]  497 1188 2127 1582  801 1969  536 1241 2225 1087  605 1208  459 2320
## [43] 2130 1540 1372 1899 2289 1728 1911 2350 2000 2355  416 2253 1006 1667
## [57]  388 1847 1717 1695 1181 2313  585 1040 1719 1001 1274  935 2303 2313
## [71] 2353 1619 2402
p_time_of_art_initiation[neg_diff_art.init_inf]
##  [1] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
## [24] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
## [47] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
## [70] -1 -1 -1 -1
summary(p_time_of_art_initiation)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -1.00   -1.00   -1.00   77.27   -1.00 2483.00

ART metrics

The ART initiation portion was simplified so that all infected's initate ART 1 year after seroconversion. At presetnt, the proportion of the infected on ART was

  on.art <- which(net%v%"art.status" == 1)
  length(on.art)

[1] 685

  length(on.art)/length(infected)

[1] 0.9870317

0.99. The total proportion of all people in the population who were on ART is 12.92%.

Biomarkers

   ## Selected upto n.to.follow infecteds at random
   n_to_follow <- 10
   infectees <- inf_event_data$infectee
   uniq_biom_pid <- unique(biom_data$p_id)

   infectee_uniq_biom_pid <- uniq_biom_pid[which(uniq_biom_pid %in% infectees)]
   new <- biom_data[which(biom_data$p_id == infectee_uniq_biom_pid[1]),]

   if (length(infectee_uniq_biom_pid > n_to_follow)){
     infectee_uniq_biom_pid <- infectee_uniq_biom_pid[1:n_to_follow]#first 10
     #infectee_uniq_biom_pid <- sample(infectee_uniq_biom_pid, 
                                      #n_to_follow, replace=FALSE)# random 10
   }

   for (i in 2:length(infectee_uniq_biom_pid)){
     new_entry <- biom_data[which(biom_data$p_id == infectee_uniq_biom_pid[i]),]
     new <- rbind(new, new_entry)
   }

   new$p_id <- as.character(new$p_id)

   par(mfrow=c(2,1))
   ggplot(new, aes(x=tick/365, y=viral_load, color=p_id))+
          geom_line()+xlab("year")

plot of chunk unnamed-chunk-18

   ggplot(new, aes(x=tick/365, y=cd4_count, color=p_id))+
          geom_line()+xlab("year")

plot of chunk unnamed-chunk-18

Example CD4 and viral load trajectories for 10 randomly selected, HIV-infected individuals are shown above. Individuals on ART are identifiable as those whose CD4 counts start to recover after decreasing initially; individuals who do not go on ART do experience a monotonic decrease.

Conclusion

The features we have in the model appear to work well so far. We are in the process of adding more (as mentioned above), and continuing to monitor emergent trends in the behavior of key parameters and outcomes.

Suggestions welcome!