Package 'MGDrivE'

Aggregate Female Output by Genotype

Description

Aggregate over male mate genotype to convert female matrix output into vector output.

Usage

aggregateFemales(
  readDir,
  writeDir = NULL,
  genotypes,
  remFile = TRUE,
  verbose = TRUE
)

Arguments

readDir	Directory to read input from
writeDir	Directory to write output to. Default is readDir
genotypes	Character vector of possible genotypes; found in driveCube$genotypesID
remFile	Boolean flag to remove original (unaggregated) file
verbose	Chatty? Default is TRUE

Examples

## Not run: 
# This example assumes user has already run MGDrivE and generated output.
#  This also assumes that the user has already split output by patch.
# See vignette for complete example.

# set read/write directory
fPath <- "path/to/data/containing/folder"

# Need genotypes from the cube run in the simulation
#  This is dependent on the simulation run
#  Using Mendelian cube for this example
cube <- cubeMendelian()

# no return value from function
aggregateFemales(readDir= fPath, writeDir = NULL, genotypes = cube$genotypesID,
                 remFile = TRUE)

## End(Not run)

---

Aggregate Output Over Landscape

Description

This function aggregates the output of a run over the entire output, i.e., all of the patches. It writes the output one level above the folder pointed to by readDir, if writeDir is NULL. Output consists of 2 csv files, one for males and one for females, "...M_LandscapeAgg_Run...csv".

Usage

aggregateOutput(readDir, writeDir=NULL)

Arguments

readDir	Directory where output was written to
writeDir	Directory to write output to. Default is one level above readDir

Examples

## Not run: 
# This assumes user has run MGDrivE and output is in fPath.
#  See vignette for examples on how to run MGDrivE

# read/write dirs
fPath <- "folder/containing/output"
oPath <- "folder/to/write/stuff"

# first, split output by patch and aggregate females by mate genotype
# remember, cube is for example and changes with simulation
#  landscape aggregation will work if females are not aggregated, but it's slower
cube <- cubeMendelian()

splitOutput(readDir = fPath, writeDir = NULL, remFile = TRUE)
aggregateFemales(readDir= fPath, writeDi = NULL, genotypes = cube$genotypesID,
                 remFile = TRUE)

# aggregate mosquitoes over entire landscape
#  no return value
aggregateOutput(readDir = fPath, writeDir = NULL)

## End(Not run)

---

Make List of Batch Migration Parameters

Description

Sets up a list containing the probability of a batch migration, the fractional amount of males/females that migrate, and the weighted probabilities for where to migrate. The default weights for migration are equal for all patches. These can be changed after running the function. This is only used in oneDay_Migration_Stochastic_Network.

Usage

basicBatchMigration(
  batchProbs = 1e-05,
  sexProbs = c(0.01, 0.01),
  numPatches = 1
)

Arguments

batchProbs	Probability of a batch migration, either 1 number or a vector of length equal to the number of patches
sexProbs	Population fraction of males and females that migrate. Either a vector c(M,F) or matrix of 2 columns
numPatches	Number of patches in the simulation

Examples

# to setup for 3 patches
batchMigration = basicBatchMigration(batchProbs = 1e-5, sexProbs = c(0.1, 0.01), numPatches = 3)

---

Make List of Modified Mosquito Releases

Description

Sets up a release schedule for a single patch, returns a list to be used in oneDay_releases_Patch or oneDay_eggReleases_Patch. This function is no longer intended to be used alone, please use the standard interface, generateReleaseVector.

Usage

basicRepeatedReleases(releaseStart, releaseEnd, releaseInterval, releaseMatrix)

Arguments

releaseStart	Day releases start
releaseEnd	Day releases end
releaseInterval	Interval between releases
releaseMatrix	Numeric matrix specifying the genotype and release amount

Examples

## Not run: 
# Setup for 3 patches but only release in the first with a defined release
#  schedule, for the cube cubeHomingDrive:

patchReleases = replicate(n = 3, expr = {
  list(maleReleases = NULL, femaleReleases = NULL, eggReleases = NULL, matedFemaleReleases = NULL)
},simplify = FALSE)

patchReleases[[1]]$femaleReleases = MGDrivE::basicRepeatedReleases(releaseStart = 5,
                                                          releaseEnd = 30,
                                                          releaseInterval = 5,
                                                          releaseMatrix = matrix(c(5,100),1,2))

patchReleases[[1]]$maleReleases = MGDrivE::basicRepeatedReleases(releaseStart = 50,
                                                        releaseEnd = 60,
                                                        releaseInterval = 1,
                                                        releaseMatrix = matrix(c(5,100),1,2))

## End(Not run)

---

Calculate Aquatic Stage Survival Probability

Description

Calculate $\theta_{st}$, density-independent survival probability, given by:

$\theta_{st}=(1-\mu_{st})^{T_{st}}$

Usage

calcAquaticStageSurvivalProbability(mortalityRate, stageDuration)

Arguments

mortalityRate	Daily mortality probability, $\mu_{st}$
stageDuration	Duration of aquatic stage, $T^{st}$

---

Calculate Average Generation Time

Description

Calculate $g$, average generation time, given by:

$g=T_e+T_l+T_p+\frac{1}{\mu_{ad}}$

Usage

calcAverageGenerationTime(stagesDuration, adultMortality)

Arguments

stagesDuration	Vector of lengths of aquatic stages, $T_{e}, T_{l}, T_{p}$
adultMortality	Adult mortality rate, $\mu_{ad}$

---

Calculate Geodesic Distance - Cosine Method

Description

This function calculates geodesic distance using the cosine method.

Usage

calcCos(latLongs, r = 6378137)

Arguments

latLongs	Two column matrix of latitudes/longitudes
r	Earth radius. Default is WGS-84 radius

Examples

# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# cosine distance formula
distMat = calcCos(latLongs = latLong)

---

Calculate Density-dependent Larval Mortality

Description

Calculate $\alpha$, the strength of density-dependent mortality during the larval stage, given by:

$\alpha=\Bigg( \frac{1/2 * \beta * \theta_e * Ad_{eq}}{R_m-1} \Bigg) * \Bigg( \frac{1-(\theta_l / R_m)}{1-(\theta_l / R_m)^{1/T_l}} \Bigg)$

Usage

calcDensityDependentDeathRate(
  fertility,
  thetaAq,
  tAq,
  adultPopSizeEquilibrium,
  populationGrowthRate
)

Arguments

fertility	Number of eggs per oviposition for wild-type females, $\beta$
thetaAq	Vector of density-independent survival probabilities of aquatic stages, $\theta_{e}, \theta_{l}$
tAq	Vector of lengths of aquatic stages, $T_{e}, T_{l}, T_{p}$
adultPopSizeEquilibrium	Adult population size at equilibrium, $Ad_{eq}$
populationGrowthRate	Population growth in absence of density-dependent mortality $R_{m}$

---

Calculate Exponential Stochastic Matrix

Description

Given a distance matrix from calcVinEll, calculate a stochastic matrix where one step movement probabilities follow an exponential density.

Usage

calcExpKernel(distMat, rate)

Arguments

distMat	Distance matrix from calcVinEll
rate	Rate parameter of Exponential distribution

Details

The distribution and density functions for the exponential kernel are given below:

$F(x)=1-e^{-\lambda x}$

$f(x)=\lambda e^{-\lambda x}$

where $\lambda$ is the rate parameter of the exponential distribution.

Examples

# setup distance matrix
# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Ellipsoid  distance formula
distMat = calcVinEll(latLongs = latLong)

# calculate exponential distribution over distances
#  rate is just for example
kernMat = calcExpKernel(distMat = distMat, rate = 10)

---

Calculate Gamma Stochastic Matrix

Description

Given a distance matrix from calcVinEll, calculate a stochastic matrix where one step movement probabilities follow a gamma density.

Usage

calcGammaKernel(distMat, shape, rate)

Arguments

distMat	Distance matrix from calcVinEll
shape	Shape parameter of GammaDist distribution
rate	Rate parameter of GammaDist distribution

Details

The distribution and density functions for the gamma kernel are given below:

$F(x)=\frac{1}{\Gamma(\alpha)}\gamma(\alpha,\beta x)$

$f(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}$

where $\Gamma(\alpha)$ is the Gamma function, $\gamma(\alpha,\beta x)$ is the lower incomplete gamma function, and $\alpha,\beta$ are the shape and rate parameters, respectively.

Examples

# setup distance matrix
# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Ellipsoid  distance formula
distMat = calcVinEll(latLongs = latLong)

# calculate gamma distribution over distances
#  shape and rate are just for example
kernMat = calcGammaKernel(distMat = distMat, shape = 1, rate = 1)

---

