Developing L Modules

Models for Aquatic Mosquito Ecology

Source: vignettes/L-How_To_Add_Model.Rmd

This vignette discusses how to add a new L module to ramp.library. As background, we strongly recommend reading the ramp.xds vignette about the L Interface.

Each L module implements a unique model family describing human/host infection dynamics and immunity that can be used by ramp.xds. Before writing a new module, please check and see if the model family is already part of the SimBA library.

If you need help, please contact Professor David L Smith

We have written this vignette to produce a new template file using knitr::purl. The template can be modified to build a new L module.

Do this:

  1. Download the markdown file that created this document “X-How_To_Add_Lodel.Rmd.”

  2. Set your working directory to the location of the downloaded copy of “X-How_To_Add_Lodel.Rmd”

  3. Execute these commands:

library(knitr)
knitr::purl("L-How_To_Add_Model.Rmd", 
            output = "Lmodule-newLname.R", 
            documentation = 0, 
            quiet=TRUE)

  1. You should see a new file called Lmodule-newLname.R in the working directory. You can now move and rename it.

  2. Open the file and do a global search for newLname and replace it with the a short string that will define the model as its Xname.

  3. Edit the functions following the directions below.

Style Guide

We recommend that new models conform to a set of style guidelines:

  • All the names of the variables by name using the indices attached to L_obj

  • All documentation is written in roxygen2

  • Functions should use with(*,{...}) to make the code easier to read.

  • A vignette should accompany each model, including citations to the relevant publications, and proposed tests of the code.

  • The functions should easy to relate to the equations defining a model, following notation defined in an accompanying vignette.

Dynamics

dLdt

The sections that need to be modified are marked # Change this

  • Replace L with the your first variable (e.g. mosquito population density), and replace F_L with the formula that computes the derivatives. Do the same with Y and add lines for each variable.

  • Lodify derivs = c(dL, dY, ...) to return the LHS of the derivative equations you just wrote.

#' @title Compute the derivatives for parasite infection dynamics in human population strata 
#' @description Implements [dLdt] for the newLname model
#' @inheritParams ramp.xds::dLdt
#' @return a [numeric] vector
#' @export
dLdt.newLname <- function(t, y, xds_obj, i) {

  # do not change this
  eta <- xds_obj$terms$eta[[i]]

  # do not change this
  # attach the variables by name 
  with(get_L_vars(y, xds_obj, i),{
    # do not change this
    # attach the parameters by name 
    with(xds_obj$L_obj[[i]], {
      # Do not change this 
      dL <- F_L(eta, L) 
      ... 
      
      # Change this 
      derivs = c(dL, ...) 
      
      # return the derivatives 
      return(derivs)
    })
  })
}

Update_Lt 

#' @title Compute the derivatives for parasite infection dynamics in human population strata 
#' @description Implements [UpdateLt] for the newLname model
#' @inheritParams ramp.xds::UpdateLt
#' @return a [numeric] vector
#' @export
Update_Lt.newLname<- function(t, y, xds_obj, s) {
  
  # do not change this
  eta <- xds_obj$terms$eta[[i]]

   # do not change this
  # attach the variables by name 
  with(get_L_vars(y, xds_obj, i),{
    # do not change this
    # attach the parameters by name 
    with(xds_obj$L_obj[[i]], {
      # Change this
      Lt <- L + Lambda + F_L(L) 
      ... 
      
      # Change this 
      states = c(Lt, ...) 
      
      # return the derivatives 
      return(states)
    })
  })
}

L Lodel Object

make_L_obj

  • In the following, replace p1 with the first parameter, p2 with the second, and so on.

  • Unless you know what you’re doing, leave the “Ports for human / host demography” unchanged.

  • If you don’t want to add ports for mass treatment, delete the line with mda and msat and leave mda=FALSE and msat=FALSE in skill_set_L.newLname. If you do want your module to have that functionality, add the effects of mass treatment to the parameters that define state transitions in dLdt or Update_Lt and set mda=TRUE and msat=TRUE in the skill set.

