EPIONCH0-IBM is a stochastic, individual-based model for onchocerciasis transmission. A mathematical description can be found at https://doi.org/10.1371/journal.pntd.0007557.

1. System requirements

EPIONCHO IBM is written in R (1) which is freely available for installation on Windows, Mac OS and Linux platforms from https://www.r-project.org. It is also recommended to install R Studio https://www.rstudio.com which provides a helpful user interface for R. EPIONCHO-IBM has been tested on Mac OS and Windows, but should be compatible with Linux and run on any desktop or laptop machine.

2. Installation guide

The model package can be downloaded and installed directly from a GitHub repository (https://github.com/mrc-ide/EPIONCHO.IBM). The remotes (2) package must be installed to do this. Installation should take no more than a few seconds on most desktop or laptop machines.

## Skipping install of 'EPIONCHO.IBM' from a github remote, the SHA1 (a2b3f1ea) has not changed since last install.
##   Use `force = TRUE` to force installation

3. Demo

3.1 Simulating an endemic equilibrium without MDA

Before simulating treatment, the model must be simulated to the endemic equilibrium. It is best to run the model to equilibrium (run_equilibrium = TRUE, give.treat = 0), save the output and use this as the starting point for MDA. Note that running the model to the equilibrium can be slow, taking in excess of 5 minutes (depending on machine performance). To run the model to equilibrium without treating in the same run, use the following code:

#length of simulation in years
timesteps = 30

#should treatment be given, when and how often
give.treat.in = 0
treat.strt = 1; treat.stp = 16
trt.int = 1

#annual biting rate, which determines infection prevalence (60% microfilarae prevalence)
ABR.in = 1082

output_equilibrium <- ep.equi.sim(time.its = timesteps,
                                  ABR = ABR.in,
                                  N.in = 440,
                                  treat.int = trt.int,
                                  treat.prob = 0.65,
                                  give.treat = give.treat.in,
                                  treat.start = treat.strt,
                                  treat.stop = treat.stp,
                                  pnc = 0.05,
                                  min.mont.age = 5,
                                  delta.hz.in = 0.186,
                                  delta.hinf.in = 0.003,
                                  c.h.in = 0.005,
                                  gam.dis.in = 0.3,
                                  run_equilibrium = TRUE,
                                  equilibrium,
                                  print_progress = TRUE)
## [1] "no MDA to be simulated"
## [1] "1 yrs; 3.3 %"
## [1] "2 yrs; 6.7 %"
## [1] "3 yrs; 10 %"
## [1] "4 yrs; 13.3 %"
## [1] "5 yrs; 16.7 %"
## [1] "6 yrs; 20 %"
## [1] "7 yrs; 23.3 %"
## [1] "8 yrs; 26.7 %"
## [1] "9 yrs; 30 %"
## [1] "10 yrs; 33.3 %"
## [1] "11 yrs; 36.7 %"
## [1] "12 yrs; 40 %"
## [1] "13 yrs; 43.3 %"
## [1] "14 yrs; 46.7 %"
## [1] "15 yrs; 50 %"
## [1] "16 yrs; 53.3 %"
## [1] "17 yrs; 56.7 %"
## [1] "18 yrs; 60 %"
## [1] "19 yrs; 63.3 %"
## [1] "20 yrs; 66.7 %"
## [1] "21 yrs; 70 %"
## [1] "22 yrs; 73.3 %"
## [1] "23 yrs; 76.7 %"
## [1] "24 yrs; 80 %"
## [1] "25 yrs; 83.3 %"
## [1] "26 yrs; 86.7 %"
## [1] "27 yrs; 90 %"
## [1] "28 yrs; 93.3 %"
## [1] "29 yrs; 96.7 %"

To stop the printing of the time as the model runs, set print_progress = FALSE. Currently, with print_progress = TRUE, when the model run reaches each year in the simulation, the year and % progress will be printed. If run_equilibrium = TRUE then no input is expected for equilibrium, which is where a previously simulated endemic equilibrium can be entered. This will be discussed later in section 3.3. Above, no MDA is not requested by setting gv.trt = 0. Note that values must be entered for treatment start, stop and interval (1 is annual, 0.5 is biannual), e.g treat.strt = 0; treat.stp = 16; trt.int = 1, even if gv.trt = 0.

It is also possible to manually change the density-dependent parameters relating to the establishment of the adult Onchocerca volvulus (delta.hz.in, delta.hinf.in and c.h.in) and individual exposure heterogeneity in humans (gam.dis.in), in the function calling the model (above). This feature enables the user to modify these parameters to capture the combinations of exposure heterogeneity parameter values and fitted density-dependent parameters in Hamley et al. 2019.

3.2 Visualising microfilarial infection prevalence through time

The output from run_EPIONCHO_IBM is a list with four elements.

names(output_equilibrium)
## [1] "mf_prev"                 "mf_intens"              
## [3] "L3"                      "all_equilibrium_outputs"

The first three elements contain the temporal dynamics (at 1 day increments) for the microfilarial prevalence (individuals age > 5), the population mean microfilarial intensity (individuals age > 5) and the mean number of L3 larvae in the black fly population. To plot these over time, use:

tme <- seq(1, 30*366-1)/366

plot(tme, output_equilibrium$mf_prev, type = 'l', xlab = 'time (years)', ylab = 'microfilarial prevalence', ylim = c(0, 1))

plot(tme, output_equilibrium$mf_intens, type = 'l', xlab = 'time (years)', ylab = 'mean microfilarial intensity')

plot(tme, output_equilibrium$L3, type = 'l', xlab = 'time (years)', ylab = 'mean L3 in fly population')

The fourth element,"all_equilibrium_outputs" is the element required to simulate MDA without having to run the model to the endemic equilibrium.

3.3 Simulating annual and biannual mass drug administration with ivermectin

3.3.1 Annual MDA

Assuming you have already run the model and saved the output as described in section 3.1, you can simulate 15 years of annual mass drug administration with ivermectin by running:

timesteps = 30 
treat.strt = 1; treat.stp = 26 #if model run stops in year 30, the last round is at the beginning of year 25
gv.trt = 1
trt.int = 1

output_treat_annual <- ep.equi.sim(time.its = timesteps,
                                   ABR = ABR.in,
                                   N.in = 440,
                                   treat.int = trt.int,
                                   treat.prob = 0.65,
                                   give.treat = gv.trt,
                                   treat.start = treat.strt,
                                   treat.stop = treat.stp,
                                   pnc = 0.05,
                                   min.mont.age = 5,
                                   delta.hz.in = 0.186,
                                   delta.hinf.in = 0.003,
                                   c.h.in = 0.005,
                                   gam.dis.in = 0.3,
                                   run_equilibrium = FALSE,
                                   equilibrium = output_equilibrium[[4]],
                                   print_progress = FALSE)
## [1] "25 MDA rounds to be given"
## [1] "MDA given at"
##  [1] "1yrs"  "2yrs"  "3yrs"  "4yrs"  "5yrs"  "6yrs"  "7yrs"  "8yrs"  "9yrs" 
## [10] "10yrs" "11yrs" "12yrs" "13yrs" "14yrs" "15yrs" "16yrs" "17yrs" "18yrs"
## [19] "19yrs" "20yrs" "21yrs" "22yrs" "23yrs" "24yrs" "25yrs"
##  [1]  367  733 1099 1465 1831 2197 2563 2929 3295 3661 4027 4393 4759 5125 5491
## [16] 5857 6223 6589 6955 7321 7687 8053 8419 8785 9151

If in the above model run print_progress = FALSE is called, but gv.trt = 1 is also specified, only the duration of MDA will be printed (no progress updates). Currently, with print_progress = TRUE, when the model run reaches each year in the simulation, the year and % progress will be printed as before. Note that timesteps is the total number of years for which the model is run. treat.start and treat.stop are the years within this duration at which MDA is given.

To visualise the infection dynamics through time during annual MDA, use:

tme <- seq(0, 30*366-1)/366

plot(tme, output_treat_annual$mf_prev, type = 'l', xlab = 'time', ylab = 'microfilarial prevalence', ylim = c(0, 1))

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

plot(tme, output_treat_annual$mf_intens, type = 'l', xlab = 'time', ylab = 'mean microfilarial intensity')

abline(v = seq(1, 25), col = 'grey', lwd = 0.1)

3.3.2 Biannual MDA

Assuming you have run the model to the endemic equilibrium (as described in secton 3.1), temporal infection trends during biannual MDA can be obtained by running:

timesteps = 30 
treat.strt = 1; treat.stp = 26 #if treatment stops in year 16, the last round is at the beginning of year 15
gv.trt = 1
trt.int = 0.5


output_treat_biannual <- ep.equi.sim(time.its = timesteps,
                                   ABR = ABR.in,
                                   N.in = 440,
                                   treat.int = trt.int,
                                   treat.prob = 0.65,
                                   give.treat = gv.trt,
                                   treat.start = treat.strt,
                                   treat.stop = treat.stp,
                                   pnc = 0.05,
                                   min.mont.age = 5,
                                   delta.hz.in = 0.186,
                                   delta.hinf.in = 0.003,
                                   c.h.in = 0.005,
                                   gam.dis.in = 0.3,
                                   run_equilibrium = FALSE,
                                   equilibrium = output_equilibrium[[4]],
                                   print_progress = TRUE)
## [1] "50 MDA rounds to be given"
## [1] "MDA given at"
##  [1] "1yrs"    "1.5yrs"  "2yrs"    "2.5yrs"  "3yrs"    "3.5yrs"  "4yrs"   
##  [8] "4.5yrs"  "5yrs"    "5.5yrs"  "6yrs"    "6.5yrs"  "7yrs"    "7.5yrs" 
## [15] "8yrs"    "8.5yrs"  "9yrs"    "9.5yrs"  "10yrs"   "10.5yrs" "11yrs"  
## [22] "11.5yrs" "12yrs"   "12.5yrs" "13yrs"   "13.5yrs" "14yrs"   "14.5yrs"
## [29] "15yrs"   "15.5yrs" "16yrs"   "16.5yrs" "17yrs"   "17.5yrs" "18yrs"  
## [36] "18.5yrs" "19yrs"   "19.5yrs" "20yrs"   "20.5yrs" "21yrs"   "21.5yrs"
## [43] "22yrs"   "22.5yrs" "23yrs"   "23.5yrs" "24yrs"   "24.5yrs" "25yrs"  
## [50] "25.5yrs"
##  [1]  367  550  733  916 1099 1282 1465 1648 1831 2014 2197 2380 2563 2746 2929
## [16] 3112 3295 3478 3661 3844 4027 4210 4393 4576 4759 4942 5125 5308 5491 5674
## [31] 5857 6040 6223 6406 6589 6772 6955 7138 7321 7504 7687 7870 8053 8236 8419
## [46] 8602 8785 8968 9151 9334
## [1] "1 yrs; 3.3 %"
## [1] "2 yrs; 6.7 %"
## [1] "3 yrs; 10 %"
## [1] "4 yrs; 13.3 %"
## [1] "5 yrs; 16.7 %"
## [1] "6 yrs; 20 %"
## [1] "7 yrs; 23.3 %"
## [1] "8 yrs; 26.7 %"
## [1] "9 yrs; 30 %"
## [1] "10 yrs; 33.3 %"
## [1] "11 yrs; 36.7 %"
## [1] "12 yrs; 40 %"
## [1] "13 yrs; 43.3 %"
## [1] "14 yrs; 46.7 %"
## [1] "15 yrs; 50 %"
## [1] "16 yrs; 53.3 %"
## [1] "17 yrs; 56.7 %"
## [1] "18 yrs; 60 %"
## [1] "19 yrs; 63.3 %"
## [1] "20 yrs; 66.7 %"
## [1] "21 yrs; 70 %"
## [1] "22 yrs; 73.3 %"
## [1] "23 yrs; 76.7 %"
## [1] "24 yrs; 80 %"
## [1] "25 yrs; 83.3 %"
## [1] "26 yrs; 86.7 %"
## [1] "27 yrs; 90 %"
## [1] "28 yrs; 93.3 %"
## [1] "29 yrs; 96.7 %"

Note than above trt.int = 0.5 (rather than trt.int = 1) to request biannual MDA. The output can be visualised using:

tme <- seq(0, 30*366-1)/366

plot(tme, output_treat_biannual$mf_prev, type = 'l', xlab = 'time', ylab = 'microfilarial prevalence', ylim = c(0, 1))

abline(v = seq(1, 25, 0.5), col = 'grey', lwd = 0.1)

plot(tme, output_treat_biannual$mf_intens, type = 'l', xlab = 'time', ylab = 'mean microfilarial intensity')

abline(v = seq(1, 25, 0.5), col = 'grey', lwd = 0.1)

3.3.3 Coverage and non-compliance

The population level coverage is a proportion and is entered via treat.prob. Note that because children <5 do not receive treatment, coverage of 1 cannot be achieved. Additionally, there is the option to set the proportion of individuals who never take treatment (pnc = 0.05).

4. References

  1. R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2018).

  2. Csárdi, G. RCurl: Package ‘remotes’: R Package Installation from Remote Repositories, Including ‘GitHub’. (2022). https://cran.r-project.org/web/packages/remotes/