Calculate Geodesic Distance - Haversine Method

Description

This function calculates geodesic distance using the Haversine method.

Usage

calcHaversine(latLongs, r = 6378137)

Arguments

latLongs	Two column matrix of latitudes/longitudes
r	Earth radius. Default is WGS-84 radius

Examples

# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Haversine distance formula
distMat = calcHaversine(latLongs = latLong)

---

Calculate Zero-inflated Exponential Stochastic Matrix

Description

Given a distance matrix from calcVinEll, calculate a stochastic matrix where one step movement probabilities follow an zero-inflated exponential density with a point mass at zero. The point mass at zero represents the first stage of a two-stage process, where mosquitoes decide to stay at their current node or leave anywhere. This parameter can be calculated from lifetime probabilities to stay at the current node with the helper function calcZeroInflation.

Usage

calcHurdleExpKernel(distMat, rate, p0)

Arguments

distMat	Distance matrix from calcVinEll
rate	Rate parameter of Exponential distribution
p0	Point mass at zero

Details

If a mosquito leaves its current node, with probability $1-p_{0}$, it then chooses a destination node according to a standard exponential density with rate parameter $rate$.

The distribution and density functions for the zero inflated exponential kernel are given below:

$F(x)=p_{0}\theta(x) + (1-p_{0})(1-e^{-\lambda x})$

$f(x)=p_{0}\delta(x)+(1-p_{0})\lambda e^{-\lambda x}$

where $\lambda$ is the rate parameter of the exponential distribution, $\theta(x)$ is the Heaviside step function and $\delta(x)$ is the Dirac delta function.

Examples

# setup distance matrix
# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Ellipsoid  distance formula
distMat = calcVinEll(latLongs = latLong)

# calculate hurdle exponential distribution over distances
#  rate and point mass are just for example
kernMat = calcHurdleExpKernel(distMat = distMat, rate = 1/1e6, p0 = 0.1)

---

Calculate Distribution of Larval Population

Description

This hidden function calculates the distribution of larvae through time by treating the larval-stage as a discrete-time Markov chain, and solving for the stationary distribution. As the only aquatic population known for initializing MGDrivE is the equilibrium larval population, this acts as an anchor from which to calculate the egg and pupae distributions from (see set_initialPopulation_Patch).

Usage

calcLarvalDist(mu, t)

Arguments

mu	Double, death rate
t	Integer, stage time

---

Calculate Equilibrium Larval Population

Description

Equilibrium larval population size to sustain adult population.

Usage

calcLarvalPopEquilibrium(alpha, Rm)

Arguments

alpha	See calcDensityDependentDeathRate
Rm	See calcPopulationGrowthRate

---

Calculate Larval Stage Mortality Rate

Description

Calculate $\mu_{l}$, the larval mortality, given by

$\mu_l=1-\Bigg( \frac{R_m * \mu_{ad}}{1/2 * \beta * (1-\mu_m)} \Bigg)^{\frac{1}{T_e+T_l+T_p}}$

Usage

calcLarvalStageMortalityRate(
  generationPopGrowthRate,
  adultMortality,
  fertility,
  aquaticStagesDuration
)

Arguments

generationPopGrowthRate	See calcPopulationGrowthRate
adultMortality	Adult mortality rate, $\mu_{ad}$
fertility	Number of eggs per oviposition for wild-type females, $\beta$
aquaticStagesDuration	Vector of lengths of aquatic stages, $T_{e}, T_{l}, T_{p}$

---

Calculate Lognormal Stochastic Matrix

Description

Given a distance matrix from calcVinEll, calculate a stochastic matrix where one step movement probabilities follow a lognormal density.

Usage

calcLognormalKernel(distMat, meanlog, sdlog)

Arguments

distMat	Distance matrix from calcVinEll
meanlog	Log mean of Lognormal distribution
sdlog	Log standard deviation of Lognormal distribution

Details

The distribution and density functions for the lognormal kernel are given below:

$F(x)=\frac{1}{2} + \frac{1}{2} \mathrm{erf}[\frac{\mathrm{ln}x-\mu}{\sqrt{2}\sigma}]$

$f(x)=\frac{1}{x\sigma\sqrt{2\pi}}\mathrm{exp}\left( -\frac{(\mathrm{ln}x-\mu)^{2}}{2\sigma^{2}} \right)$

where $\mu$ is the mean on the log scale, and $\sigma$ is the standard deviation on the log scale.

Examples

# setup distance matrix
# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Ellipsoid  distance formula
distMat = calcVinEll(latLongs = latLong)

# calculate lognormal distribution over distances
#  mean and standard deviation are just for example
kernMat = calcLognormalKernel(distMat = distMat, meanlog = 100, sdlog = 10)

---

Solve for Omega (additional genotype-specific mortality)

Description

Solves for root of equation of geometrically-distributed lifespan for value of omega.

Usage

calcOmega(mu, lifespanReduction)

Arguments

mu	Daily mortality probability (discrete-time hazard, called muAd in code)
lifespanReduction	Target reduced lifespan, between 0 and 1 (target average lifespan will be $\frac{1}{\mu_{Ad}} \times lifespanReduction$)

Examples

# reduce lifespan by 10%
#  Example mu is an average for Aedes
newOmega <- calcOmega(mu = 0.11, lifespanReduction = 0.90)

---

Calculate Generational Population Growth Rate

Description

Calculate $R_{m}$, population growth in absence of density-dependent mortality, given by:

$(r_{m})^{g}$

Usage

calcPopulationGrowthRate(dailyPopGrowthRate, averageGenerationTime)

Arguments

dailyPopGrowthRate	Daily population growth rate, $r_{m}$
averageGenerationTime	See calcAverageGenerationTime

---

Summary Statistics for Stochastic MGDrivE

Description

This function reads in all repetitions for each patch and calculates either the mean, quantiles, or both. User chooses the quantiles, up to 4 decimal places, and enters them as a vector. Quantiles are calculated empirically. (order does not matter)

Usage

calcQuantiles(readDir, writeDir, mean = TRUE, quantiles = NULL, verbose = TRUE)

Arguments

readDir	Directory to find repetition folders in
writeDir	Directory to write output
mean	Boolean, calculate mean or not. Default is TRUE
quantiles	Vector of quantiles to calculate. Default is NULL
verbose	Chatty? Default is TRUE

Details

Given the readDir, this function assumes the follow file structure:

• readDir

  • repetition 1

    • patch 1

    • patch 2

    • patch 3

  • repetition 2

    • patch 1

    • patch 2

    • patch 3

  • repetition 3

  • repetition 4

  • ...
    

Output files are *.csv contain the mean or quantile in the file name, i.e. M/FMean(patchNum).csv and M/FQuantile(quantNum)_(patchNum).csv.

Value

Writes output to files in writeDir

Examples

## Not run: 
# This function assumes network$multRun() has been performed, or several
#  network$oneRun() have been performed and all of the data has been split
#  and aggregated.

# read/write paths
fPath <- "path/to/folder/ofFolders/with/data"
oPath <- "my/path/output"

# here, only calculate mean, no quantiles
#  no return value
calcQuantiles(readDir = fPath, writeDir = oPath, mean = TRUE,
              quantiles = NULL)

# here, calculate 2.5% and 97.5% quantiles
calcQuantiles(readDir = fPath, writeDir = oPath, mean = FALSE,
              quantiles = c(0.025, 0.975))

## End(Not run)

---

Calculate Geodesic Distance - Vincenty Ellipsoid Method

Description

This function calculates geodesic distance using the original Vincenty Ellipsoid method.

Usage

calcVinEll(
  latLongs,
  a = 6378137,
  b = 6356752.3142,
  f = 1/298.257223563,
  eps = 1e-12,
  iter = 100
)

Arguments

latLongs	Two column matrix of latitudes/longitudes
a	Equatorial radius of the earth, default is WGS-84 radius
b	Polar radius of the earth, default is WGS-84 radius
f	Flattening or inverse eccentricity, default eccentricity is WGS-84
eps	Convergence criteria
iter	Maximum number of iterations to attempt convergence

Examples

# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Ellipsoid  distance formula
distMat = calcVinEll(latLongs = latLong)

---

Calculate Geodesic Distance - Vincenty Sphere Method

Description

This function calculates geodesic distance using the Vincenty sphere method.

Usage

calcVinSph(latLongs, r = 6378137)

Arguments

latLongs	Two column matrix of latitudes/longitudes
r	Earth radius. Default is WGS-84 radius

Examples

# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Sphere  distance formula
distMat = calcVinSph(latLongs = latLong)

---

Calculates the zero-inflation part of a hurdle exponential kernel.

Description