#' @title Lake parameters for newLname human model, with defaults
#' @param nPatches is the number of population strata
#' @param options a [list] that could overwrite defaults
#' @param p1 the first parameter 
#' @param p2 the second parameter 
#' @param p3 the third parameter 
#' @return a [list]
#' @export
make_L_obj_newLname = function(nPatches, options=list(),
                         p1=1, p2=2, p3=3){
  with(options,{
    L_obj = list()
    class(L_obj) <- c("newLname")

    # Change this
    L_obj$p1 = checkIt(p1, nPatches)
    L_obj$p2 = checkIt(p2, nPatches)
    L_obj$p3 = checkIt(p3, nPatches)

    # Don't change this
     L_obj <- setup_K_obj(nPatches, L_obj) 
    
    return(L_obj)
  })}

setup_L_obj

If you did a global search and replace on newLname, this function will not need to be modified

#' @title Setup L_obj.newLname
#' @description Implements [setup_L_obj] for the newLname model
#' @inheritParams ramp.xds::setup_L_obj
#' @return a [list] vector
#' @export
setup_L_obj.newLname = function(newLname, xds_obj, i, options=list()){
  xds_obj$L_obj[[i]] = make_L_obj_newLname(xds_obj$nPatches, options)
  return(xds_obj)
}

Parameters

Variables

setup_L_ix

In the following, each variable in your model gets its own pair of lines. Replace L in L_ix with the first variable name, Y in Y_ix with the second, and add lines for each variable by copying and pasting and editing the pair of lines. Be sure the replacement gets applied to all four instances.

#' @title Add indices for human population to parameter list
#' @description Implements [setup_L_ix] for the newLname model.
#' @inheritParams ramp.xds::setup_L_ix
#' @return none
#' @importFrom utils tail
#' @export
setup_L_ix.newLname <- function(xds_obj, i) {with(xds_obj,{
  
  L_ix <- seq(from = max_ix+1, length.out=xds_obj$nPatches[i])
  max_ix <- tail(L_ix, 1)

  Y_ix <- seq(from = max_ix+1, length.out=xds_obj$nPatches[i])
  max_ix <- tail(Y_ix, 1)

  ... 
  
  xds_obj$max_ix = max_ix
  xds_obj$L_obj[[i]]$ix = list(L_ix=L_ix, Y_ix=Y_ix, ...)
  return(xds_obj)
})}

get_L_vars

In the following, each variable in your model gets its own pair of lines.

  • Leave H as-is

  • Replace L in L and L_ix with the first variable name.

  • Replace Y in Y and Y_ix with the second variable name.

  • Copy the form for every other variable

  • The variable X0 is the class that gets newborns. If you use it or if you want to have it available, compute it and add it to the returned list.

#' @title Return the variables as a list
#' @description This method dispatches on the type of `xds_obj$L_obj`
#' @inheritParams ramp.xds::get_L_vars
#' @return a [list]
#' @export
get_L_vars.newLname <- function(y, xds_obj, i) {
  with(xds_obj$L_obj[[i]]$ix,{
      H = y[H_ix]
      L = y[L_ix]
      Y = y[Y_ix]
      ...
      X0 = H-L-Y 
      return(list(H=H, L=L, Y=Y,...,X0=X0))})
}

parse_L_orbits

In the following, each variable in your model gets its own pair of lines.

  • Leave H as-is

  • Replace L in L and L_ix with the first variable name.

  • Replace Y in Y and Y_ix with the second variable name.

  • Copy the form for every other variable

  • The variable X0 is the class that gets newborns. If you use it or if you want to have it available, compute it and add it to the returned list.

#' @title parse the output of deSolve and return variables for the newLname model
#' @description Implements [parse_L_orbits] for the newLname model
#' @inheritParams ramp.xds::parse_L_orbits
#' @return a named list  
#' @export
parse_L_orbits.newLname <- function(outputs, xds_obj, i) {
  with(xds_obj$L_obj[[i]]$ix,{
    H <- outputs[,H_ix]
    L <- outputs[,L_ix]
    Y <- outputs[,Y_ix]
    X0 <- H-L-Y
    vars <- list(H=H, L=L, Y=Y, X0=X0)
    return(vars)
})}

Initial Values

setup_L_inits

If you did the global search and replace, you won’t need to change this function.

#' @title Setup initial values for *newLname*
#' 
#' @inheritParams ramp.xds::setup_L_inits
#' 
#' @return a **`ramp.xds`** object 
#' 
#' @export
setup_L_inits.newLname = function(xds_obj, H, i, options=list()){
  xds_obj$L_obj[[i]]$inits = make_L_inits_newLname(xds_obj$nPatches[i], H, options)
  return(xds_obj)
}

