## ---- include = FALSE--------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", echo = TRUE, hold = TRUE, fig.width = 7, fig.height = 4, eval = TRUE ) # set seed for reproducibility set.seed(seed = 42) ## ----------------------------------------------------------------------------- # setup movement matrix for 1 node moveMat <- matrix(data = 1, nrow = 1, ncol = 1) moveMat ## ----------------------------------------------------------------------------- # setup movement matrix for 2 nodes #################### # 2 nodes, no migration #################### moveMat <- matrix(data = c(1,0,0,1), nrow = 2, ncol = 2, byrow = TRUE) moveMat #################### # 2 nodes, with migration #################### # 5% migration per day from population 1 # 10% migraton per day from population 2 # Notice that the rows sum to 1 moveMat <- matrix(data = c(0.95, 0.05, 0.10, 0.90), nrow = 2, ncol = 2, byrow = TRUE) moveMat ## ----------------------------------------------------------------------------- # setup random movement matrix for 5 nodes #################### # 5 nodes #################### nNodes <- 5 # fill with random data moveMat <- matrix(data = runif(n = nNodes*nNodes), nrow = nNodes, ncol = nNodes) # normalize moveMat <- moveMat/rowSums(x = moveMat) moveMat ## ----------------------------------------------------------------------------- # setup line with 10 nodes #################### # 10 nodes in a line #################### nNodes <- 10 # define function for use triDiag <- function(upper, lower){ # return matrix retMat <- matrix(data = 0, nrow = length(upper) + 1, ncol = length(upper) + 1) # set index values for upper/lower triangles indx <- 1:length(upper) # set forward/backward migration using matrix access retMat[cbind(indx+1,indx)] <- lower retMat[cbind(indx,indx+1)] <- upper # set stay probs diag(x = retMat) <- 1-rowSums(x = retMat) return(retMat) } # fill movement matrix # Remember, rows need to sum to 1. moveMat <- triDiag(upper = rep.int(x = 0.05, times = nNodes-1), lower = rep.int(x = 0.05, times = nNodes-1)) moveMat ## ----------------------------------------------------------------------------- # realistic landscape # matrix of coordinates as latitude/longitude pairs lat_longs <- matrix(data = c(37.873507, -122.268181, 37.873578, -122.254430, 37.869806, -122.267639), nrow = 3, ncol = 2, byrow = TRUE, dimnames = list(NULL, c('Lat','Lon'))) # calculate distance matrix between points # dmat <- MGDrivE::calcHaversine(latLongs = lat_longs) # dmat <- MGDrivE::calcVinSph(latLongs = lat_longs) distMat <- MGDrivE::calcVinEll(latLongs = lat_longs) # calculate a zero-inflated movement kernal over the distances # p0 is the probability, per day, that a mosquito doesn't move. # This is the value used in Code sample 1 from the paper, and in the examples in our # github repository. # rate is the average migration rate per day, implying 1/rate is the average # migration distance. The average distance was estimated as ~55.5 meters per day, # which is the value used in Code sample 1 and in the examples on github. p0 <- 0.991 rate <- 1/55.5 moveMat <- MGDrivE::calcHurdleExpKernel(distMat = distMat, rate = rate, p0 = p0) moveMat ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### outFolder <- "mgdrive" #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365 # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 1-node network where mosquitoes do not leave moveMat <- matrix(data = 1, nrow = 1, ncol = 1) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) parameters cube <- cubeMendelian() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur every day, for 1 day. # There are 10 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=1, releasesInterval=0, releaseProportion=10) # generate release vector releasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- releasesVector patchReleases[[1]]$femaleReleases <- releasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run deterministic setupMGDrivE(stochasticityON = FALSE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we havee a network of 1 population. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, sampTime = 1, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=outFolder, verbose=FALSE) # run simulation MGDrivESim$oneRun(verbose = FALSE) #################### # Post Analysis #################### # split output by patch # Required for plotting later splitOutput(readDir = outFolder, remFile = TRUE, verbose = FALSE) # aggregate females by their mate choice # This reduces the female file to have the same columns as the male file aggregateFemales(readDir = outFolder, genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) # plot output to see effect plotMGDrivESingle(readDir = outFolder, totalPop = TRUE, lwd = 3.5, alpha = 1) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Theory Comparison #################### # read in simulation files totPop <- read.csv(file = file.path(outFolder, "M_Run001_Patch001.csv"), header = TRUE, sep = ",")[ ,-1] + read.csv(file = file.path(outFolder, "F_Aggregate_Run001_Patch001.csv"), header = TRUE, sep = ",")[ ,-1] ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### outFolder <- "mgdrive" #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365*2 # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 1-node network where mosquitoes do not leave moveMat <- as.matrix(1) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) patchPops <- rep(adultPopEquilibrium,sitesNumber) #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) parameters # This time, lets put fitness cost on the homozygotes, giving the heterozygote # an advantage # These genotypes correspond to ones in the cube. Look at a base cube first, # then set this. # homozygotes are 60% as fit as heterozygote over their entire lifetime # Since omega is the adult daily death rate, we use the built-in function to # calculate our desired lifetime cost as applied daily dayOmega <- calcOmega(mu = bioParameters$muAd, lifespanReduction = 0.60) omegaNew <- c("AA"=dayOmega, "aa"=dayOmega) # setup cube cube <- cubeMendelian(omega = omegaNew) #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 100, occur every day, for 5 days. # There are 50 mosquitoes released every time. releasesParameters <- list(releasesStart=100, releasesNumber=5, releasesInterval=0, releaseProportion=50) # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run deterministic setupMGDrivE(stochasticityON = FALSE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 1 population. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=patchPops, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=outFolder, verbose=FALSE) # run simulation MGDrivESim$oneRun(verbose = FALSE) #################### # Post Analysis #################### # split output by patch # Required for plotting later splitOutput(readDir = outFolder, remFile = TRUE, verbose = FALSE) # aggregate females by their mate choice # This reduces the female file to have the same columns as the male file aggregateFemales(readDir = outFolder, genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) # plot output to see effect plotMGDrivESingle(readDir = outFolder, totalPop = TRUE, lwd = 3.5, alpha = 1) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Theory Comparison #################### # read in simulation files totPop <- read.csv(file = file.path(outFolder, "M_Run001_Patch001.csv"), header = TRUE, sep = ",")[ ,-1] + read.csv(file = file.path(outFolder, "F_Aggregate_Run001_Patch001.csv"), header = TRUE, sep = ",")[ ,-1] ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### # directory # This is slightly obtuse for vignette building reasons # Really, all you need is a base directory, then the repetitions in folders inside that. # Here, our base directory is "mgdrive", and the repetition folders are "001","002", etc. # So, the final structure is "mgdrive/001","mgdrive/002", etc. outFolder <- "mgdrive" dir.create(path = outFolder) #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365*2 # number of Monte Carlo iterations nRep <- 50 # each Monte Carlo iteration gets its own folder folderNames <- file.path(outFolder, formatC(x = 1:nRep, width = 3, format = "d", flag = "0")) # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 1-node network where mosquitoes do not leave moveMat <- as.matrix(1) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) parameters # This time, lets put fitness cost on the homozygotes, giving the heterozygote # an advantage # These genotypes correspond to ones in the cube. Look at a base cube first, # then set this. # homozygotes are 60% as fit as heterozygote over their entire lifetime # Since omega is the adult daily death rate, we use the built-in function to # calculate our desired lifetime cost as applied daily dayOmega <- calcOmega(mu = bioParameters$muAd, lifespanReduction = 0.60) omegaNew <- c("AA"=dayOmega, "aa"=dayOmega) # setup cube cube <- cubeMendelian(omega = omegaNew) #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 100, occur every day, for 5 days. # There are 50 mosquitoes released every time. releasesParameters <- list(releasesStart=100, releasesNumber=5, releasesInterval=0, releaseProportion=50) # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run stochastic setupMGDrivE(stochasticityON = TRUE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 1 population. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, sampTime = 5, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=folderNames, verbose = FALSE) # run simulation MGDrivESim$multRun(verbose = FALSE) #################### # Post Analysis #################### # First, split output by patch # Second, aggregate females by their mate choice for(i in 1:nRep){ splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE) aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) } # plot output of first run to see effect plotMGDrivESingle(readDir = folderNames[1], totalPop = TRUE, lwd = 3.5, alpha = 1) # plot all 50 repetitions together plotMGDrivEMult(readDir = outFolder, lwd = 0.35, alpha = 0.75) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### outFolder <- "mgdrive" #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365 # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 2-node network where mosquitoes do not leave moveMat <- matrix(data = c(1,0,0,1), nrow = 2, ncol = 2) moveMat # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) parameters cube <- cubeMendelian() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur every day, for 5 days. # There are 50 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=5, releasesInterval=0, releaseProportion=50) # generate release vector releasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- releasesVector patchReleases[[1]]$femaleReleases <- releasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run deterministic setupMGDrivE(stochasticityON = FALSE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 2 populations. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=outFolder, verbose=FALSE) # run simulation MGDrivESim$oneRun(verbose = FALSE) #################### # Post Analysis #################### # split output by patch # Required for plotting later splitOutput(readDir = outFolder, verbose = FALSE, remFile = TRUE) # aggregate females by their mate choice # This reduces the female file to have the same columns as the male file aggregateFemales(readDir = outFolder, genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) # plot output to see effect plotMGDrivESingle(readDir = outFolder, lwd = 3.5, alpha = 1) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### outFolder <- "mgdrive" #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365 # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 2-node network with 1% per day migration rate moveMat <- matrix(data = c(0.99,0.01,0.01,0.99), nrow = 2, ncol = 2) moveMat # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) patchPops <- rep(adultPopEquilibrium,sitesNumber) #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) parameters cube <- cubeMendelian() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur every day, for 5 days. # There are 50 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=5, releasesInterval=0, releaseProportion=50) # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run deterministic setupMGDrivE(stochasticityON = FALSE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 2 populations. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=patchPops, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=outFolder, verbose = FALSE) # run simulation MGDrivESim$oneRun(verbose = FALSE) #################### # Post Analysis #################### # split output by patch # Required for plotting later splitOutput(readDir = outFolder, verbose = FALSE, remFile = TRUE) # aggregate females by their mate choice # This reduces the female file to have the same columns as the male file aggregateFemales(readDir = outFolder, genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) # plot output to see effect plotMGDrivESingle(readDir = outFolder, totalPop = TRUE, lwd = 3.5, alpha = 1) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### # This is slightly obtuse for vignette building reasons # Really, all you need is a base directory, then the repetitions in folders inside that. # Here, our base directory is "mgdrive", and the repetition folders are "001","002", etc. # So, the final structure is "mgdrive/001","mgdrive/002", etc. outFolder <- "mgdrive" dir.create(path = outFolder) #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365 # number of Monte Carlo iterations nRep <- 25 # each Monte Carlo iteration gets its own folder folderNames <- file.path(outFolder, formatC(x = 1:nRep, width = 3, format = "d", flag = "0")) # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 2-node network with 1% per day migration rate moveMat <- matrix(data = c(0.99,0.01,0.01,0.99), nrow = 2, ncol = 2) moveMat # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Mendelian cube with standard (default) parameters cube <- cubeMendelian() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur every day, for 5 days. # There are 50 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=5, releasesInterval=0, releaseProportion=50) # generate release vector releasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- releasesVector patchReleases[[1]]$femaleReleases <- releasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run stochastic setupMGDrivE(stochasticityON = TRUE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 2 populations. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=folderNames, verbose = FALSE) # run simulation MGDrivESim$multRun(verbose = FALSE) #################### # Post Analysis #################### # First, split output by patch # Second, aggregate females by their mate choice for(i in 1:nRep){ splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE) aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) } # plot output of first run to see effect plotMGDrivESingle(readDir = folderNames[1], totalPop = TRUE, lwd = 3.5, alpha = 1) # plot all 25 repetitions together plotMGDrivEMult(readDir = outFolder, lwd = 0.35, alpha = 0.75) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### # This is slightly obtuse for vignette building reasons # Really, all you need is a base directory, then the repetitions in folders inside that. # Here, our base directory is "mgdrive", and the repetition folders are "001","002", etc. # So, the final structure is "mgdrive/001","mgdrive/002", etc. outFolder <- "mgdrive" dir.create(path = outFolder) #################### # Simulation Parameters #################### # days to run the simulation tMax <- 365*2 # number of Monte Carlo iterations nRep <- 25 # each Monte Carlo iteration gets its own folder folderNames <- file.path(outFolder, formatC(x = 1:nRep, width = 3, format = "d", flag = "0")) # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 2-node network with 1% per day migration rate moveMat <- matrix(data = c(0.99,0.01,0.01,0.99), nrow = 2, ncol = 2) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Mendelian cube # This time, lets put fitness cost on the homozygotes, giving the heterozygote # an advantage # These genotypes correspond to ones in the cube. Look at a base cube first, # then set this. # Homozygotes are 60% as fit as heterozygote over their entire lifetime # Since omega is a daily death rate, we use the built-in function to calculate # our desired lifetime cost as applied daily dayOmega <- calcOmega(mu = bioParameters$muAd, lifespanReduction = 0.60) omegaNew <- c("AA"=dayOmega, "aa"=dayOmega) # setup cube cube <- cubeMendelian(omega = omegaNew) #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 100, occur every day, for 5 days. # There are 50 mosquitoes released every time. releasesParameters <- list(releasesStart=100, releasesNumber=5, releasesInterval=0, releaseProportion=50) # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run stochastic setupMGDrivE(stochasticityON = TRUE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 2 populations. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=folderNames, verbose = FALSE) # run simulation MGDrivESim$multRun(verbose = FALSE) #################### # Post Analysis #################### # First, split output by patch # Second, aggregate females by their mate choice for(i in 1:nRep){ splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE) aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) } # plot output of first run to see effect # per the structure above, we are reading "mgdrive/001" for the single plot plotMGDrivESingle(readDir = folderNames[1], totalPop = TRUE, lwd = 3.5, alpha = 1) # plot all 50 repetitions together # Here, we feed the function "mgdrive/", and it finds all repetition folders # inside that. plotMGDrivEMult(readDir = outFolder, lwd = 0.35, alpha = 0.75) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### outFolder <- "mgdrive" #################### # Simulation Parameters #################### # days to run the simulation, 2 years tMax <- 365*2 # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 1-node network where mosquitoes do not leave moveMat <- as.matrix(1) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Reciprocal translocation cube with standard (default) parameters cube <- cubeReciprocalTranslocations() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur once a week, for 5 weeks. # There are 100 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=5, releasesInterval=7, releaseProportion=100) # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run deterministic setupMGDrivE(stochasticityON = FALSE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 1 population. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=outFolder, verbose = FALSE) # run simulation MGDrivESim$oneRun(verbose = FALSE) #################### # Post Analysis #################### # split output by patch # Required for plotting later splitOutput(readDir = outFolder, remFile = TRUE, verbose = FALSE) # aggregate females by their mate choice # This reduces the female file to have the same columns as the male file aggregateFemales(readDir = outFolder, genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) # plot output to see effect plotMGDrivESingle(readDir = outFolder, totalPop = TRUE, lwd = 3.5, alpha = 1) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### outFolder <- "mgdrive" #################### # Simulation Parameters #################### # days to run the simulation, 2 years tMax <- 365*2 # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 1-node network where mosquitoes do not leave moveMat <- as.matrix(1) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) patchPops <- rep(adultPopEquilibrium,sitesNumber) #################### # Basic Inheritance pattern #################### # Reciprocal translocation cube with standard (default) parameters cube <- cubeReciprocalTranslocations() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur once a week, for 6 weeks. # There are 100 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=6, releasesInterval=7, releaseProportion=100) # generate release vector releasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- releasesVector patchReleases[[1]]$femaleReleases <- releasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run deterministic setupMGDrivE(stochasticityON = FALSE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 1 population. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=patchPops, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=outFolder, verbose = FALSE) # run simulation MGDrivESim$oneRun(verbose = FALSE) #################### # Post Analysis #################### # split output by patch # Required for plotting later splitOutput(readDir = outFolder, remFile = TRUE, verbose = FALSE) # aggregate females by their mate choice # This reduces the female file to have the same columns as the male file aggregateFemales(readDir = outFolder, genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) # plot output to see effect plotMGDrivESingle(readDir = outFolder, totalPop = TRUE, lwd = 3.5, alpha = 1) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls()) ## ----------------------------------------------------------------------------- #################### # Load libraries #################### library(MGDrivE) #################### # Output Folder #################### # This is slightly obtuse for vignette building reasons # Really, all you need is a base directory, then the repetitions in folders inside that. # Here, our base directory is "mgdrive", and the repetition folders are "001","002", etc. # So, the final structure is "mgdrive/001","mgdrive/002", etc. outFolder <- "mgdrive" dir.create(path = outFolder) #################### # Simulation Parameters #################### # days to run the simulation, 3 years tMax <- 365*3 # number of Monte Carlo iterations nRep <- 20 # each Monte Carlo iteration gets its own folder folderNames <- file.path(outFolder, formatC(x = 1:nRep, width = 3, format = "d", flag = "0")) # entomological parameters bioParameters <- list(betaK=20, tEgg=5, tLarva=6, tPupa=4, popGrowth=1.175, muAd=0.09) # a 1-node network where mosquitoes do not leave moveMat <- as.matrix(1) # parameters of the population equilibrium adultPopEquilibrium <- 500 sitesNumber <- nrow(moveMat) #################### # Basic Inheritance pattern #################### # Reciprocal translocation cube with standard (default) parameters cube <- cubeReciprocalTranslocations() #################### # Setup releases and batch migration #################### # set up the empty release vector # MGDrivE pulls things out by name patchReleases <- replicate(n=sitesNumber, expr={list(maleReleases=NULL,femaleReleases=NULL, eggReleases=NULL,matedFemaleReleases=NULL)}, simplify=FALSE) # choose release parameters # Releases start at time 25, occur once a week, for 6 weeks. # There are 100 mosquitoes released every time. releasesParameters <- list(releasesStart=25, releasesNumber=6, releasesInterval=7, releaseProportion=100) # generate male release vector maleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # generate female release vector femaleReleasesVector <- generateReleaseVector(driveCube=cube, releasesParameters=releasesParameters) # put releases into the proper place in the release list patchReleases[[1]]$maleReleases <- maleReleasesVector patchReleases[[1]]$femaleReleases <- femaleReleasesVector # batch migration is disabled by setting the probability to 0 # This is required because of the stochastic simulations, but doesn't make sense # in a deterministic simulation. batchMigration <- basicBatchMigration(batchProbs=0, sexProbs=c(.5,.5), numPatches=sitesNumber) #################### # Combine parameters and run! #################### # set MGDrivE to run stochastic setupMGDrivE(stochasticityON = TRUE, verbose = FALSE) # setup parameters for the network. This builds a list of parameters required for # every population in the network. In ths case, we have a network of 1 population. netPar <- parameterizeMGDrivE(runID=1, simTime=tMax, nPatch=sitesNumber, beta=bioParameters$betaK, muAd=bioParameters$muAd, popGrowth=bioParameters$popGrowth, tEgg=bioParameters$tEgg, tLarva=bioParameters$tLarva, tPupa=bioParameters$tPupa, AdPopEQ=adultPopEquilibrium, inheritanceCube = cube) # build network prior to run MGDrivESim <- Network$new(params=netPar, driveCube=cube, patchReleases=patchReleases, migrationMale=moveMat, migrationFemale=moveMat, migrationBatch=batchMigration, directory=folderNames, verbose = TRUE) # run simulation MGDrivESim$multRun(verbose = FALSE) #################### # Post Analysis #################### # First, split output by patch # Second, aggregate females by their mate choice for(i in 1:nRep){ splitOutput(readDir = folderNames[i], remFile = TRUE, verbose = FALSE) aggregateFemales(readDir = folderNames[i], genotypes = cube$genotypesID, remFile = TRUE, verbose = FALSE) } # plot output of first run to see effect plotMGDrivESingle(readDir = folderNames[1], lwd = 3.5, alpha = 1) # plot all 50 repetitions together plotMGDrivEMult(readDir = outFolder, lwd = 0.35, alpha = 0.75) ## ---- echo=FALSE-------------------------------------------------------------- #################### # Theory Comparison #################### # list male files mFiles <- list.files(path = outFolder, recursive = TRUE, pattern = "^M.*.csv$", full.names = TRUE) # read in simulation files successCount <- 0 for(f in mFiles){ # read in files, one by one hold <- matrix(data = scan(file = f, what = double(), sep = ",", skip = 1, quiet = TRUE), nrow = tMax, ncol = 1 + cube$genotypesN, byrow = TRUE) # check if the simulation was successful successCount <- successCount + (hold[tMax,10]!=0) } # percentage of success sCPerc <- successCount / nRep * 100 ## ---- echo=FALSE-------------------------------------------------------------- #################### # Cleanup before next run #################### unlink(x = outFolder, recursive = TRUE) rm(list=ls())