Given the probability of an adult mosquito to stay in the same patch throughout its whole lifespan, and its mortality, it calculates the height of the pulse-density part of the hurdle kernel.

Usage

calcZeroInflation(stayThroughLifespanProbability, adultMortality)

Arguments

stayThroughLifespanProbability	Probability of a mosquito to spend its whole lifespan in the same node
adultMortality	Adult mortality rate

Examples

# setup distance matrix
# two-column matrix with latitude/longitude, in degrees
latLong = cbind(runif(n = 5, min = 0, max = 90),
                runif(n = 5, min = 0, max = 180))

# Vincenty Ellipsoid  distance formula
distMat = calcVinEll(latLongs = latLong)

# get hurdle height
# Lets assume 80% stay probs and adult mortality of 0.1
hHeight <- calcZeroInflation(stayThroughLifespanProbability = 0.80,
                             adultMortality = 0.1)

# calculate hurdle exponential distribution over distances
kernMat = calcHurdleExpKernel(distMat = distMat, rate = 10, p0 = hHeight)

---

Export a Cube to .csv

Description

Export a cube as multiple .csv files (one for each genotype; slices of z-axis). This function will create the directory if it doesn't exist. Files are stored as slice_(z-slice)_(genotype).csv

Usage

cube2csv(cube, directory, digits = 3)

Arguments

cube	A cube object (see MGDrivE-Cube for options)
directory	Directory to write .csv files to
digits	Number of significant digits to retain in .csv output

Examples

## Not run: 
# output directory
oPath <- "path/to/write/output"

# setup inheritance cube for export, using Mendelian as the example
cube <- cubeMendelian()

# write out
cube2csv(cube = cube, directory = oPath, digits = 3)

## End(Not run)

---

Inheritance Cube: ClvR (Cleave and Rescue)

Description

Based on the Cleave-and-Rescue system of Oberhofer, this is a 2-locus Cas9-based toxin-antidote system. The first locus carries the Cas9, gRNAs, and a recoded copy of an essential gene. The second locus is the targeted essential gene. This gene can be completely haplosufficient (hSuf = 1) or completely haploinsufficient (hSuf = 0). It is assumed that having 2 copies of the gene (be it wild-type at the second locus or recoded at the first) confers complete viability.

Usage