make_L_inits 

#' @title Lake initial values for the Xname human model, with defaults
#' @param nPatches the number of strata in the model
#' @param options a [list] to overwrite defaults
#' @param L0 the initial value for L 
#' @param Y0 the initial value for L 
#' @return a [list]
#' @export
make_L_inits_newLname = function(nPatches, H, options = list(), L0=NULL, Y0=1){with(options,{
  stopifnot(is.numeric(L0))
  stopifnot(is.numeric(Y0))
  L = checkIt(L0, nPatches)
  Y = checkIt(Y0, nPatches)
  return(list(L=L, Y=Y, ...))
})}

change_L_inits 

#' @title Return the parameters as a list
#' 
#' @inheritParams ramp.xds::change_L_inits 
#' 
#' @return an **`xds`** object
#' @export
change_L_inits.SIS <- function(xds_obj, i=1, options=list()) {
  with(get_L_inits(xds_obj, i), with(options,{
    xds_obj = change_H(H, xds_obj, i)
    xds_obj$Xinits[[i]]$L = L 
    xds_obj$Xinits[[i]]$Y = Y 
    ... 
    return(xds_obj)
  }))}

Dynamical Terms

F_emerge

In the following, replace … with a formula to compute the effective density of infectious humans: it is population density by strata weighted by the probability of transmitting.

#' @title Size of effective infectious human population
#' @description Implements [F_emerge] for the SIS model.
#' @inheritParams ramp.xds::F_emerge
#' @return a [numeric] vector of length `nPatches`
#' @export
F_emerge.newLname <- function(y, xds_obj, i) {
  with(get_L_vars(y, xds_obj, i), 
    with(xds_obj$L_obj[[i]], {
      emerge = ...  
      return(emerge)
}))}

F_L

In the following, replace … with a formula to compute the effective density of infectious humans: it is population density by strata weighted by the probability of transmitting.

#' @title Size of effective infectious human population
#' @description Implements [F_L] for the SIS model.
#' @inheritParams ramp.xds::F_L
#' @return a [numeric] vector of length `nPatches`
#' @export
F_L.newLname <- function(y, xds_obj, i) {
  with(get_L_vars(y, xds_obj, i), 
    with(xds_obj$L_obj[[i]], {
      L = ...  
      return(L)
}))}

Other Outputs

To make ramp.xds easy to use, we provide some functions that parse the outputs and compute some standard metrics.

F_capacity

Replace L/H with the formula to compute prevalence.

#' @title Compute the true prevalence of infection / parasite rate
#' @description Implements F_pr for the newLname model.
#' @inheritParams F_pr
#' @return a [numeric] vector of length `nPatches`
#' @export
F_capacity.newLname <- function(vars, L_obj) {
  with(vars, 
    with(L_obj,
       return(...)
  ))
}

Internal Consistency

Optional Functions

Plotting

16. xds_lines_X

xds_lines_X_SIS 

#' Add lines for the density of infected individuals for the SIS model
#'
#' @param L a list with the outputs of parse_deout_X_SIS
#' @param nPatches the number of population strata
#' @param clrs a vector of colors
#' @param llty an integer (or integers) to set the `lty` for plotting
#'
#' @export
xds_lines_X_SIS = function(L, nPatches, clrs=c("darkblue","darkred"), llty=1){
  with(L,{
    if(nPatches==1) {
      lines(time, S, col=clrs[1], lty = llty[1])
      lines(time, I, col=clrs[2], lty = llty[1])
    }
    if(nPatches>1){
      if (length(clrs)==2) clrs=matrix(clrs, 2, nPatches)
      if (length(llty)==1) llty=rep(llty, nPatches)

      for(i in 1:nPatches){
        lines(time, S[,i], col=clrs[1,i], lty = llty[i])
        lines(time, I[,i], col=clrs[2,i], lty = llty[i])
      }
    }
  })}

xds_lines_X_newLname 

