EPIONCHO-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

1.1 Software

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.

1.2 Run time

When run on laptop machine using Windows, a single run of EPIONCHO-IBM with the activated morbidity module takes at least 30 minutes to locally run on most machines.

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 (18fe7ea0) has not changed since last install.
##   Use `force = TRUE` to force installation

In the updated model (version 1.1.0), morbidity is included, capturing onchocerciasis related skin and ocular disease.

EPIONCHO-IBM has previously been extended to include a OAE module, to simulate OAE prevalence and incidence. Below, we present a flow chart detailing the process overview for the the morbidity module which now features skin and ocular disease, including how this module is parameterised, and how this module connects to the main transmission model.

Skin disease Flowchart Ocular disease Flowchart

We also provide details on how onchocerciasis associated epilepsy (OAE) has also been previously incorporated (see Stapley et al. 2024 Nat Commun), with a detailed practical guide/demo for running the onchocerciasis-associated epilepsy (OAE) module, please see the Running OAE in EPIONCHO-IBM vignette

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). Turning on the epilepsy module (epilepsy module = YES) outputs each individual’s OAE status, assigned, when aged between 3 and 15 years, according to an onset probability derive from [Chesnais et al. 2018]. To run the model to equilibrium without treating in the same run, use the following code:

#length of simulation in years
timesteps = 100

#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
ABR.in = 615 # 50%
output_equilibrium <- ep.equi.sim(time.its = timesteps,
                                  ABR = ABR.in,
                                  N.in = 400,
                                  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,
                                  vector.control.strt = NA,
                                  gam.dis.in = 0.3,
                                  run_equilibrium = TRUE,
                                  morbidity_module = "YES",
                                  equilibrium = NA,
                                  print_progress = TRUE)