cubeClvR(
  cF = 1,
  crF = 0,
  ccF = cF,
  ccrF = crF,
  cM = 1,
  crM = 0,
  ccM = cM,
  ccrM = crM,
  dW = 0,
  drW = 0,
  ddW = dW,
  ddrW = drW,
  hSuf = 1,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cF	Female cutting rate, one ClvR allele
crF	Female functional resistance rate, one ClvR allele
ccF	Female cutting rate, two ClvR alleles
ccrF	Female functional resistance rate, two ClvR alleles
cM	Male cutting rate, one ClvR allele
crM	Male functional resistance rate, one ClvR allele
ccM	Male cutting rate, two ClvR alleles
ccrM	Male functional resistance rate, two ClvR alleles
dW	Female deposition cutting rate
drW	Female deposition functional resistance rate
ddW	Female deposition (HH) cutting rate
ddrW	Female deposition (HH) functional resistance rate
hSuf	Haplosufficiency level, default is completely sufficient
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: 2-Locus ClvR (Cleave and Rescue)

Description

Based on the Cleave-and-Rescue system of Oberhofer, this is a 3-locus Cas9-based toxin-antidote system. The first locus carries the Cas9, the second locus carries the gRNAs, and a recoded copy of an essential gene. The third locus is the targeted essential gene. This gene can be completely haplosufficient (hSuf = 1) or completely haploinsufficient (hSuf = 0). It is assumed that having 2 copies of the gene (be it wild-type at the second locus or recoded at the first) confers complete viability. Additionally, loci 1 and 2 can be linked, given crM and crF, imitating the original 2-locus ClvR system. For this construct, the first locus will have 2 alleles, the second will have 2 alleles, and the third will have 3 alleles:

• Locus 1

  • W: Wild-type

  • C: Cas9

• Locus 2

  • W: Wild-type

  • G: gRNAs and recoded essential gene

• Locus 3

  • W: Wild-type

  • R: Functional resistant

  • B: Non-functional resistant

Usage

cubeClvR2(
  cF = 1,
  crF = 0,
  ccF = cF,
  ccrF = crF,
  cM = 1,
  crM = 0,
  ccM = cM,
  ccrM = crM,
  dW = 0,
  drW = 0,
  ddW = dW,
  ddrW = drW,
  hSuf = 1,
  crossF = 0.5,
  crossM = 0.5,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cF	Female cutting rate, one ClvR allele
crF	Female functional resistance rate, one ClvR allele
ccF	Female cutting rate, two ClvR alleles
ccrF	Female functional resistance rate, two ClvR alleles
cM	Male cutting rate, one ClvR allele
crM	Male functional resistance rate, one ClvR allele
ccM	Male cutting rate, two ClvR alleles
ccrM	Male functional resistance rate, two ClvR alleles
dW	Female deposition cutting rate
drW	Female deposition functional resistance rate
ddW	Female deposition (HH) cutting rate
ddrW	Female deposition (HH) functional resistance rate
hSuf	Haplosufficiency level, default is completely sufficient
crossF	Crossover rate in females, 0 is completely linked, 0.5 is unlinked, 1.0 is complete divergence
crossM	Crossover rate in males, 0 is completely linked, 0.5 is unlinked, 1.0 is complete divergence
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Details

Female deposition is implemented incorrectly. Right now, it is performed on male alleles prior to zygote formation - it should happen post-zygote formation. Since this construct doesn't have HDR, this should be fine.
Additionally, it is assumed that deposition requries loaded Cas9-RNP complexes from the mother, having Cas9 and no maternal gRNA, even in the presence of paternal gRNA, will not result in maternal deposition mediated cleavage.

Copy-number dependent rates are based on Cas9, not gRNA. The assumption is that RNA is easier to produce, and therefore won't limit cleavage by Cas9.

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: ECHACR

Description

This function creates an ECHACR construct, it has 5 alleles at the first locus and 4 alleles at the second.

• W: Wild-type

• H: Homing allele

• E: Eraser allele

• R: No-cost resistance allele

• B: Detrimental resistance allele

• cHW: Rate of homing from H, W -> H transition

• cEH: Rate of homing from E, H -> E transition

• cEW: Rate of homing from E, W -> E transition

Usage

cubeECHACR(
  cHW = 1,
  cEW = 1,
  cEH = 1,
  chHW = 0,
  crHW = 0,
  ceEW = 0,
  crEW = 0,
  ceEH = 0,
  crEH = 0,
  d1 = 0,
  d2 = 0,
  d3 = 0,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cHW	Cutting efficiency of drive allele at locus 1
cEW	Cutting efficiency of ECHACR element into W
cEH	Cutting efficiency of ECHACR element into H
chHW	Homing efficiency of drive allele at locus 1
crHW	Resistance allele efficiency of drive allele at locus 1
ceEW	Homing efficiency of ECHACR element into W
crEW	Resistance allele efficiency of ECHACR element into W
ceEH	Homing efficiency of ECHACR element into H
crEH	Resistance allele efficiency of ECHACR element into H
d1	Background mutation rate from W into R allele
d2	Background mutation rate from H into R allele
d3	Background mutation rate from E into R allele
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Details

This inheritance pattern corresponds to the Active Genetic Neutralizing Elements for Halting or Deleting Gene Drives publication.

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: ECHACRX

Description

This function creates an X-linked ECHACR construct, it has 5 alleles at the first locus and 4 alleles at the second.

• W: Wild-type

• H: Homing allele

• E: Eraser allele

• R: No-cost resistance allele

• B: Detrimental resistance allele

• cHW: Rate of homing from H, W -> H transition

• cEH: Rate of homing from E, H -> E transition

• cEW2: Rate of homing from E, W -> E transition

Usage

cubeECHACRX(
  cHW = 1,
  cEHW = 1,
  cEW1 = 1,
  cEW2 = 1,
  cEH = 1,
  chHW = 0,
  crHW = 0,
  chEHW = 0,
  crEHW = 0,
  ceEW1 = 0,
  crEW1 = 0,
  ceEW2 = 0,
  crEW2 = 0,
  ceEH = 0,
  crEH = 0,
  d1 = 0,
  d2 = 0,
  d3 = 0,
  dHW = 0,
  dEH = 0,
  dEW = 0,
  drHW = 0,
  drEH = 0,
  drEW = 0,
  crossF = 0,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cHW	Cutting efficiency of drive allele at locus 1
cEHW	Cutting efficiency of drive allele, in the presence of ECHACR element, at locus 1
cEW1	Cutting efficiency of ECHACR element into W at locus 1
cEW2	Cutting efficiency of ECHACR element into W at locus 2
cEH	Cutting efficiency of ECHACR element into H
chHW	Homing efficiency of drive allele at locus 1
crHW	Resistance allele efficiency of drive allele at locus 1
chEHW	Homing efficiency of drive allele, in the presence of ECHACR element, at locus 1
crEHW	Resistance allele efficiency of drive allele, in the presence of ECHACR element, at locus 1
ceEW1	Homing efficiency of ECHACR element into W at locus 1
crEW1	Resistance allele efficiency of ECHACR element into W at locus 1
ceEW2	Homing efficiency of ECHACR element into W at locus 2
crEW2	Resistance allele efficiency of ECHACR element into W at locus 2
ceEH	Homing efficiency of ECHACR element into H
crEH	Resistance allele efficiency of ECHACR element into H
d1	Background mutation rate from W into R allele
d2	Background mutation rate from H into R allele
d3	Background mutation rate from E into R allele
dHW	Female H deposition rate against W
dEH	Female E deposition rate against H
dEW	Female E deposition rate against W
drHW	Female resistance generation rate, from H allele
drEH	Female resistance generation rate, from E allele
drEW	Female resistance generation rate, from E allele
crossF	Female crossover rate. 0 is fully linked, 0.5 is unlinked, 1 is negatively linked
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Details

This inheritance pattern corresponds to the Active Genetic Neutralizing Elements for Halting or Deleting Gene Drives publication.

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: Homing Drive with 1 Resistance Allele

Description

This function creates an inheritance cube to model a homing gene drive (such as a CRISPR-Cas9 system) that creates 1 type of resistance allele. It assumes no sex-specific inheritance patterns and the construct is on an autosome.

Usage

cubeHoming1RA(
  c = 1,
  ch = 0,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

c	Cutting rate
ch	Successful homing rate rate
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: CRISPR (Clustered Regularly Interspaced Short Palindromic Repeats) with 2 Resistance Alleles and maternal deposition

Description

This is a sex-specific version of the original cube cubeHoming1RA. It assumes that the construct is on an autosome and there can be different male/female homing rates. It also has maternal deposition, i.e., when the male provides a W allele to a female with a H allele, some portion are cut during oogenesis. If the maternal deposition parameters are zero (d* parameters), this is a normal CRISPR drive.

Usage

cubeHomingDrive(
  cM = 1,
  cF = 1,
  dF = 0,
  chM = 0,
  crM = 0,
  chF = 0,
  crF = 0,
  dhF = 0,
  drF = 0,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cM	Male homing rate
cF	Female homing rate
dF	Female deposition homing rate
chM	Male correct homing rate
crM	Male resistance generating rate
chF	Female correct homing rate
crF	Female resistance generating rate
dhF	Female correct deposition rate
drF	Female resistance deposition rate
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: Killer-Rescue System

Description

This function creates an inheritance cube to model a Killer-Rescue system. Killer-Rescue is a 2-locus system: one locus has a toxin and the other locus contains the antidote. The loci are assumed independent and are non-homing.
This drive has 3 alleles at locus 1 and 2 alleles and locus 2:

• Locus 1

  • T: Wild-type allele

  • K: "Killer" toxin allele

  • R: Broken toxin allele

• Locus 2

  • W: Wild-type allele

  • A: Antidote allele

Usage

cubeKillerRescue(
  eR = 0,
  Keff = 1,
  Aeff = 1,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

eR	Conversion of K allele to R allele, a basal mutation rate
Keff	Toxin efficacy
Aeff	Antidote efficacy
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: MEDEA (Maternal Effect Dominant Embryonic Arrest)

Description

This function creates an inheritance cube to model a MEDEA drive system. This system was first discovered in flour beetles. It biases inheritance by expressing a maternal toxin such that offspring die unless they express a zygotic antidote.
This drive has 3 alleles at 1 locus:

• W: Wild-type allele

• M: MEDEA allele

• R: Resistance allele

Usage

cubeMEDEA(
  rM = 0,
  rW = 0,
  Teff = 1,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

rM	Breakdown of MEDEA allele, no homing/toxin/antidote, M -> R conversion
rW	De novo resistance generation, W -> R conversion
Teff	Efficacy of the toxin
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: Mendelian

Description

This function creates a Mendelian Inheritance Cube. It only handles simple, alphabetic genotypes.
The default is 3 alleles at 1 locus, this can be extended to however many alleles one is interested in, but only at 1 locus.

Usage

cubeMendelian(
  gtype = c("AA", "Aa", "aa"),
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

gtype	Vector of genotypes, with the wild-type in the first position
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Generate and Modify Default Genotype-specific Parameters

Description

This is an internal function for cubes.

Usage

cubeModifiers(
  gtype,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

gtype	character vector of genotypes
eta	genotype-specific mating fitness, handles assortative mating as well
phi	genotype-specific sex ratio at emergence
omega	genotype-specific multiplicative modifier of adult mortality
xiF	genotype-specific female pupatory success
xiM	genotype-specific male pupatory success
s	genotype-specific fractional reduction(increase) in fertility

---

Inheritance Cube: 1 Locus Maternal-Toxin/Zygotic-Antidote System

Description

This function creates a 1 locus maternal-toxin/zygotic-antidote system. This is similar to the construct called UDmel. There is no resistance generation in this model.
This drive has 3 alleles at 1 locus:

• A: Maternal-toxin 1, zygotic-antidote 2

• B: Maternal-toxin 2, zygotic-antidote 1

• W: Wild-type allele

Usage

cubeOneLocusTA(
  TAEfficacy = 1,
  TBEfficacy = 1,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

TAEfficacy	Maternal toxin A efficacy
TBEfficacy	Maternal toxin B efficacy
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: Reciprocal Translocation

Description

This function creates an inheritance cube to model a reciprocal translocation. This technology was the original form of underdominant system. It involves 2 chromosomes, each with two alleles.
This drive has 4 alleles at 2 loci:

• a: Wild-type at locus A

• A: Translocation at locus A

• b: Wile-type at locus B

• B: Translocation at locus B

Usage

cubeReciprocalTranslocations(
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: RIDL (Release of Insects with Dominant Lethality)

Description

This function creates a RIDL system. RIDL (Release of Insects with Dominant Lethality), is a form of SIT. Created by Oxitec, this is based on a positive feedback loop using the toxic tTAV gene, controlled under lab conditions by the TetO promoter. This has 2 alleles at 1 locus

• W: Wild-type allele

• R: OX513 RIDL allele

Usage

cubeRIDL(
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: Split CRISPR drive with 2 Resistance Alleles and male/female specific homing

Description

This is a sex-specific version of a split CRISPR drive. At one locus is the Cas9, inherited in a Mendelian fashion. At a second, unlinked, locus are the gRNAs. When the two loci occur together, the gRNAs drive, with potential damaged alleles, but the Cas9 remains Mendelian. It is assumed that this is an autosomal drive. This drive corresponds to the confinable gene drive system developed by the Akbari lab.

Usage

cubeSplitDrive(
  cM = 1,
  chM = 0,
  crM = 0,
  ccM = cM,
  cchM = chM,
  ccrM = crM,
  cF = 1,
  chF = 0,
  crF = 0,
  ccF = cF,
  cchF = chF,
  ccrF = crF,
  dW = 0,
  dhW = 0,
  drW = 0,
  ddW = dW,
  ddhW = dhW,
  ddrW = drW,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cM	Cutting efficiency in males, one Cas9 allele
chM	Homing efficiency in males, one Cas9 allele
crM	Resistance efficiency in males, one Cas9 allele
ccM	Cutting efficiency in males, two Cas9 alleles
cchM	Homing efficiency in males, two Cas9 alleles
ccrM	Resistance efficiency in males, two Cas9 alleles
cF	Cutting efficiency in females, one Cas9 allele
chF	Homing efficiency in females, one Cas9 allele
crF	Resistance efficiency in females, one Cas9 allele
ccF	Cutting efficiency in females, two Cas9 alleles
cchF	Homing efficiency in females, two Cas9 alleles
ccrF	Resistance efficiency in females, two Cas9 alleles
dW	Maternal deposition cutting, one Cas9 allele
dhW	Maternal deposition homing, one Cas9 allele
drW	Maternal deposition resistance, one Cas9 allele
ddW	Maternal deposition cutting, two Cas9 alleles
ddhW	Maternal deposition homing, two Cas9 alleles
ddrW	Maternal deposition resistance, two Cas9 alleles
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: tGD

Description

The trans-complementing Gene Drive (tGD) is a 1-locus, 2 target site drive. The first target site corresponds to the Cas protein, the second to an effector gene and the gRNAs. There are two sets of gRNAs, because each target site may have different cutting/homing/resistance rates, and each sex can have different rates for all of those things. Additionally, the parent that you receive your Cas from dictates its efficiency. Therefor, this construct has 5 alleles at the first locus and 4 alleles at the second.

• Locus 1

  • W: Wild-type

  • P: Paternal Cas9

  • M: Maternal Cas9

  • R: Resistant allele 1

  • B: Resistant allele 2

• Locus 2

  • W: Wild-type

  • G: gRNAs

  • R: Resistant 1

  • B: Resistant 2

Usage

cubeTGD(
  cM1 = 0,
  cM2 = 0,
  cP1 = 0,
  cP2 = 0,
  hM1 = 0,
  hM2 = 0,
  hP1 = 0,
  hP2 = 0,
  rM1 = 0,
  rM2 = 0,
  rP1 = 0,
  rP2 = 0,
  crM = 0,
  crP = 0,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cM1	Maternally inherited Cas9 cutting rate at locus 1
cM2	Maternally inherited Cas9 cutting rate at locus 2
cP1	Paternally inherited Cas9 cutting rate at locus 1
cP2	Paternally inherited Cas9 cutting rate at locus 2
hM1	Maternally inherited Cas9 homing efficiency at locus 1
hM2	Maternally inherited Cas9 homing efficiency at locus 2
hP1	Paternally inherited Cas9 homing efficiency at locus 1
hP2	Paternally inherited Cas9 homing efficiency at locus 2
rM1	Maternally inherited Cas9 resistance efficiency at locus 1
rM2	Maternally inherited Cas9 resistance efficiency at locus 2
rP1	Paternally inherited Cas9 resistance efficiency at locus 1
rP2	Paternally inherited Cas9 resistance efficiency at locus 2
crM	Maternal crossover rate, 0 is completely linked, 0.5 is unlinked, 1.0 is complete divergence
crP	Paternal crossover rate, 0 is completely linked, 0.5 is unlinked, 1.0 is complete divergence
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Details

This drive corresponds to the transcomplementing gene drive developed by the Gantz and Bier lab.

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: tGDX

Description

The trans-complementing Gene Drive (tGD) is a 1-locus, 2 target site drive. The first target site corresponds to the Cas protein, the second to an effector gene and the gRNAs. There are two sets of gRNAs, because each target site may have different cutting/homing/resistance rates, and each sex can have different rates for all of those things. Additionally, the parent that you receive your Cas from dictates its efficiency. Therefor, this construct has 6 alleles at the first locus and 5 alleles at the second.

• Locus 1

  • W: Wild-type

  • P: Paternal Cas9

  • M: Maternal Cas9

  • R: Resistant allele 1

  • B: Resistant allele 2

  • Y: Y allele

• Locus 2

  • W: Wild-type

  • G: gRNAs

  • R: Resistant 1

  • B: Resistant 2

  • Y: Y allele

Usage

cubeTGDX(
  cM1 = 0,
  cM2 = 0,
  cP1 = 0,
  cP2 = 0,
  hM1 = 0,
  hM2 = 0,
  hP1 = 0,
  hP2 = 0,
  rM1 = 0,
  rM2 = 0,
  rP1 = 0,
  rP2 = 0,
  crM = 0,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

cM1	Maternally inherited Cas9 cutting rate at locus 1
cM2	Maternally inherited Cas9 cutting rate at locus 2
cP1	Paternally inherited Cas9 cutting rate at locus 1
cP2	Paternally inherited Cas9 cutting rate at locus 2
hM1	Maternally inherited Cas9 homing efficiency at locus 1
hM2	Maternally inherited Cas9 homing efficiency at locus 2
hP1	Paternally inherited Cas9 homing efficiency at locus 1
hP2	Paternally inherited Cas9 homing efficiency at locus 2
rM1	Maternally inherited Cas9 resistance efficiency at locus 1
rM2	Maternally inherited Cas9 resistance efficiency at locus 2
rP1	Paternally inherited Cas9 resistance efficiency at locus 1
rP2	Paternally inherited Cas9 resistance efficiency at locus 2
crM	Maternal crossover rate, 0 is completely linked, 0.5 is unlinked, 1.0 is complete divergence
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Details

This drive corresponds to the transcomplementing gene drive developed by the Gantz and Bier lab.

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: 2 Locus Maternal-Toxin/Zygotic-Antidote System

Description

This function creates a 2 locus maternal-toxin/zygotic-antidote system. This is similar to the construct called UDmel. There is no resistance generation in this model.
This drive has 2 unlinked alleles, 1 allele each at 2 loci:

• A: Maternal-toxin 1, zygotic-antidote 2

• a: Wild-type at locus A

• B: Maternal-toxin 2, zygotic-antidote 1

• b: Wild-type at locus B

Usage

cubeTwoLocusTA(
  TAEfficacy = 1,
  TBEfficacy = 1,
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

TAEfficacy	Maternal toxin A efficacy
TBEfficacy	Maternal toxin B efficacy
eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Inheritance Cube: Wolbachia

Description

This function creates an inheritance cube to model a Wolbachia infection. Wolbachia is a parasite that can infect mosquitoes. It biases its inheritance through cytoplasmic incompatibility.
This drive has 2 alleles at 1 locus:

• W: has Wolbachia

• w: does not have Wolbachia

Usage

cubeWolbachia(
  eta = NULL,
  phi = NULL,
  omega = NULL,
  xiF = NULL,
  xiM = NULL,
  s = NULL
)

Arguments

eta	Genotype-specific mating fitness
phi	Genotype-specific sex ratio at emergence
omega	Genotype-specific multiplicative modifier of adult mortality
xiF	Genotype-specific female pupatory success
xiM	Genotype-specific male pupatory success
s	Genotype-specific fractional reduction(increase) in fertility

Details

Cytoplasmic Incompatibility:

• male W cross female w -> all offspring die (complete penetrance)

• male w cross female W -> all offspring inherit Wolbachia

Value

Named list containing the inheritance cube, transition matrix, genotypes, wild-type allele, and all genotype-specific parameters.

---

Erase all files in a directory

Description

Given a directory path, check that it exists, and if so, delete all its contents.

Usage

eraseDirectory(directory, verbose = TRUE)

Arguments

directory	Directory whose contents will be deleted
verbose	Chatty? Default is TRUE

Examples

## Not run: 
# Path to directory, can tilde expand
myPath <- "~/path/to/write/output"

# Erase directory
#  No return value
eraseDirectory(directory = myPath)

## End(Not run)

---

Make List of Modified Mosquito Releases

Description

Sets up a release schedule for a single patch, calls basicRepeatedReleases internally.

Usage

generateReleaseVector(driveCube, releasesParameters, nameGenotypes = NULL)

Arguments

driveCube	Gene-drive cube
releasesParameters	A list containing the releasesStart, releasesNumber releasesInterval, and releaseProportion named values.
nameGenotypes	Optional list to specify different genotypes for egg/male/female releases. This is required for mated female releases. This parameter overrides the default release type.

Examples

# setup a drive cube, using Mendelian as the example
cube <- cubeMendelian()

# setup release parameter list
#  releasesStart is the time of first release
#  releasesNumber is the number of releases
#  releasesInterval is the number of days between releases
#  releaseProportion is the number of mosquitoes released
relParams <- list(releasesStart = 25, releasesNumber = 1,
                  releasesInterval = 0, releaseProportion = 10)

# generate male releases
mRelVec <- generateReleaseVector(driveCube = cube,
                                 releasesParameters = relParams)

# generate mated female releases
fRelVec <- generateReleaseVector(driveCube = cube,
                                 releasesParameters = relParams,
                                 nameGenotypes = list(c("AA","AA", 10),
                                                      c("AA","aa", 10)))

---

Get alpha

Description

Return density dependent mortality, see calcDensityDependentDeathRate

Usage

get_alpha_Network(ix)

Arguments

ix	Index of patch

---

Get beta

Description

Return size of wild-type egg batch

Usage

get_beta_Network()

---

Get conADF

Description

Return connection where adult female dynamics are written to

Usage

get_conF_Network()

---

Get conADM

Description

Return connection where adult male dynamics are written to

Usage

get_conM_Network()

---

Get Element(s) of Drive Cube by Index

Description

Return elements or slices of drive cube. If all NULL return entire cube.

Usage

get_drivecubeindex_Network(fG = NULL, mG = NULL, oG = NULL)

Arguments

fG	Female genotype index
mG	Male genotype index
oG	Offspring genotype index

---

Get eta

Description

Get eta

Usage

get_eta_Network(fIdx)

Arguments

fIdx	Index of female genotype to pull Return genotype-specific mating fitness

---

Get female Population

Description

Return females (nGenotypes X nGenotypes matrix)

Usage

get_femalePop_Patch()

---

Get genotypesID

Description

Return character vector of possible genotypes

Usage

get_genotypesID_Network()

---

Get genotypesN

Description

Return number of possible genotypes

Usage

get_genotypesN_Network()

---

Get male Population

Description

Return males (nGenotypes vector)

Usage

get_malePop_Patch()

---

Get muAd

Description

Return adult mortality

Usage

get_muAd_Network()

---

Get muAq

Description

Return larval mortality, see calcLarvalStageMortalityRate

Usage

get_muAq_Network()

---

Get nPatch

Description

Return number of patches

Usage

get_nPatch_Network()

---

Get omega

Description

Return genotype-specific multiplicative modifier of adult mortality

Usage

get_omega_Network()

---

Get Patch Release Schedule

Description

Return the release schedule for a patch for male or female

Usage

get_patchReleases_Network(patch, sex = "M")

Arguments

patch	Index of patch
sex	Character in 'M', 'F', 'Egg', 'mF'

---

Get phi

Description

Return genotype-specific sex ratio at emergence

Usage

get_phi_Network()

---

Get s

Description

Return genotype-specific fractional reduction(increase) in fertility

Usage

get_s_Network()

---

Get Female Viability Mask (tau)

Description

Get Female Viability Mask (tau)

Usage

get_tau_Network(fG = NULL, mG = NULL, oG = NULL)

Arguments

fG	Number for which female genotype to get
mG	Number for which male genotype to get
oG	Number for which offspring genotype to get Return matrix

---

Get timeAq

Description

Return duration of aquatic stages.

Usage

get_timeAq_Network(stage = NULL)

Arguments

stage

Character in 'E', 'L', 'P'; if NULL return total duration

---

Get tNow

Description

Return current simulation time

Usage

get_tNow_Network()

---

Get xiF

Description

Return genotype-specific female pupatory success

Usage

get_xiF_Network()

---

Get xiM

Description

Return genotype-specific male pupatory success

Usage

get_xiM_Network()

---

Utility to Imitate ggplot2 Colors

Description

Sample at equally spaced intervals along the color wheel

Usage

ggColUtility(n, alpha = 0.75)

Arguments

n	Number of colors
alpha	Transparency

---

Kernels Parameters

Description

A named list containing maximum likelihood fitted parameter values from mosquito dispersal estimates.

Usage

data(kernels)

Format

named list with 5 elements:

lnorm_mean

log mean of log-normal density

lnorm_sd

log standard deviation of log-normal density

gamma_shape

shape parameter of gamma density

gamma_sd

rate parameter of gamma density

exp_rate

rate parameter of exponential density

---

MGDrivE: Mosquito Gene Drive Explorer

Description

MGDrivE: Mosquito Gene Drive Explorer

Introduction

Recent developments of CRISPR-Cas9 based homing endonuclease gene drive systems, for the suppression or replacement of mosquito populations, have generated much interest in their use for control of mosquito-borne diseases (such as dengue, malaria, Chikungunya and Zika). This is because genetic control of pathogen transmission may complement or even substitute traditional vector-control interventions, which have had limited success in bringing the spread of these diseases to a halt. Despite excitement for the use of gene drives for mosquito control, current modeling efforts have analyzed only a handful of these new approaches (usually studying just one per framework). Moreover, these models usually consider well-mixed populations with no explicit spatial dynamics. To this end, we are developing MGDrivE (Mosquito Gene DRIVe Explorer), in cooperation with the 'UCI Malaria Elimination Initiative', as a flexible modeling framework to evaluate a variety of drive systems in spatial networks of mosquito populations. This framework provides a reliable testbed to evaluate and optimize the efficacy of gene drive mosquito releases. What separates MGDrivE from other models is the incorporation of mathematical and computational mechanisms to simulate a wide array of inheritance-based technologies within the same, coherent set of equations. We do this by treating the population dynamics, genetic inheritance operations, and migration between habitats as separate processes coupled together through the use of mathematical tensor operations. This way we can conveniently swap inheritance patterns whilst still making use of the same set of population dynamics equations. This is a crucial advantage of our system, as it allows other research groups to test their ideas without developing new models and without the need to spend time adapting other frameworks to suit their needs.

Brief Description

MGDrivE is based on the idea that we can decouple the genotype inheritance process from the population dynamics equations. This allows the system to be treated and developed in three semi-independent modules that come together to form the system. The way this is done will be described later in this document but a reference diagram is shown here.

Previous Work

The original version of this model was based on work by (Deredec et al. 2011; Hancock and Godfray 2007) and adapted to accommodate CRISPR homing dynamics in a previous publication by our team (Marshall et al. 2017). As it was described, we extended this framework to be able to handle a variable number of genotypes, and migration across spatial scenarios. We accomplish this by adapting the equations to work in a tensor-oriented manner, where each genotype can have different processes affecting their particular strain (death rates, mating fitness, sex-ratio bias, et cetera).

Notation and Conventions

Before beginning the full description of the model we will define some of the conventions we followed for the notation of the written description of the system.

• Overlines are used to denote the dimension of a tensor.

• Subscript brackets are used to indicate an element in time. For example: $L_{[t-1]}$ is the larval population at time: $t-1$.

• Parentheses are used to indicate the parameter(s) of a function. For example: $\overline{O(T_{e}+T_{l})}$ represents the function $O$ evaluated with the parameter: $T_{e}+T_{l}$

• Matrices follow a 'row-first' indexing order (i: row, j: column)

In the case of one dimensional tensors, each slot represents a genotype of the population. For example, the male population is stored in the following way:

$\overline{Am} = \left(\begin{array}{c} g_1 \\ g_2 \\ g_3 \\ \vdots \\ g_n \end{array}\right) _{i}$

All the processes that affect mosquitoes in a genotype-specific way are defined and stored in this way within the framework.

There are two tensors of squared dimensionality in the model: the adult females matrix, and the genotype-specific male-mating ability ($\overline{\eta}$) In the case of the former the rows represent the females' genotype, whilst the columns represent the genotype of the male they mated with:

$\overline{\overline{Af}} = \left(\begin{array}{ccccc} g_{11} & g_{12} & g_{13} & \cdots & g_{1n}\\ g_{21} & g_{22} & g_{23} & \cdots & g_{2n}\\ g_{31} & g_{32} & g_{33} & \cdots & g_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ g_{n1} & g_{n2} & g_{n3} & \cdots & g_{nn} \end{array}\right) _{ij}$

The genotype-specific male mating ability, on the other hand, stores the females' genotype in the rows, and the male genotypes in the columns of the matrix.

References

Deredec A, Godfray HCJ, Burt A (2011). “Requirements for effective malaria control with homing endonuclease genes.” Proceedings of the National Academy of Sciences of the United States of America, 108(43), E874–80. ISSN 1091-6490, doi:10.1073/pnas.1110717108, https://www.pnas.org/content/108/43/E874.

Hancock PA, Godfray HCJ (2007). “Application of the lumped age-class technique to studying the dynamics of malaria-mosquito-human interactions.” Malaria journal, 6, 98. ISSN 1475-2875, doi:10.1186/1475-2875-6-98, https://malariajournal.biomedcentral.com/articles/10.1186/1475-2875-6-98.

Marshall J, Buchman A, C. HMS, Akbari OS (2017). “Overcoming evolved resistance to population-suppressing homing-based gene drives.” Nature Scientific Reports, 1–46. ISSN 2045-2322, doi:10.1038/s41598-017-02744-7, https://www.nature.com/articles/s41598-017-02744-7.

---

MGDrivE: Inheritance Cube

Description

To model an arbitrary number of genotypes efficiently in the same mathematical framework, we use a 3-dimensional array structure (cube) where each axis represents the following information:

• x: female adult mate genotype

• y: male adult mate genotype

• z: proportion of the offspring that inherits a given genotype (layer)

Details

The cube structure gives us the flexibility to apply tensor operations to the elements within our equations, so that we can calculate the stratified population dynamics rapidly; and within a readable, flexible computational framework. This becomes apparent when we define the equation we use for the computation of eggs laid at any given point in time:

$\overline{O(T_x)} = \sum_{j=1}^{n} \Bigg( \bigg( (\beta*\overline{s} * \overline{ \overline{Af_{[t-T_x]}}}) * \overline{\overline{\overline{Ih}}} \bigg) * \Lambda \Bigg)^{\top}_{ij}$

In this equation, the matrix containing the number of mated adult females ($\overline{\overline{Af}}$) is multiplied element-wise with each one of the layers containing the eggs genotypes proportions expected from this cross ($\overline{\overline{\overline{Ih}}}$). The resulting matrix is then multiplied by a binary 'viability mask' ($\Lambda$) that filters out female-parent to offspring genetic combinations that are not viable due to biological impediments (such as cytoplasmic incompatibility). The summation of the transposed resulting matrix returns us the total fraction of eggs resulting from all the male to female genotype crosses ($\overline{O(T_x)}$).

Note: For inheritance operations to be consistent within the framework the summation of each element in the z-axis (this is, the proportions of each one of the offspring's genotypes) must be equal to one.

Drive-specific Cubes

An inheritance cube in an array object that specifies inheritance probabilities (offspring genotype probability) stratified by male and female parent genotypes. MGDrivE provides the following cubes to model different gene drive systems:

• cubeOneLocusTA: 1 Locus Maternal-Toxin/Zygotic-Antidote System

• cubeTwoLocusTA: 2 Locus Maternal-Toxin/Zygotic-Antidote System

• cubeClvR: 1 Locus Cleave and Rescue (ClvR)

• cubeClvR2: 2 Locus Cleave and Rescue (ClvR)

• cubeHoming1RA: Homing Drive with 1 Resistance Allele

• cubeHomingDrive: CRISPR (Clustered Regularly Interspaced Short Palindromic Repeats) with 2 Resistance Allele

• cubeECHACR: ERACR-CHACR

• cubeECHACRX: ECHACR, X-Linked

• cubeKillerRescue: Killer-Rescue System

• cubeMEDEA: MEDEA (Maternal Effect Dominant Embryonic Arrest)

• cubeReciprocalTranslocations: Reciprocal Translocation

• cubeRIDL: RIDL (Release of Insects with Dominant Lethality)

• cubeMendelian: Mendelian

• cubeSplitDrive: Split CRISPR drive

• cubeTGD: trans-complementing Gene Drive

• cubeTGDX: trans-complementing Gene Drive, X-Linked

• cubeWolbachia: Wolbachia

Functions for Cubes

We provide one auxiliary function to operate on cube objects.

• cube2csv: Export slices of a cube to .csv format

---

MGDrivE: Model's Mathematical Description

Description

The original version of this model was based on work by (Deredec et al. 2011; Hancock and Godfray 2007) and adapted to accommodate CRISPR homing dynamics in a previous publication by our team (Marshall et al. 2017). As it was described, we extended this framework to be able to handle a variable number of genotypes, and migration across spatial scenarios. We did this by adapting the equations to work in a tensor-oriented manner, where each genotype can have different processes affecting their particular strain (death rates, mating fitness, sex-ratio bias, et cetera).

Inheritance Cube and Oviposition

To allow the extension of the framework to an arbitrary number of genotypes, we transformed traditional inheritance matrices into inheritance cubes, where each of the axis represents the following information:

• x: female adult mate genotype

• y: male adult mate genotype

• z: proportion of the offspring that inherits a given genotype (slice)

The 'cube' structure gives us the flexibility to apply tensor operations to the elements within our equations, so that we can calculate the stratified population dynamics rapidly; and within a readable, flexible computational framework. This becomes apparent when we define the equation we use for the computation of eggs laid at any given point in time:

$\overline{O(T_x)} = \sum_{j=1}^{n} \Bigg( \bigg( (\beta*\overline{s} * \overline{ \overline{Af_{[t-T_x]}}}) * \overline{\overline{\overline{Ih}}} \bigg) * \Lambda \Bigg)^{\top}_{ij}$

In this equation, the matrix containing the number of mated adult females ($\overline{\overline{Af}}$) is multiplied element-wise with each one of the slices containing the eggs genotypes proportions expected from this cross ($\overline{\overline{\overline{Ih}}}$). The resulting matrix is then multiplied by a binary 'viability mask' ($\Lambda$) that filters out female-parent to offspring genetic combinations that are not viable due to biological impediments (such as cytoplasmic incompatibility). The summation of the transposed resulting matrix returns us the total fraction of eggs resulting from all the male to female genotype crosses ($\overline{O(T_{x})}$).

Note: For inheritance operations to be consistent within the framework, the summation of each element in the 'z' axis (this is, the proportions of each one of the offspring's genotypes) must be equal to one.

Population Dynamics

During the three aquatic stages, a density-independent mortality process takes place:

$\theta_{st}=(1-\mu_{st})^{T_{st}}$

Along with a density dependent process dependent on the number of larvae in the environment:

$F(L[t])=\Bigg(\frac{\alpha}{\alpha+\sum{\overline{L[t]}}}\Bigg)^{1/T_l}$

where $\alpha$ represents the strength of the density-dependent process. This parameter is calculated with:

$\alpha=\Bigg( \frac{1/2 * \beta * \theta_e * Ad_{eq}}{R_m-1} \Bigg) * \Bigg( \frac{1-(\theta_l / R_m)}{1-(\theta_l / R_m)^{1/T_l}} \Bigg)$

in which $\beta$ is the species' fertility in the absence of gene-drives, $Ad_{eq}$ is the adult mosquito population equilibrium size, and $R_{m}$ is the population growth in the absence of density-dependent mortality. This population growth is calculated with the average generation time ($g$), the adult mortality rate ($\mu_{ad}$), and the daily population growth rate ($r_{m}$):

$g=T_{e}+T_{l}+T_{p}+\frac{1}{\mu_{ad}}\\R_{m}=(r_{m})^{g}$

Larval Stages

The computation of the larval stage in the population is crucial to the model because the density dependent processes necessary for equilibrium trajectories to be calculated occur here. This calculation is performed with the following equation:

$D(\theta_l,T_x) = \begin{array}{ll} \theta_{l[0]}^{'}=\theta_l & \quad i = 0 \\ \theta_{l[i+1]}^{'} = \theta_{l[i]}^{'} *F(\overline{L_{[t-i-T_x]}}) & \quad i \leq T_l \end{array}$

In addition to this, we need the larval mortality ($\mu_{l}$):

$\mu_{l}=1-\Bigg( \frac{R_{m} * \mu_{ad}}{1/2 * \beta * (1-\mu_{m})} \Bigg)^{\frac{1}{T_{e}+T_{l}+T_{p}}}$

With these mortality processes, we are now able to calculate the larval population:

$\overline{L_{[t]}}= \overline{L_{[t-1]}} * (1-\mu_{l}) * F(\overline{L_{[t-1]})}\\ +\overline{O(T_{e})}* \theta_{e} \\ - \overline{O(T_{e}+T_{l})} * \theta_{e} * D(\theta_{l},0)$

where the first term accounts for larvae surviving one day to the other; the second term accounts for the eggs that have hatched within the same period of time; and the last term computes the number of larvae that have transformed into pupae.

Adult Stages

We are ultimately interested in calculating how many adults of each genotype exist at any given point in time. For this, we first calculate the number of eggs that are laid and survive to the adult stages with the equation:

$\overline{E^{'}}= \overline{O(T_{e}+T_{l}+T_{p})} \\ * \bigg(\overline{\xi_{m}} * (\theta_{e} * \theta_{p}) * (1-\mu_{ad}) * D(\theta_{l},T_{p}) \bigg)$

With this information we can calculate the current number of male adults in the population by computing the following equation:

$\overline{Am_{[t]}}= \overline{Am_{[t-1]}} * (1-\mu_{ad})*\overline{\omega_{m}}\\ + (1-\overline{\phi}) * \overline{E^{'}}\\ + \overline{\nu m_{[t-1]}}$

in which the first term represents the number of males surviving from one day to the next; the second one, the fraction of males that survive to adulthood ($\overline{E'}$) and emerge as males ($1-\phi$); the last one is used to add males into the population as part of gene-drive release campaigns.

Female adult populations are calculated in a similar way:

$\overline{\overline{Af_{[t]}}}= \overline{\overline{Af_{[t-1]}}} * (1-\mu_{ad}) * \overline{\omega_{f}}\\ + \bigg( \overline{\phi} * \overline{E^{'}}+\overline{\nu f_{[t-1]}}\bigg)^{\top} * \bigg( \frac{\overline{\overline{\eta}}*\overline{Am_{[t-1]}}}{\sum{\overline{Am_{[t-1]}}}} \bigg)$

where we first compute the surviving female adults from one day to the next; and then we calculate the mating composition of the female fraction emerging from pupa stage. To do this, we obtain the surviving fraction of eggs that survive to adulthood ($\overline{E'}$) and emerge as females ($\phi$), we then add the new females added as a result of gene-drive releases ($\overline{\nu f_{[t-1]}}$). After doing this, we calculate the proportion of males that are allocated to each female genotype, taking into account their respective mating fitnesses ($\overline{\overline{\eta}}$) so that we can introduce the new adult females into the population pool.

Gene Drive Releases and Effects

As it was briefly mentioned before, we are including the option to release both male and/or female individuals into the populations. Another important t hing to emphasize is that we allow flexible releases sizes and schedules. Our ] model handles releases internally as lists of populations compositions so, it is possible to have releases performed at irregular intervals and with different numbers of mosquito genetic compositions as long as no new genotypes are introduced (which have not been previously defined in the inheritance cube).

$\overline{\nu} = \bigg\{ \left(\begin{array}{c} g_1 \\ g_2 \\ g_3 \\ \vdots \\ g_n \end{array}\right)_{t=1} , \left(\begin{array}{c} g_1 \\ g_2 \\ g_3 \\ \vdots \\ g_n \end{array}\right)_{t=2} , \cdots , \left(\begin{array}{c} g_1 \\ g_2 \\ g_3 \\ \vdots \\ g_n \end{array}\right)_{t=x} \bigg\}$

So far, however, we have not described the way in which the effects of these gene-drives are included into the mosquito populations dynamics. This is done through the use of various modifiers included in the equations:

• $\overline{\omega}$: Relative increase in mortality (zero being full mortality effects and one no mortality effect)

• $\overline{\phi}$: Relative shift in the sex of the pupating mosquitoes (zero biases the sex ratio towards males, whilst 1 biases the ratio towards females).

• $\overline{\overline{\eta}}$: Standardized mating fitness (zero being complete fitness ineptitude, and one being regular mating skills).

• $\overline{\beta}$: Fecundity (average number of eggs laid).

• $\overline{\xi}$: Pupation success (zero being full mortality and one full pupation success).

Migration

To simulate migration within our framework we are considering patches (or nodes) of fully-mixed populations in a network structure. This allows us to handle mosquito movement across spatially-distributed populations with a transitions matrix, which is calculated with the tensor outer product of the genotypes populations tensors and the transitions matrix of the network as follows:

$\overline{Am_{(t)}^{i}}= \sum{\overline{A_{m}^j} \otimes \overline{\overline{\tau m_{[t-1]}}}} \\ \overline{\overline{Af_{(t)}^{i}}}= \sum{\overline{\overline{A_{f}^j}} \otimes \overline{\overline{\tau f_{[t-1]}}}}$

In these equations the new population of the patch $i$ is calculated by summing the migrating mosquitoes of all the $j$ patches across the network defined by the transitions matrix $\tau$, which stores the mosquito migration probabilities from patch to patch. It is worth noting that the migration probabilities matrices can be different for males and females; and that there's no inherent need for them to be static (the migration probabilities may vary over time to accommodate wind changes due to seasonality).

Parameters

This table compiles all the parameters required to run MGDrivE clustered in six categories:

• Life Stages: These deal with the structure of mosquito population.

• Bionomics: This set of parameters is related to the behavior of the specific mosquito species being modeled.

• Gene Drive: Genotype-specific vectors of parameters that affect how each gene-drive modifies the responses of populations to them.

• Releases: List of vectors that control the release of genetically-modified mosquitoes.

• Population: General mosquito-population parameters that control environmentally-determined variables.

• Network: Related to migration between nodes of population units

Stochasticity

MGDrivE allows stochasticity to be included in the dynamics of various processes; in an effort to simulate processes that affect various stages of mosquitoes lives. In the next section, we will describe all the stochastic processes that can be activated in the program. It should be noted that all of these can be turned on and off independently from one another as required by the researcher.

Mosquito Biology

Oviposition

Stochastic egg laying by female/male pairs is separated into two steps: calculating the number of eggs laid by the females and then distributing laid eggs according to their genotypes. The number of eggs laid follows a Poisson distribution conditioned on the number of female/male pairs and the fertility of each female.

$Poisson( \lambda = numFemales*Fertility)$

Multinomial sampling, conditioned on the number of offspring and the relative viability of each genotype, determines the genotypes of the offspring.

$Multinomial \left(numOffspring, p_1, p_2\dots p_b \right)=\frac{numOffspring!}{p_1!\,p_2\,\dots p_n}p_1^{n_1}p_2^{n_2}\dots p_n^{n_n}$

Sex Determination

Sex of the offspring is determined by multinomial sampling. This is conditioned on the number of eggs that live to hatching and a probability of being female, allowing the user to design systems that skew the sex ratio of the offspring through reproductive mechanisms.

$Multinomial(numHatchingEggs, p_{female}, p_{female})$

Mating Stochastic mating is determined by multinomial sampling conditioned on the number of males and their fitness. It is assumed that females mate only once in their life, therefore each female will sample from the available males and be done, while the males are free to potentially mate with multiple females. The males' ability to mate is modulated with a fitness term, thereby allowing some genotypes to be less fit than others (as seen often with lab releases).

$Multinomial(numFemales, p_1f_1, p_2f_2, \dots p_nf_n)$

Hatching

Other Stochastic Processes All remaining stochastic processes (larval survival, hatching , pupating, surviving to adult hood) are determined by multinomial sampling conditioned on factors affecting the current life stage. These factors are determined empirically from mosquito population data.

Migration

Variance of stochastic movement (not used in diffusion model of migration).

References

Deredec A, Godfray HCJ, Burt A (2011). “Requirements for effective malaria control with homing endonuclease genes.” Proceedings of the National Academy of Sciences of the United States of America, 108(43), E874–80. ISSN 1091-6490, doi:10.1073/pnas.1110717108, https://www.pnas.org/content/108/43/E874.

Hancock PA, Godfray HCJ (2007). “Application of the lumped age-class technique to studying the dynamics of malaria-mosquito-human interactions.” Malaria journal, 6, 98. ISSN 1475-2875, doi:10.1186/1475-2875-6-98, https://malariajournal.biomedcentral.com/articles/10.1186/1475-2875-6-98.

Marshall J, Buchman A, C. HMS, Akbari OS (2017). “Overcoming evolved resistance to population-suppressing homing-based gene drives.” Nature Scientific Reports, 1–46. ISSN 2045-2322, doi:10.1038/s41598-017-02744-7, https://www.nature.com/articles/s41598-017-02744-7.

---

Movement Matrix: All 2

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatAll2)

Format

A matrix with 3 rows and 3 columns:

Patches 1 and 3 are sources for patch 2, which is a sink.

---

Movement Matrix: Cascade 3

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatCascade3)

Format

A matrix with 3 rows and 3 columns:

Mosquitoes in patch 1 have equal probability to stay or move to 2; mosquitoes in patch 2 have equal probability to stay or move to 3; mosquitoes in patch 3 stay there.

---

Movement Matrix: Diagonal

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatDiag)

Format

A matrix with 3 rows and 3 columns:

3 independent patches.

---

Movement Matrix: Diagonal One City

Description

A movement matrix for simulation with 1 patch.

Usage

data(moveMatDiagOneCity)

Format

A matrix with 1 rows and 1 columns:

A 1 by 1 matrix with entry 1.

---

Movement Matrix: Die

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatDie)

Format

A matrix with 3 rows and 3 columns:

All entries of matrix are 0 for testing that all mosquitoes will be killed.

---

Movement Matrix: Independent 3

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatIndependent3)

Format

A matrix with 3 rows and 3 columns:

Mosquitoes in patch 1 stay with probability 0.975, move to patch 2 with probability 0.025, mosquitoes in patch 2 and 3 stay in their patches.

---

Movement Matrix: Mixed Spill

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatMixedSpil)

Format

A matrix with 3 rows and 3 columns:

Mosquitoes in patch 1 stay with probability 0.999, move to patch 2 with probability 0.001, mosquitoes in patch 2 and 3 stay in their patches.

---

Movement Matrix: Tale of Two Cities

Description

A movement matrix for simulation with 2 patches.

Usage

data(moveMatTaleOfTwoCities)

Format

A matrix with 2 rows and 2 columns:

Mosquitoes do not move between the two patches.

---

Movement Matrix: Tri-diagonal

Description

A movement matrix for simulation with 12 patches.

Usage

data(moveMatTriDiagonal)

Format

A matrix with 12 rows and 12 columns:

Tri-diagonal matrix with approximately 0.985 probability on diagonal and rest of probability mass on k-1 and k+1 off-diagonal elements.

---

Movement Matrix: Triple

Description

A movement matrix for simulation with 3 patches.

Usage

data(moveMatTriple)

Format

A matrix with 3 rows and 3 columns:

All entries of matrix are 1 for testing that mosquitoes will be produced.

---

Run Simulation

Description

Run multiple simulations on this network

Usage

multRun_Network(verbose = TRUE)

Arguments

verbose

Chatty? Default is TRUE

---