#' Add lines for the density of infected individuals for the newLname model
#'
#' @param L a list with the outputs of parse_deout_X_newLname
#' @param nPatches the number of population strata
#' @param clrs a vector of colors
#' @param llty an integer (or integers) to set the `lty` for plotting
#'
#' @export
xds_lines_X_newLname = function(L, nPatches, clrs=c("darkblue","darkred"), llty=1){
  with(L,{
    if(nPatches==1) {
      lines(time, S, col=clrs[1], lty = llty[1])
      lines(time, I, col=clrs[2], lty = llty[1])
    }
    if(nPatches>1){
      if (length(clrs)==2) clrs=matrix(clrs, 2, nPatches)
      if (length(llty)==1) llty=rep(llty, nPatches)

      for(i in 1:nPatches){
        lines(time, S[,i], col=clrs[1,i], lty = llty[i])
        lines(time, I[,i], col=clrs[2,i], lty = llty[i])
      }
    }
  })}

17. xds_plot_X

We provide a basic plotting function. The function xds_lines_X_* does not have an S3 version. Instead, it is called by the S3 function xds_plot_X

Generic 

#' Basic plotting for epidemiological models
#'
#' @param xds_obj a list that defines an `ramp.xds` model (*e.g.*,  generated by `xde_setup()`)
#' @param i the host species index
#' @param clrs a vector of colors
#' @param llty an integer (or integers) to set the `lty` for plotting
#' @param stable a logical: set to FALSE for `orbits` and TRUE for `stable_orbits`
#' @param add_axes a logical: plot axes only if TRUE
#'
#' @export
xds_plot_X = function(xds_obj, i=1, clrs="black", llty=1, stable=FALSE, add_axes=TRUE){
  UseLethod("xds_plot_X", xds_obj$L_obj[[i]])
}

xds_plot_X.SIS 

#' Plot the density of infected individuals for the newLname model
#'
#' @inheritParams ramp.xds::xds_plot_X
#' @export
xds_plot_X.SIS = function(xds_obj, i=1, clrs=c("darkblue","darkred"), llty=1, stable=FALSE, add_axes=TRUE){
  vars=with(xds_obj$outputs,if(stable==TRUE){stable_orbits}else{orbits})

  if(add_axes==TRUE)
    with(vars$L[[i]],
         plot(time, 0*time, type = "n", ylim = c(0, max(H)),
              ylab = "# Infected", xlab = "Time"))


  xds_lines_X_SIS(vars$L[[i]], xds_obj$Hpar[[i]]$nPatches, clrs, llty)
}

xds_plot_X.newLname 

#' Plot the density of infected individuals for the newLname model
#'
#' @inheritParams ramp.xds::xds_plot_X
#' @export
xds_plot_X.newLname = function(xds_obj, i=1, clrs=c("darkblue","darkred"), llty=1, stable=FALSE, add_axes=TRUE){
  vars=with(xds_obj$outputs,if(stable==TRUE){stable_orbits}else{orbits})

  if(add_axes==TRUE)
    with(vars$L[[i]],
         plot(time, 0*time, type = "n", ylim = c(0, max(H)),
              ylab = "# Infected", xlab = "Time"))

  xds_lines_X_newLname(vars$L[[i]], xds_obj$Hpar[[i]]$nPatches, clrs, llty)
}

Reusable Code

ramp.xds has a modular design. The modular design makes it possible to break down published models to serve as dynamical components in new models for malaria, a kind of plug-and-play functionality for two kinds of systems:

  • systems of ordinary and delay differential equations models for malaria and other mosquito-transmitted pathogens.

  • stochastic and deterministic discrete time systems for malaria and other mosquito-transmitted pathogens.

This package – ramp.library – is a collection of models to implement a large number of dynamical model families and other algorithms taken from the literature describing malaria and other mosquito-transmitted pathogens (see Reiner, et al. 2013)1. Those models have been implemented here as code that has undergone rigorous testing and that can be used to build models for simulation-based analytics.

The mathematics supporting this framework were explained in Spatial Dynamics of Lalaria Transmission 2

This vignette was written to help add new modules to ramp.library.

  1. Reiner RC Jr, Perkins TA, Barker CL, Niu T, Chaves LF, Ellis AL, et al. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. J R Soc Interface. 2013;10: 20120921.↩︎

  2. Wu SL, Henry JL, Citron DT, Lbabazi Ssebuliba D, Nakakawa Nsumba J, Sánchez C. HM, et al. (2023) Spatial dynamics of malaria transmission. PLoS Comput Biol 19(6): e1010684. https://doi.org/10.1371/journal.pcbi.1010684↩︎