## [1] "Ov16 full population store times"
## [1] 36599
## [1] "no MDA to be simulated"
## [1] "1 yrs; 1 %"
## [1] "2 yrs; 2 %"
## [1] "3 yrs; 3 %"
## [1] "4 yrs; 4 %"
## [1] "5 yrs; 5 %"
## [1] "6 yrs; 6 %"
## [1] "7 yrs; 7 %"
## [1] "8 yrs; 8 %"
## [1] "9 yrs; 9 %"
## [1] "10 yrs; 10 %"
## [1] "11 yrs; 11 %"
## [1] "12 yrs; 12 %"
## [1] "13 yrs; 13 %"
## [1] "14 yrs; 14 %"
## [1] "15 yrs; 15 %"
## [1] "16 yrs; 16 %"
## [1] "17 yrs; 17 %"
## [1] "18 yrs; 18 %"
## [1] "19 yrs; 19 %"
## [1] "20 yrs; 20 %"
## [1] "21 yrs; 21 %"
## [1] "22 yrs; 22 %"
## [1] "23 yrs; 23 %"
## [1] "24 yrs; 24 %"
## [1] "25 yrs; 25 %"
## [1] "26 yrs; 26 %"
## [1] "27 yrs; 27 %"
## [1] "28 yrs; 28 %"
## [1] "29 yrs; 29 %"
## [1] "30 yrs; 30 %"
## [1] "31 yrs; 31 %"
## [1] "32 yrs; 32 %"
## [1] "33 yrs; 33 %"
## [1] "34 yrs; 34 %"
## [1] "35 yrs; 35 %"
## [1] "36 yrs; 36 %"
## [1] "37 yrs; 37 %"
## [1] "38 yrs; 38 %"
## [1] "39 yrs; 39 %"
## [1] "40 yrs; 40 %"
## [1] "41 yrs; 41 %"
## [1] "42 yrs; 42 %"
## [1] "43 yrs; 43 %"
## [1] "44 yrs; 44 %"
## [1] "45 yrs; 45 %"
## [1] "46 yrs; 46 %"
## [1] "47 yrs; 47 %"
## [1] "48 yrs; 48 %"
## [1] "49 yrs; 49 %"
## [1] "50 yrs; 50 %"
## [1] "51 yrs; 51 %"
## [1] "52 yrs; 52 %"
## [1] "53 yrs; 53 %"
## [1] "54 yrs; 54 %"
## [1] "55 yrs; 55 %"
## [1] "56 yrs; 56 %"
## [1] "57 yrs; 57 %"
## [1] "58 yrs; 58 %"
## [1] "59 yrs; 59 %"
## [1] "60 yrs; 60 %"
## [1] "61 yrs; 61 %"
## [1] "62 yrs; 62 %"
## [1] "63 yrs; 63 %"
## [1] "64 yrs; 64 %"
## [1] "65 yrs; 65 %"
## [1] "66 yrs; 66 %"
## [1] "67 yrs; 67 %"
## [1] "68 yrs; 68 %"
## [1] "69 yrs; 69 %"
## [1] "70 yrs; 70 %"
## [1] "71 yrs; 71 %"
## [1] "72 yrs; 72 %"
## [1] "73 yrs; 73 %"
## [1] "74 yrs; 74 %"
## [1] "75 yrs; 75 %"
## [1] "76 yrs; 76 %"
## [1] "77 yrs; 77 %"
## [1] "78 yrs; 78 %"
## [1] "79 yrs; 79 %"
## [1] "80 yrs; 80 %"
## [1] "81 yrs; 81 %"
## [1] "82 yrs; 82 %"
## [1] "83 yrs; 83 %"
## [1] "84 yrs; 84 %"
## [1] "85 yrs; 85 %"
## [1] "86 yrs; 86 %"
## [1] "87 yrs; 87 %"
## [1] "88 yrs; 88 %"
## [1] "89 yrs; 89 %"
## [1] "90 yrs; 90 %"
## [1] "91 yrs; 91 %"
## [1] "92 yrs; 92 %"
## [1] "93 yrs; 93 %"
## [1] "94 yrs; 94 %"
## [1] "95 yrs; 95 %"
## [1] "96 yrs; 96 %"
## [1] "97 yrs; 97 %"
## [1] "98 yrs; 98 %"
## [1] "99 yrs; 99 %"

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). However, if you use a individual exposure heterogeneity of 0.2, 0.3, or 0.4, the density dependent parameters will be set automatically, to their respectively fitted density-dependent parameters in Hamley et al. 2019.

3.2 Visualising microfilarial infection and OAE prevalence through time

The output from run_EPIONCHO_IBM_with_OAE is a list with 20 elements.

names(output_equilibrium)
##  [1] "mf_prev"                                 
##  [2] "mf_intens"                               
##  [3] "ov16_seroprevalence_no_seroreversion"    
##  [4] "ov16_seroprevalence_finite_seroreversion"
##  [5] "L3"                                      
##  [6] "ABR"                                     
##  [7] "all_infection_burdens"                   
##  [8] "Ke"                                      
##  [9] "blackfly_l3_prevalence"                  
## [10] "years"                                   
## [11] "all_mf_prevalence_age_grouped"           
## [12] "all_mf_intensity_age_grouped"            
## [13] "ov16_indiv_matrix"                       
## [14] "ov16_timetrend_outputs"                  
## [15] "ov16_timetrend_outputs_adj"              
## [16] "worm_burden_outputs"                     
## [17] "ABR_recorded"                            
## [18] "coverage.recorded"                       
## [19] "percent_never_treated"                   
## [20] "oae_incidence_outputs"                   
## [21] "all_morbidity_prevalence_outputs"        
## [22] "all_equilibrium_outputs"

The first five elements (mf_prev, mf_intens, ov16_seroprevalence_no_seroreversion, ov16_seroprevalence_finite_seroreversion and L3) contain the temporal dynamics (at 1 day increments) for the microfilarial prevalence (individuals age > 5), the population mean microfilarial intensity (individuals age > 5), the temporal dynamics for seroprevalence (w/o and w/ seroreverstion), and the mean number of L3 larvae in the black fly population. To plot these over time, use:

