R/parameters.R

In pahanc/malariasimulation_Ace_params: An individual based model for malaria

#' @title Get model parameters
#' @description
#' get_parameters creates a named list of parameters for use in the model. These
#' parameters are passed to process functions. These parameters are explained in
#' "The US President's Malaria Initiative, Plasmodium falciparum transmission
#' and mortality: A modelling study."
#'
#' @param overrides a named list of parameter values to use instead of defaults
#' The parameters are defined below.
#'
#' fixed state transitions:
#'
#' * dd - the delay for humans to move from state D to A; default = 5
#' * dt - the delay for humans to move from state Tr to Ph; default = 5
#' * da - the delay for humans to move from state A to U; default = 200
#' * du - the delay for humans to move from state U to S; default = 110
#' * del - the delay for mosquitoes to move from state E to L; default = 6.64
#' * dl - the delay for mosquitoes to move from state L to P; default = 3.72
#' * dpl - the delay mosquitoes to move from state P to Sm; default = 0.643
#' * mup - the rate at which pupal mosquitoes die; default = 0.249
#' * mum - the rate at which developed mosquitoes die; default = 0.1253333
#'
#' immunity decay rates:
#'
#' * rm - decay rate for maternal immunity to clinical disease; default = 67.6952
#' * rvm - decay rate for maternal immunity to severe disease; default = 76.8365
#' * rb - decay rate for acquired pre-erythrocytic immunity; default = 3650
#' * rc - decay rate for acquired immunity to clinical disease; default = 10950
#' * rva - decay rate for acquired immunity to severe disease; default = 10950
#' * rid - decay rate for acquired immunity to detectability; default = 3650
#'
#' probability of pre-erythrocytic infection:
#'
#' * b0 - maximum probability due to no immunity; default = 0.59
#' * b1 - maximum reduction due to immunity; default = 0.5
#' * ib0 - scale parameter; default = 43.9
#' * kb - shape parameter; default = 2.16
#'
#' probability of clinical infection:
#'
#' * phi0 - maximum probability due to no immunity; default = 0.792
#' * phi1 - maximum reduction due to immunity; default = 0.00074
#' * ic0 - scale parameter; default = 18.02366
#' * kc - shape parameter; default = 2.36949
#'
#' probability of severe infection:
#'
#' * theta0 - maximum probability due to no immunity; default = 0.0749886
#' * theta1 - maximum reduction due to immunity; default = 0.0001191
#' * iv0 - scale parameter; default = 1.09629
#' * kv - shape parameter; default = 2.00048
#' * fv0 - age dependent modifier; default = 0.141195
#' * av - age dependent modifier; default = 2493.41
#' * gammav - age dependent modifier; default = 2.91282
#'
#' immunity reducing probability of detection:
#'
#' * fd0 - time-scale at which immunity changes with age; default = 0.007055
#' * ad - scale parameter relating age to immunity; default = 7993.5
#' * gammad - shape parameter relating age to immunity; default = 4.8183
#' * d1 - minimum probability due to immunity; default = 0.160527
#' * id0 - scale parameter; default = 1.577533
#' * kd - shape parameter; default = 0.476614
#'
#' immunity boost grace periods:
#'
#' * ub - period in which pre-erythrocytic immunity cannot be boosted; default = 7.2
#' * uc - period in which clinical immunity cannot be boosted; default = 6.06
#' * uv - period in which severe immunity cannot be boosted; default = 11.4321
#' * ud - period in which immunity to detectability cannot be boosted; default = 9.44512
#'
#' infectivity towards mosquitoes:
#'
#' * cd - infectivity of clinically diseased humans towards mosquitoes; default = 0.068
#' * gamma1 - parameter for infectivity of asymptomatic humans; default = 1.82425
#' * cu - infectivity of sub-patent infection; default = 0.0062
#' * ct - infectivity of treated infection; default = 0.021896
#'
#' unique biting rate:
#'
#' * a0 - age dependent biting parameter; default = 2920
#' * rho - age dependent biting parameter; default = 0.85
#' * sigma_squared - heterogeneity parameter; default = 1.67
#' * n_heterogeneity_groups - number discretised groups for heterogeneity, used
#' for sampling mothers; default = 5
#'
#' mortality parameters:
#'
#' * average_age - the average age of humans (in timesteps), this is only used 
#' if custom_demography is FALSE; default = 7663
#' * pcm - new-born clinical immunity relative to mother's; default = 0.774368
#' * pvm - new-born severe immunity relative to mother's; default = 0.195768
#' * me - early stage larval mortality rate; default = 0.0338
#' * ml - late stage larval mortality rate; default = 0.0348
#'
#' carrying capacity parameters:
#'
#' * model_seasonality - boolean switch TRUE iff the simulation models seasonal rainfall; default = FALSE 
#' * g0 - rainfall fourier parameter; default = 2
#' * g - rainfall fourier parameter; default = 0.3, 0.6, 0.9
#' * h - rainfall fourier parameters; default = 0.1, 0.4, 0.7
#' * gamma - effect of density dependence on late instars relative to early
#' instars; default = 13.25
#' * rainfall_floor - the minimum rainfall value (must be above 0); default 0.001
#'
#' initial state proportions:
#'
#' * s_proportion - the proportion of `human_population` that begin as susceptible; default = 0.420433246
#' * d_proportion - the proportion of `human_population` that begin with
#' clinical disease; default = 0.007215064
#' * a_proportion - the proportion of `human_population` that begin as
#' asymptomatic; default = 0.439323667
#' * u_proportion - the proportion of `human_population` that begin as
#' subpatents; default = 0.133028023
#' * t_proportion - the proportion of `human_population` that begin treated; default = 0
#'
#' initial immunity values:
#'
#' * init_icm - the immunity from clinical disease at birth; default = 0
#' * init_ivm - the immunity from severe disease at birth; default = 0
#' * init_ib  - the initial pre-erythrocitic immunity; default = 0
#' * init_ica - the initial acquired immunity from clinical disease; default = 0
#' * init_iva - the initial acquired immunity from severe disease; default = 0
#' * init_id  - the initial acquired immunity to detectability; default = 0
#'
#' incubation periods:
#'
#' * de - Duration of the human latent period of infection; default = 12
#' * delay_gam - Lag from parasites to infectious gametocytes; default = 12.5
#' * dem - Extrinsic incubation period in mosquito population model; default = 10
#'
#' vector biology:
#' species specific values are vectors
#'
#' * beta - the average number of eggs laid per female mosquito per day; default = 21.2
#' * total_M - the initial number of adult mosquitos in the simulation; default = 1000
#' * init_foim - the FOIM used to calculate the equilibrium state for mosquitoes; default = 0
#' * species - names of the species in the simulation; default = "All" 
#' * species_proportions - the relative proportions of each species; default = 1
#' * blood_meal_rates - the blood meal rates for each species; default = 0.3333333333
#' * Q0 - proportion of blood meals taken on humans; default = 0.92
#' * foraging_time - time spent taking blood meals; default = 0.69
#'
#' feeding cycle: 
#' please set vector control strategies using `set_betnets` and `set_spraying`
#'
#' * bednets - boolean for if bednets are enabled; default = FALSE 
#' * phi_bednets - proportion of bites taken in bed; default = 0.85
#' * k0 - proportion of females bloodfed with no net; default = 0.699
#' * spraying - boolean for if indoor spraying is enabled; default = FALSE 
#' * phi_indoors - proportion of bites taken indoors; default = 0.90
#'
#' treatment parameters:
#' please set treatment parameters with the convenience functions in
#' `drug_parameters.R`
#'
#' * drug_efficacy - a vector of efficacies for available drugs; default = turned off
#' * drug_rel_c - a vector of relative onward infectiousness values for drugs; default = turned off
#' * drug_prophylaxis_shape - a vector of shape parameters for weibull curves to
#' model prophylaxis for each drug; default = turned off
#' * drug_prophylaxis_scale - a vector of scale parameters for weibull curves to
#' model prophylaxis for each drug; default = turned off
#' * clinical_treatment_drugs - a vector of drugs that are available for
#' clinically diseased (these values refer to the index in drug_* parameters); default = NULL, NULL, NULL
#' * clinical_treatment_coverage - a vector of coverage values for each drug; default = NULL, NULL, NULL
#'
#' RTS,S paramters: 
#' please set RTS,S strategies with the convenience functions in
#' `vaccine_parameters.R:set_rtss_epi`
#' `vaccine_parameters.R:set_mass_rtss`
#'
#' * rtss_doses - the dosing schedule before the vaccine takes effect; default =
#' c(0, 1.5 * 30, 3 * 30)
#' default = 365
#' * rtss_vmax - the maximum efficacy of the vaccine; default = 0.93
#' * rtss_alpha - shape parameter for the vaccine efficacy model; default = 0.74
#' * rtss_beta - scale parameter for the vaccine efficacy model; default = 99.4
#' * rtss_cs - peak parameters for the antibody model (mean and std. dev); default = 6.37008, 0.35
#' * rtss_cs_boost - peak parameters for the antibody model for booster rounds (mean and std. dev); default = 5.56277, 0.35
#' * rtss_rho - delay parameters for the antibody model (mean and std. dev); default = 2.37832, 1.00813
#' * rtss_rho_boost - delay parameters for the antibody model for booster rounds (mean and std. dev); default = 1.03431, 1.02735
#' * rtss_ds - delay parameters for the antibody model, short-term weaning (mean and std. dev); default = 3.74502, 0.341185 (White MT et al. 2015 Lancet ID)
#' * rtss_dl - delay parameters for the antibody model, long-term weaning (mean and std. dev); default = 6.30365, 0.396515 (White MT et al. 2015 Lancet ID)
#'
#' MDA and SMC parameters:
#' please set these parameters with the convenience functions in `mda_parameters.R`
#'
#' TBV parameters: 
#' please set TBV parameters with the convenience functions in
#' `vaccine_parameters.R:set_tbv`
#'
#' * tbv_mt - effect on treated infectiousness; default = 35
#' * tbv_md - effect on diseased infectiousness; default = 46.7
#' * tbv_ma - effect on asymptomatic infectiousness; default = 3.6
#' * tbv_mu - effect on subpatent infectiousness; default = 0.8
#' * tbv_k  - scale parameter for effect on infectiousness; default = 0.9
#' * tbv_tau - peak antibody parameter; default = 22
#' * tbv_rho - antibody component parameter; default = 0.7
#' * tbv_ds - antibody short-term delay parameter; default = 45
#' * tbv_dl - antibody long-term delay parameter; default = 591
#' * tbv_tra_mu - transmission reduction parameter; default = 12.63
#' * tbv_gamma1 - transmission reduction parameter; default = 2.5
#' * tbv_gamma2 - transmission reduction parameter; default = 0.06
#'
#' rendering:
#' All values are in timesteps and all ranges are inclusive
#'
#' * prevalence_rendering_min_ages - the minimum ages for clinical prevalence
#' outputs; default = 730
#' * prevalence_rendering_max_ages - the corresponding max ages; default = 3650
#' * incidence_rendering_min_ages - the minimum ages for incidence
#' outputs (includes asymptomatic microscopy +); default = turned off
#' * incidence_rendering_max_ages - the corresponding max ages; default = turned off 
#' * clinical_incidence_rendering_min_ages - the minimum ages for clinical incidence outputs (symptomatic); default = 0
#' * clinical_incidence_rendering_max_ages - the corresponding max ages; default = 1825
#' * severe_incidence_rendering_min_ages - the minimum ages for severe incidence
#' outputs; default = turned off
#' * severe_incidence_rendering_max_ages - the corresponding max ages; default = turned off
#'
#' miscellaneous:
#'
#' * human_population - the initial number of humans to model; default = 100
#' * human_population_timesteps - the timesteps at which the population should
#' change; default = 0
#' * mosquito_limit - the maximum number of mosquitoes to allow for in the
#' simulation; default = 1.00E+05
#' * individual_mosquitoes - boolean whether adult mosquitoes are modeled
#' individually or compartmentally; default = TRUE 
#' * r_tol - the relative tolerance for the ode solver; default = 1e-4
#' * a_tol - the absolute tolerance for the ode solver; default = 1e-4
#' * ode_max_steps - the max number of steps for the solver; default = 1e6
#' * enable_heterogeneity - boolean whether to include heterogeneity in biting
#' rates; default = TRUE 
#'
#' @export
get_parameters <- function(overrides = list(),square_number=square_number,run_number=run_number) {
  
  supp_filename_gamb<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Suppress/mosq_supp_gamb",square_number,"_",run_number,".csv",sep="")
  mosq_suppression_gamb<-unlist(read.csv(supp_filename_gamb,header=F,colClasses="numeric"))
  #mosq_suppression<-unlist(read.csv("/Imperial March 2021/SEA project/mosq_suppression.csv",header=F,colClasses="numeric"))
  dimnames(mosq_suppression_gamb)<-NULL
  mosq_suppression_gamb<-as.vector(mosq_suppression_gamb)
  ##No Drive:
  #mosq_suppression_gamb<-rep(1,length(mosq_suppression_gamb))

  supp_filename_arab<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Suppress/mosq_supp_arab",square_number,"_",run_number,".csv",sep="")
  mosq_suppression_arab<-unlist(read.csv(supp_filename_arab,header=F,colClasses="numeric"))
  dimnames(mosq_suppression_arab)<-NULL
  mosq_suppression_arab<-as.vector(mosq_suppression_arab)
  # #No Drive:
  #mosq_suppression_arab<-rep(1,length(mosq_suppression_arab))


  #supp_filename_fun<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/mosq_supp_het_stronger_fun",square_number,"_",run_number,".csv",sep="")
  #mosq_suppression_fun<-unlist(read.csv(supp_filename_fun,header=F,colClasses="numeric"))
  #dimnames(mosq_suppression_fun)<-NULL
  #mosq_suppression_fun<-as.vector(mosq_suppression_fun)
  # #No Drive:
  mosq_suppression_fun<-rep(1,length(mosq_suppression_gamb))#NO FUNESTUS IN SOUROU
  
  supp_filename_new<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Suppress/mosq_supp_arab",square_number,"_",run_number,".csv",sep="")
  mosq_suppression_new<-unlist(read.csv(supp_filename_arab,header=F,colClasses="numeric"))
  dimnames(mosq_suppression_new)<-NULL
  mosq_suppression_new<-as.vector(mosq_suppression_new)
  # #No Drive:
  mosq_suppression_new<-rep(1,length(mosq_suppression_new))

  mosq_supp_lst<-list()
  mosq_supp_lst[[1]]<-mosq_suppression_gamb
  mosq_supp_lst[[2]]<-mosq_suppression_arab
  mosq_supp_lst[[3]]<-mosq_suppression_fun
  mosq_supp_lst[[4]]<-mosq_suppression_new
  

  seas_filename_gamb<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Emerge/mosq_seasonality_gamb",square_number,"_",run_number,".csv",sep = "")
  mosq_seasonality_gamb<-unlist(read.csv(seas_filename_gamb,header=F,colClasses="numeric"))
  dimnames(mosq_seasonality_gamb)<-NULL
  mosq_seasonality_gamb<-as.vector(mosq_seasonality_gamb)

  seas_filename_arab<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Emerge/mosq_seasonality_arab",square_number,"_",run_number,".csv",sep = "")
  mosq_seasonality_arab<-unlist(read.csv(seas_filename_arab,header=F,colClasses="numeric"))
  dimnames(mosq_seasonality_arab)<-NULL
  mosq_seasonality_arab<-as.vector(mosq_seasonality_arab)

  #seas_filename_fun<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Emerge/mosq_seasonality_fun",square_number,"_",run_number,".csv",sep = "")
  #mosq_seasonality_fun<-unlist(read.csv(seas_filename_fun,header=F,colClasses="numeric"))
  #dimnames(mosq_seasonality_fun)<-NULL
  mosq_seasonality_fun<-as.vector(mosq_seasonality_gamb)#NO FUNESTUS IN SOUROU
  
  seas_filename_new<-paste("Q:\\for_hpc\\Seas and supp species specific/With funestus suppression/Mar 24 Sourou/Mar 24 Sourou Emerge/mosq_seasonality_arab",square_number,"_",run_number,".csv",sep = "")
  mosq_seasonality_new<-unlist(read.csv(seas_filename_new,header=F,colClasses="numeric"))
  dimnames(mosq_seasonality_new)<-NULL
  mosq_seasonality_new<-as.vector(mosq_seasonality_new)

  mosq_seas_lst<-list()
  mosq_seas_lst[[1]]<-mosq_seasonality_gamb
  mosq_seas_lst[[2]]<-mosq_seasonality_arab
  mosq_seas_lst[[3]]<-mosq_seasonality_fun
  mosq_seas_lst[[4]]<-mosq_seasonality_new


  parameters <- list(
    use_Ace_mosq = FALSE,
    mosq_suppression = mosq_supp_lst,
    mosq_seasonality = mosq_seas_lst,
    emergence = 0,
    dens_indep = FALSE,
    dd    = 5,
    dt    = 5,
    da    = 200,
    du    = 110,
    del   = 6.64,
    dl    = 3.72,
    dpl   = .643,
    mup   = .249,
    mum   = .1253333,
    mum_atsb = 0.09,
    sigma_squared   = 1.67,
    n_heterogeneity_groups = 5,
    # immunity decay rates
    rm    = 67.6952,
    rvm   = 76.8365,
    rb    = 10 * 365,
    rc    = 30 * 365,
    rva   = 30 * 365,
    rid   = 10 * 365,
    # blood immunity parameters
    b0    = 0.59,
    b1    = 0.5,
    ib0   = 43.9,
    kb    = 2.16,
    # immunity boost grace periods
    ub    = 7.2,
    uc    = 6.06,
    uv    = 11.4321,
    ud    = 9.44512,
    # infectivity towards mosquitos
    cd    = 0.068,
    gamma1= 1.82425,
    cu    = 0.0062,
    ct    = 0.021896,
    # unique biting rate
    a0    = 8 * 365,
    rho   = .85,
    # clinical immunity parameters
    phi0 = .792,
    phi1 = .00074,
    ic0   = 18.02366,
    kc    = 2.36949,
    # severe disease immunity parameters
    theta0  = .0749886,
    theta1  = .0001191,
    kv      = 2.00048,
    fv0     = 0.141195,
    av      = 2493.41,
    gammav  = 2.91282,
    iv0     = 1.09629,
    # delay for infection
    de        = 12,
    delay_gam = 12.5,
    dem       = 10,
    # asymptomatic immunity parameters
    fd0   = 0.007055,
    ad    = 21.9 * 365,
    gammad= 4.8183,
    d1    = 0.160527,
    id0   = 1.577533,
    kd    = .476614,
    # mortality parameters
    average_age = 8030,
    pcm   = .774368,
    pvm   = .195768,
    # carrying capacity parameters
    g0    = 2,
    g     = c(.3, .6, .9),
    h     = c(.1, .4, .7),
    gamma = 13.25,
    model_seasonality = FALSE,
    rainfall_floor = 0.001,
    # larval mortality rates
    me    = .0338,
    ml    = .0348,
    # initial state proportions
    s_proportion = 0.420433246,
    d_proportion = 0.007215064,
    a_proportion = 0.439323667,
    u_proportion = 0.133028023,
    t_proportion = 0,
    # initial immunities
    init_ica = 0,
    init_iva = 0,
    init_icm = 0,
    init_ivm = 0,
    init_id  = 0,
    init_ib  = 0,
    # vector biology
    beta     = 21.2,
    total_M  = 1000,
    init_foim= 0,
    # order of species: An gambiae s.s, An arabiensis, An funestus
    species             = 'All',
    species_proportions = 1,
    blood_meal_rates    = 1/3,
    Q0                  = .92,
    foraging_time       = .69,
    # atsb
    atsb = FALSE,
    # bed nets
    bednets = FALSE,
    phi_bednets = .85,
    k0 = .699,
    # indoor spraying
    spraying = FALSE,
    phi_indoors = .90,
    # treatment
    drug_efficacy          = numeric(0),
    drug_rel_c             = numeric(0),
    drug_prophylaxis_shape = numeric(0),
    drug_prophylaxis_scale = numeric(0),
    clinical_treatment_drugs     = list(),
    clinical_treatment_timesteps = list(),
    clinical_treatment_coverages = list(),
    # rts,s
    rtss           = FALSE,
    rtss_doses     = c(0, 1.5 * 30, 3 * 30),
    rtss_vmax      = .93,
    rtss_alpha     = .74,
    rtss_beta      = 99.4,
    rtss_cs        = c(6.37008, 0.35),
    rtss_cs_boost  = c(5.56277, 0.35),
    rtss_rho       = c(2.37832, 1.00813),
    rtss_rho_boost = c(1.03431, 1.02735),
    rtss_ds        = c(3.74502, 0.341185),
    rtss_dl        = c(6.30365, 0.396515),
    # MDA
    mda = FALSE,
    mda_drug = 0,
    mda_timesteps = NULL,
    mda_coverages = NULL,
    mda_min_age = -1,
    mda_max_age = -1,
    smc = FALSE,
    smc_drug = 0,
    smc_timesteps = NULL,
    smc_coverages = NULL,
    smc_min_age = -1,
    smc_max_age = -1,
    # tbv
    tbv = FALSE,
    tbv_mt = 35,
    tbv_md = 46.7,
    tbv_ma = 3.6,
    tbv_mu = 0.8,
    tbv_k = 0.9,
    tbv_tau = 22,
    tbv_rho = .7,
    tbv_ds = 45,
    tbv_dl = 591,
    tbv_tra_mu = 12.63,
    tbv_gamma1 = 2.5,
    tbv_gamma2 = .06,
    tbv_timesteps = NULL,
    tbv_coverages = NULL,
    tbv_ages = NULL,
    # rendering
    prevalence_rendering_min_ages = 2 * 365,
    prevalence_rendering_max_ages = 10 * 365,
    incidence_rendering_min_ages = numeric(0),
    incidence_rendering_max_ages = numeric(0),
    clinical_incidence_rendering_min_ages = numeric(0),
    clinical_incidence_rendering_max_ages = 5 * 365,
    severe_prevalence_rendering_min_ages = numeric(0),
    severe_prevalence_rendering_max_ages = numeric(0),
    severe_incidence_rendering_min_ages = numeric(0),
    severe_incidence_rendering_max_ages = numeric(0),
    # misc
    custom_demography = FALSE,
    human_population = 100,
    human_population_timesteps = 0,
    mosquito_limit   = 100 * 1000,
    individual_mosquitoes = TRUE,
    enable_heterogeneity = TRUE,
    r_tol = 1e-4,
    a_tol = 1e-4,
    ode_max_steps = 1e6
  )

 #specify a fixed total_M in the absense of interventions
  #parameters$total_M_orig<-parameters$total_M*mosq_seasonality_gamb[1:365]

  # Override parameters with any client specified ones
  if (!is.list(overrides)) {
    stop('overrides must be a list')
  }

  for (name in names(overrides)) {
    if (!(name %in% names(parameters))) {
      stop(paste('unknown parameter', name, sep=' '))
    }
    parameters[[name]] <- overrides[[name]]
  }

  props <- c(
    parameters$s_proportion,
    parameters$d_proportion,
    parameters$a_proportion,
    parameters$u_proportion,
    parameters$t_proportion
  )

  if (!approx_sum(props, 1)) {
    stop("Starting proportions do not sum to 1")
  }

  parameters
}

