MGDrivE 3: A decoupled vector-human framework for epidemiological simulation of mosquito genetic control tools and their surveillance

Novel mosquito genetic control tools, such as CRISPR-based gene drives, hold great promise in reducing the global burden of vector-borne diseases. As these technologies advance through the research and development pipeline, there is a growing need for modeling frameworks incorporating increasing levels of entomological and epidemiological detail in order to address questions regarding logistics and biosafety. Epidemiological predictions are becoming increasingly relevant to the development of target product profiles and the design of field trials and interventions, while entomological surveillance is becoming increasingly important to regulation and biosafety. We present MGDrivE 3 (Mosquito Gene Drive Explorer 3), a new version of a previously-developed framework, MGDrivE 2, that investigates the spatial population dynamics of mosquito genetic control systems and their epidemiological implications. The new framework incorporates three major developments: i) a decoupled sampling algorithm allowing the vector portion of the MGDrivE framework to be paired with a more detailed epidemiological framework, ii) a version of the Imperial College London malaria transmission model, which incorporates age structure, various forms of immunity, and human and vector interventions, and iii) a surveillance module that tracks mosquitoes captured by traps throughout the simulation. Example MGDrivE 3 simulations are presented demonstrating the application of the framework to a CRISPR-based homing gene drive linked to dual disease-refractory genes and their potential to interrupt local malaria transmission. Simulations are also presented demonstrating surveillance of such a system by a network of mosquito traps. MGDrivE 3 is freely available as an open-source R package on CRAN (https://cran.r-project.org/package=MGDrivE2) (version 2.1.0), and extensive examples and vignettes are provided. We intend the software to aid in understanding of human health impacts and biosafety of mosquito genetic control tools, and continue to iterate per feedback from the genetic control community.


Equilibrium distribution based on baseline entomological innoculation rate (EIR)
and age structure of the population; and 3. Population-level immunity functions which modulate various infection probabilities.
The state space is modeled as a set of partial differential equations (PDEs).The infection states are: susceptible (S), treated symptomatic disease (T), untreated symptomatic disease (D), asymptomatic infection that is detectable by rapid diagnostic test (RDT) (A), sub-patent infection that is undetectable by RDT (U), and post-treatment prophylaxis (P).The force of infection on humans (which depends on the EIR) is denoted Λ, the probability that symptoms develop after an infectious challenge is denoted Φ, and the fraction of clinical cases that receive effective treatment is denoted f T .The set of human state PDEs is shown below, with a representing age and t representing time.
April 23, 2024 1/7 Here, d i indicates the mean duration of state i.Additionally, the model includes four forms of population-level immunity: • Pre-erythrocytic immunity, I B , reduces the probability of infection if bitten by an infectious mosquito; • Acquired and maternal clinical immunity, I CA and I CM , represent the effects of blood stage immunity in reducing the probability of developing clinical symptoms; and • Detection immunity, I D , represents the effect of blood stage immunity in reducing the detectability of an infection and onward transmission to mosquitoes.
The PDEs describing immunity are below.Note that ε represents the EIR, u i limits the rate at which immunity can be boosted at high exposure for immunity state i, and d i determines the duration of immunity for immunity state i.
gives the rate equation for the susceptible (S) state for age category i where η i gives the aging rate from S i −→ S i+1 and similarly η i−1 gives the aging rate from S i−1 −→ S i .
For the youngest age group, the η i−1 S i−1 term would be left out, and for the oldest age group, the η i S i term would be left out.
One implementation note is that the model assumes a fixed latent period of 12 days after an infectious challenge from a mosquito, after which either symptoms develop or an asymptomatic infection proceeds.Because of this fixed delay, the equations are technically formulated as "delay differential equations," where the current state depends on the previous state.
To initialize the distribution of disease and immunity states, the model takes as input the baseline EIR, the age structure of the population, the proportion of treated cases, and baseline entomological parameters.Some of the mosquito life cycle parameters will vary in the presence of interventions, which will be described in the next section.
Finally, an important novel contribution of this work incorporates a model for genotype-specific transmission probabilities.In the gene drive context, it is important to understand how mosquitoes modified with a certain allele can affect disease transmission.In the traditional context [5], force of infection on humans (λ H ) is proportional to the EIR (ε) and the probability of successful infection upon biting (b).
In the ICL malaria model, the force of infection is expanded to include a term corresponding to pre-erythrocytic immunity (I B ): In our adapted model, we allow for varying transmission probabilities depending on the genotype distribution of circulating mosquitoes in the model.For example, we may consider a gene drive system in which the homozygous transgenic mosquito (denoted "HH") confers perfect infection blocking such that b HH = 0, and where wildtype mosquitoes (denoted "WW") have an infection probability of 0.55 (b W W = 0.55).Then, to calculate the total infection probability for any time point, t, we take the weighted average over the circulating proportion of infectious mosquitoes of each genotype (p HH and p W W ) at time t, i.e.: More generally, for a genotype set G, we have the total infection probability as the weighted average: where p g (t) represents the population frequency of mosquitoes having genotype g at time t, and b g represents the infection probability for genotype g.With all of the above components in place, the epidemiology model is fully specified.

Epidemiological outcomes
In this modeling framework, it is important to specify the epidemiological outcomes of interest.Generally, we are interested in clinical incidence of disease, which refers to the April 23, 2024 3/7 number of new symptomatic cases per day, and P. falciparum pathogen prevalence (PfPr ), which refers to the proportion of the population harboring the malaria pathogen, whether symptomatic or asymptomatic.Because our model is age-structured, we can consider these outcomes for different age categories.Malaria outcomes are often more severe for younger individuals [6], so it makes sense to consider incidence and prevalence for younger age groups (e.g., 0-2 years) in order to better understand how gene drive interventions will mitigate these cases.Additionally, we can consider all-ages prevalence and incidence to understand how the intervention will perform in the entire population.
Mathematically, we define PfPr as the sum of all individuals in infectious disease states: symptomatic and treated (T), symptomatic and untreated (D), asymptomatic patent infection (A), and asymptomatic subpatent infection (U).Therefore, the all-ages (often denoted by the subscript 0-99 to denote the entire lifespan in years) pathogen prevalence at a given time point, t, is given by: where A is the set of all age compartments.Similarly, the 0-2 years PfPr is given by: As for clinical incidence, we first define some parameters: • ϕ: the probability of acquiring clinical disease upon infection (proportional to immunity levels via a Hill function); • λ H : the force of infection on humans (linearly proportional to the EIR, ε); and • Y : the sum of non-clinical disease states, susceptible (S), asymptomatic patent infection (A), and subpatent infection (U).
Then we can define the all-ages clinical incidence as: and the 0-2 years clinical incidence as: Generally, we are interested in these outcomes with respect to their baseline or pre-intervention values.In our analyses, we will calculate the reduction in prevalence and clinical incidence as our outcomes of interest.As we will be running many stochastic repetitions of the simulation for a given parameter set, the mean reduction over the repetition set and simulation timespan will be used.Note that in this formulation, each disease state is a proportion, with all disease states summing to 1.If instead we wish to model a population of N H humans, then we would simply divide each outcome by N H to obtain the proportional value.

Additional interventions
Here we show the full derivation of how indoor residual spraying (IRS), long-lasting insecticide-treated nets (LLIN), and artemisisin-based combination therapy (ACT) April 23, 2024 4/7 interventions modify baseline mosquito life cycle parameters.This derivation is adapted from previous work [4,7].
First, we assume that, at baseline, we have three proportions of active vector control interventions, {χ IRS , χ LLIN , χ ACT }, which represent the proportion of humans in the model covered by the given intervention.Then, χ ACT corresponds to the proportion of symptomatically infected humans that are treated upon infection, f T .
Then, {χ IRS , χ LLIN } jointly modify various mosquito life cycle parameters.First, we model the impact of LLINs and IRS on the length of the mosquito gonotrophic cycle (i.e., the time taken for a mosquito to take a blood meal and lay eggs before seeking its next blood meal).This time can be divided into τ 1 (the time spent foraging) and τ 2 (the time spent ovipositing and resting).Then, the length of the gonotrophic cycle in the presence of vector control is given by: where τ 1 (0, 0) represents the time time spent foraging with LLIN and IRS coverages of zero, and: Here, Q 0 represents the human blood index, θ B represents the proportion of bites taken on a person in bed, θ I represents the proportion of bites taken on a person outdoors, r IRS represents the probability of repeating a feeding attempt in the presence of IRS, r IRS,LLIN represents the probability of repeating a feeding attempt in the presence of IRS and LLINs, and: Then, with the modified gonotrophic cycle calculated (δ C ), we can model the impact of LLINs and IRS on the adult mosquito death rate.We express the mortality rate in the presence of vector control as: where p represents the probability of an adult mosquito surviving one day.Then we can break down p into two components p 1 (the probability of surviving the mosquito stage) and p 2 (the probability of surviving the blood meal stage): where: z is the same as above and w gives the probability that a mosquito successfully feeds and survives a single feeding attempt: Here, s LLIN and s IRS represent the probability of feeding and surviving in the presence of LLINs and IRS, respectively.The non-intervention survival probabilities are given by: Now, we have mathematical expressions for the gonotrophic cycle length and adult mortality rate (δ c and µ V,c respectively).We can finally model the impact of LLINs and IRS on the egg laying rate of the adult female mosquito.In the absence of vector control, the egg laying rate is given by: where ε is the number of viable eggs laid per oviposition cycle.Then, with the previously defined parameters, the egg laying rate in the presence of vector control interventions is simply: Finally, we can modify the human biting rate per mosquito in the presence of LLINs and IRS.In the absence of interventions, the biting rate is given by: The biting index under intervention is given by: where w is the calculated probability from above.Then, using the modified gonotrophic cycle length previously derived (δ c ), the modified biting rate is thus: With these definitions in place, we have fully specified the impact of vector control interventions on mosquito life cycle parameters.