Given prevalence of disease in humans (modeled as an SEIR: Susceptible-Latent-Infected-Recovered process with birth and death) and entomological parameters of transmission, this function calculates the quasi-stationary distribution of adult female mosquitoes across SEI (Susceptible-Exposed-Infectious) stages, allowing for Erlang distributed E stage.

  node_list = "b",
  NF = NULL,
  phi = 0.5,
  NH = NULL,
  log_dd = TRUE,
  pop_ratio_Aq = NULL,
  pop_ratio_F = NULL,
  pop_ratio_M = NULL,
  pop_ratio_H = c(1, 0, 0, 0),



a named list of parameters (see details)


a character vector specifying what type of nodes to create; (m = a node with only mosquitoes, h = a node with only humans, b = a node with both humans and mosquitoes)


vector of female mosquitoes at equilibrium, for mosquito-only nodes


sex ratio of mosquitoes at emergence


vector of humans at equilibrium, for human-only nodes


Boolean: TRUE implies logistic density dependence, FALSE implies Lotka-Volterra model


the set of places (P) (see details)


May be empty; if not, a named vector or matrix. (see details)


May be empty; if not, a named vector or matrix. (see details)


May be empty; if not, a named vector or matrix. (see details)


Prevalence in human-only nodes, default is all susceptible


an inheritance cube from the MGDrivE package (e.g. cubeMendelian)


a vector of the equilibrium number of females in each SEI stage


This function handles 3 types of nodes: Human only, mosquito only, and nodes with both. These nodes are set using the node_list parameter. Mosquito-only node equilibrium calls equilibrium_lifeycle, which follows one of two models: classic logistic dynamics or the Lotka-Volterra competition model. This is determined by the parameter log_dd, and it changes elements of the return list: K is returned for logistic dynamics, or gamma is returned for Lotka-Volterra dynamics. This is parameterized with the NF parameter to define the adult female numbers. This parameter only needs to be supplied if there are mosquito-only nodes.

Human-only nodes don't require any equilibrium calculations. These nodes use the NH and pop_ratio_H to set adult human populations and infection rates in nodes. These two parameters only need to be supplied if there are human-only nodes. pop_ratio_H needs to be a matrix with the number of rows equal to the number of human-only patches, and 4 columns. The columns are assumed to be fractions of the population in "S", "E", "I", or "R" states, and every row must sum to 1.

For human and mosquito nodes, this function calls make_Q_SEI to construct the infinitesimal generator matrix which is used to solve for the quasi-stationary (stochastic) or equilibrium (deterministic) distribution of mosquitoes over stages. Parameters are provided by params.

For information on the method used to solve this distribution, see section "3.1.3 Nonsingularity of the Subintensity Matrix" of:

  • Bladt, Mogens, and Bo Friis Nielsen. Matrix-exponential distributions in applied probability. Vol. 81. New York: Springer, 2017.

The places (spn_P) object is generated from one of the following: spn_P_lifecycle_node, spn_P_lifecycle_network, spn_P_epiSIS_node, spn_P_epiSIS_network, spn_P_epiSEIR_node, or spn_P_epiSEIR_network.

The initial population genotype ratios are set by supplying the pop_ratio_Aq, pop_ratio_F, and pop_ratio_M values. The default value is NULL, and the function will use the wild-type alleles provided in the cube object. However, one can supply several different objects to set the initial genotype ratios. All genotypes provided must exist in the cube (this is checked by the function). If a single, named vector is provided, then all patches will be initialized with the same ratios. If a matrix is provided, with the number of columns (and column names) giving the initial genotypes, and a row for each patch, each patch can be set to a different initial ratio. The three parameters do not need to match each other.

The params argument supplies all of the ecological and epidemiological parameters necessary to calculate equilibrium values. This is used to set the initial population distribution and during the simulation to maintain equilibrium. This params must include the following named parameters, noted as being the same as lifecycle parameters, or new for the epidemiological equilibrium

  • (Lifecycle parameters)

    • qE: inverse of mean duration of egg stage

    • nE: shape parameter of Erlang-distributed egg stage

    • qL: inverse of mean duration of larval stage

    • nL: shape parameter of Erlang-distributed larval stage

    • qP: inverse of mean duration of pupal stage

    • nP: shape parameter of Erlang-distributed pupal stage

    • muE: egg mortality

    • muL: density-independent larvae mortality

    • muP: pupae mortality

    • muF: adult female mortality

    • muM: adult male mortality

    • beta: egg-laying rate, daily

    • nu: mating rate of unmated females

  • (Epidemiological parameters)

    • NH: number of humans, can be a vector

    • X: SEIR prevalence in humans, can be a vector of length 4 for 1 node, or a matrix for many nodes

    • NFX: number of female mosquitoes, only required if any prevalence (X) is zero

    • b: mosquito to human transmission efficiency, can be a vector

    • c: human to mosquito transmission efficiency, can be a vector

    • r: rate of recovery in humans (1/duration of infectiousness)

    • muH: death rate of humans (1/avg lifespan)

    • f: rate of blood feeding

    • Q: human blood index

    • qEIP: related to scale parameter of Gamma distributed EIP (1/qEIP is mean length of EIP)

    • nEIP: shape parameter of Gamma distributed EIP

    • delta: inverse duration of the latent stage (E)

The return list contains all of the parameters necessary later in the simulations.

For equilibrium without epidemiological parameters, see equilibrium_lifeycle. For equilibrium without latent humans (SIS dynamics), see equilibrium_SEI_SIS.