#' @title Parameterise total_M and carrying capacity for mosquitos from EIR
#'
#' @description NOTE: the inital EIR is likely to change unless the rest of the
#' model is in equilibrium. NOTE: please set seasonality first, since the mosquito_limit
#' will estimate an upper bound from the peak season.
#'
#' max_total_M is calculated using the equilibrium solution from "Modelling the
#' impact of vector control interventions on Anopheles gambiae population
#' dynamics"
#'
#' @param parameters the parameters to modify
#' @param EIR to work from
#' @export
parameterise_mosquito_equilibrium <- function(parameters, EIR) {
  parameterise_total_M(parameters, equilibrium_total_M(parameters, EIR))
}

#' @title Parameterise total_M
#'
#' @description Sets total_M and an upper bound for the number of mosquitoes in
#' the simulation. NOTE: please set seasonality first, since the mosquito_limit
#' will estimate an upper bound from the peak season.
#'
#' @param parameters the parameters to modify
#' @param total_M the initial adult mosquitoes in the simulation
#' @export
parameterise_total_M <- function(parameters, total_M) {
  parameters$total_M <- total_M
  if (!parameters$individual_mosquitoes) {
    return(parameters)
  }
  max_total_M <- 0
  for (i in seq_along(parameters$species)) {
    species_M <- total_M * parameters$species_proportions[[i]]
    K0 <- calculate_carrying_capacity(parameters, species_M, i)
    R_bar <- calculate_R_bar(parameters)
    max_K <- max(vnapply(seq(365), function(t) {
      carrying_capacity(
        t,
        parameters$model_seasonality,
        parameters$g0,
        parameters$g,
        parameters$h,
        K0,
        R_bar,
        parameters$rainfall_floor
      )
    }))
    omega <- calculate_omega(parameters, i)
    mum <- weighted.mean(parameters$mum, parameters$species_proportions)
    max_total_M <- max_total_M + max_K * (
      1 / (
        2 * parameters$dl * mum * (
          1 + parameters$dpl * parameters$mup
        )
      )
    ) * (
      1 / (
        parameters$gamma * (omega + 1)
      )
    ) * (
      omega / (parameters$ml * parameters$del) - (
        1 / (parameters$ml * parameters$dl)
      ) - 1
    )
  }
  parameters$mosquito_limit <- ceiling(max_total_M * 5) #Allow for random fluctuations
  parameters
}