years <- output_equilibrium$years

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

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

plot(years, output_equilibrium$ov16_seroprevalence_no_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ no seroreversion', ylim = c(0, 1))

plot(years, output_equilibrium$ov16_seroprevalence_finite_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ finite seroreversion', ylim = c(0, 1))

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

The ABR attribute lets you know what the input ABR used was. The all_infection_burdens attribute is a matrix that contains the detailed infection burden of each individual at the final timestep. The ov16_indiv_matrix attribute is a matrix/dataframe that contains the detailed ov16 serostatus for each individual at various timesteps. See the Ov16 vignette for more detail. The years attribute lets you know what the timestep is in terms of years from timepoint 0. The ABR_recorded attribute is a temporal output showing the ABR at every year/timestep. The coverage.recorded attribute is a temporal output showing the MDA coverage at every year/timestep. It will have NA values if no treatment was given. The percent_never_treated attribute is a temporal output, showing how many individuals at a given timestep had never been treated. Without MDA, this will be all NA values. The last element,"all_equilibrium_outputs" is the element required to simulate MDA without having to run the model to the endemic equilibrium. Note: The Running_EPIONCHO_IBM.Rmd vignette will have the most up-to-date output parameters.

The all_mf_prevalence_age_grouped, all_mf_intensity_age_grouped, ov16_timetrend_outputs, ov16_timetrend_outputs_adj, oae_incidence_outputs, and all_morbidity_prevalence_outputs elements contain the temporal dynamics (at 1 day increments) for multiple age groups for mf prevalence, intensity, ov16 seroprevalence (w/ and without seroreversion), ov16 seroprevalence adjusted, OAE incidence, and morbidity prevalence (including the OAE prevalence) respectively. The prevalence outputs are grouped by [5, 80], [0-5), [0-5), [5-10), [10-15), [15-20), [20-30), [30-50), [50-80]. Other than 5-80, the rest of the age groups can be customized by editing the output_age_groups input parameter. The OAE incidence outputs are grouped by [5, 80], [0-5), [5-10], [11-15], Males, and Females. These cannot be customized.

The morbidity outputs are shown below, using ggplot2, dplyr, and tidyr:

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.3
process_morbidity_outputs <- function(data) {
  morbidity_outputs_tmp <- data$all_morbidity_prevalence_outputs
  morbidity_outputs_tmp_2 <- morbidity_outputs_tmp[, c(which(endsWith(colnames(morbidity_outputs_tmp), "prev") | colnames(morbidity_outputs_tmp) == "ABR" | colnames(morbidity_outputs_tmp) == "year"))] %>% as.data.frame() 
  years <- data$years
  morbidity_outputs_tmp_2[, "year"] <- years
  morbidity_outputs_return <- morbidity_outputs_tmp_2 %>% pivot_longer(cols=c(which(grepl("prev", colnames(morbidity_outputs_tmp_2)))), names_to="measure", values_to="mean_prevalence")
  return(morbidity_outputs_return)
}
morbidity_outputs <- process_morbidity_outputs(output_equilibrium)



ggplot(morbidity_outputs) + 
  geom_line(
    aes(x=year, y=mean_prevalence * 100)
  ) +
  theme_bw() +
  facet_wrap(~measure) +
  xlab("Years") +
  scale_y_continuous("Prevalence (%)")

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 = 25
treat.strt = 1; treat.stp = 25 #if treatment stops in year 25, the last round is at the beginning of year 24
gv.trt = 1
trt.int = 1

output_treat_annual <- ep.equi.sim(time.its = timesteps,
                                   ABR = ABR.in,
                                   N.in = 400,
                                   treat.timing = NA,
                                   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,
                                   vector.control.strt = NA,
                                   gam.dis.in = 0.3,
                                   run_equilibrium = FALSE,
                                   equilibrium = output_equilibrium$all_equilibrium_outputs,
                                   morbidity_module = "YES",
                                   print_progress = FALSE)
## [1] "Ov16 full population store times"
## [1]  366 9150 9333 9516
## [1] "24 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"
##  [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
## [1] "Target coverage at each round"
## [1] "65%"
## [1] "Treatment to be given at each round"
## [1] "IVM only"
## [1] "ABR is"
## [1] "615 bites person-1 yr-1 at endemic equilibria"

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 and morbidity dynamics through time during annual MDA, use:

years <- output_treat_annual$years

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

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

plot(years, output_treat_annual$ov16_seroprevalence_no_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ no seroreversion', ylim = c(0, 1))

plot(years, output_treat_annual$ov16_seroprevalence_finite_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ finite seroreversion', ylim = c(0, 1))

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

morbidity_outputs_annual <- process_morbidity_outputs(output_treat_annual)
ggplot(morbidity_outputs_annual) + 
  geom_line(
    aes(x=year, y=mean_prevalence * 100)
  ) +
  theme_bw() +
  facet_wrap(~measure) +
  xlab("Years") +
  scale_y_continuous("Prevalence (%)")

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 = 25
treat.strt = 1; treat.stp = 25 #if treatment stops in year 25, the last round is at the beginning of year 24
gv.trt = 1
trt.int = 0.5


output_treat_biannual <- ep.equi.sim(time.its = timesteps,
                                   ABR = ABR.in,
                                   N.in = 400,
                                   treat.timing = NA,
                                   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,
                                   vector.control.strt = NA,
                                   gam.dis.in = 0.3,
                                   run_equilibrium = FALSE,
                                   morbidity_module = "YES",
                                   equilibrium = output_equilibrium$all_equilibrium_outputs,
                                   print_progress = TRUE)
## [1] "Ov16 full population store times"
## [1]  366 9150 9333 9516
## [1] "48 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"
##  [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
## [1] "Target coverage at each round"
## [1] "65%"
## [1] "Treatment to be given at each round"
## [1] "IVM only"
## [1] "ABR is"
## [1] "615 bites person-1 yr-1 at endemic equilibria"
## [1] "1 yrs; 4 %"
## [1] "2 yrs; 8 %"
## [1] "3 yrs; 12 %"
## [1] "4 yrs; 16 %"
## [1] "5 yrs; 20 %"
## [1] "6 yrs; 24 %"
## [1] "7 yrs; 28 %"
## [1] "8 yrs; 32 %"
## [1] "9 yrs; 36 %"
## [1] "10 yrs; 40 %"
## [1] "11 yrs; 44 %"
## [1] "12 yrs; 48 %"
## [1] "13 yrs; 52 %"
## [1] "14 yrs; 56 %"
## [1] "15 yrs; 60 %"
## [1] "16 yrs; 64 %"
## [1] "17 yrs; 68 %"
## [1] "18 yrs; 72 %"
## [1] "19 yrs; 76 %"
## [1] "20 yrs; 80 %"
## [1] "21 yrs; 84 %"
## [1] "22 yrs; 88 %"
## [1] "23 yrs; 92 %"
## [1] "24 yrs; 96 %"

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

years <- output_treat_biannual$years

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

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

plot(years, output_treat_biannual$ov16_seroprevalence_no_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ no seroreversion', ylim = c(0, 1))

plot(years, output_treat_biannual$ov16_seroprevalence_finite_seroreversion, type = 'l', xlab = 'time (years)', ylab = 'seroprevalence w/ finite seroreversion', ylim = c(0, 1))

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

morbidity_outputs_annual <- process_morbidity_outputs(output_treat_biannual)
ggplot(morbidity_outputs_annual) + 
  geom_line(
    aes(x=year, y=mean_prevalence * 100)
  ) +
  theme_bw() +
  facet_wrap(~measure) +
  xlab("Years") +
  scale_y_continuous("Prevalence (%)")

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), or using a correlation parameter for a more dynamic compliance structure (comp.correlation = 0.3). More information on correlation can be seen in the complex MDA vignette.

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/