Adding new X models to `ramp.xds`

This vignette discusses how to add a new model family to ramp.library that can be used by ramp.xds. This vignette covers models of parasite infection dynamics and immunity in humans or other hosts – the dynamical components we have denoted $\cal X$ and $\cal H.$

We have written this vignette to produce a new template file using knitr::purl, which can be modified to build a new model. Simply execute the following command, to create Xmod-newXname.R. Then, rename Xmod-newXname.R by replacing newXname in the filename. Then open the file and do a global search for newXname and replace it with your new name.

knitr::purl("X-How_To_Add_Model.Rmd", 
            output = "Xmod-newXname.R", 
            documentation = 0, 
            quiet=TRUE)

This is a good start, but there’s more work to do. For more details, read the text that follows.

This basic introduction is lacking some information that typical users might find useful. In particular, we have defined best practices encoding the $\cal H$ component of models. A vignette is planned.

Reusable Code

This package – ramp.library – contains reusable code that has been rigorously tested and that implements 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)¹. The supporting code was designed to be modular, and plug-and-play. The modular design makes it possible to break down published models to serve as the dynamical components in new models for malaria. The model families in ramp.library supports nimble model building that can be used by:

ramp.xds implements a modular, flexible, and extensible framework for building systems of ordinary and delay differential equations models for malaria and other mosquito-transmitted pathogens. The mathematic framework for ramp.xds was explained in Spatial Dynamics of Malaria Transmission ²
ramp.dts implements a modular, flexible, and extensible framework for building discrete time systems, including both deterministic and stochastic difference equations for malaria and other mosquito-transmitted pathogens.

This vignette was written to help add new models to ramp.xds.

In planned companion to this vignette, we show how to add models of $\cal XH$ to ramp.dts.

Functions

A new model family of class $\cal XH$ – human infection dynamics, including immunity – has 17 required functions encoded using a standard template. Some of these required functions are S3 class functions in human-interface.R, and several examples are posted in the ramp.xds github repository, including the SIS model human-SIS.R.

To configure a new model, you must choose a string called newXname. The only rule for choosing newXname is that it must be unique. No other ${\cal XH}$ class model in ramp.xds or ramp.library can be using that string already. All the S3 functions get dispatched by the string that replaces newXname.

The 17 required functions are:

Derivatives

dXdt.newXname - differential equations are defined by a system of differential equations. in ramp.xds these are encoded in a function called dXdt that computes and returns the derivatives. The function is set up so that the system can be solved by deSolve::ode or deSolve::dede.

Variables - A set of functions is set up that document the variables by name and assign each one an index, to store, assign and update initial values. These functions streamline setup, and guarantee consistency.

create_Xinits_newXname - each model requires a function that documents the variables by name and assigns initial values
make_Xinits.newXname - is a wrapper, that gets called by ramp.xds::xde_setup. It gets called by create_Xinits_newXname, and for the $i^{th}$ host, the initial values are stored as pars$Xinits[[i]].
make_indices_X.newXname - is the function that assigns an index to each variable in the model, and stores it as pars$ix$X[[i]]. The indices are returned as a named list.
list_Xvars.newXname - retrieves the value of variables at a point in time and returns the values by name in a list; the function gets called by dXdt and by update_Xinits and it can be useful in other contexts.
get_inits_X.newXname - retrieves the stored initial values in the same order that they are returned by dXdt.newXname
update_Xinits.newXname - a utility that makes it possible to extract and store a new set of initial values for newXname from a vector $y$ holding all the variables in a model.

Parameters

create_Xpar_newXname returns a formatted list:
- the parameter values are stored by name
- class(pars) = newXname
make_Xpar.newXname is a wrapper that calls create_Xpar_newXname

Dynamical Terms are Computed using a Standard Set of Functions

F_X.newXname - compute the effective infective density of the vertebrate hosts
F_H.newXname - compute the human population density for the host strata
F_b.newXname - compute the probability a most will become infected after receiving an infective bite, which combines elements of this $\cal X$ model with a model of environmental heterogeneity in Exposure.

Outputs

parse_deout_X.newXname - this function parses the output of deSolve and returns a named list that is either a vector or a matrix holding the values of the variables at each point in time. It works like list_Xvars, which accepts a vector holding $y$ , the values of the variables at a single point in time.
F_pr.newXname - this function takes
HTC.newXname

Outputs

xds_lines_X.newXname
xds_plot_X.newXname is a wrapper that calls xds_lines_X

Style

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 Xpar
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.

Derivatives

1. dXdt

The function dXdt is established as a generic function that dispatches on Xpar. Since there might be many host species, it actually dispatches on the $i^{th}$ element in a list, pars$Xpar[[i]].

#' @title Derivatives for human population
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param t current simulation time
#' @param y state vector
#' @param pars a list
#' @param i the host species index
#' @return a [numeric] vector
#' @export
dXdt <- function(t, y, pars, i) {
  UseMethod("dXdt", pars$Xpar[[i]])
}

dXdt.SIS

This is a version of the SIS model in ramp.xds with demographics. If a model has any demographic processes (e.g. aging, migration, births and deaths), then each variable’s derivative must be transformed by the same matrix as the human population density, H. Any disease induced mortality should be included in a matrix that computes changes in human population density, dHdt.

In the following, we use the with() function to expose parameter names attached to lists to make the code easy to read.

#' @title Derivatives for human population
#' @description Implements [dXdt] for the SIS model, no demography.
#' @inheritParams dXdt
#' @return a [numeric] vector
#' @export
dXdt.SIS <- function(t, y, pars, i) {

  foi <- pars$FoI[[i]]

  with(list_Xvars(y, pars, i),{
    H <- F_H(t, y, pars, i)
    with(pars$Xpar[[i]], {
      dS <- Births(t, H, pars, i) - foi*S + r*I + dHdt(t, S, pars, i)
      dI <- foi*S - r*I + dHdt(t, I, pars, i)
      return(c(dS, dI))
    })
  })
}

dXdt.newXname

This should be straightforward, but there’s one very important wrinkle. If you want to be able to use the ramp.xds capabilities for cohort dynamics and vital dynamics, then the dynamics must be extended carefully to call the demographic matrix function dHdt. This can become tricky, so we handle it elsewhere. When writing the code:

The foi gets computed in Exposure and stored as pars$FoI[[i]] for the $i^{th}$ host species.
In the next step, the variables are extracted from the variable vector, y, using the indices, which are stored in pars$ix$X[[i]]
We compute the derivatives.
We return the variables in order. Be sure to establish and then maintain a convention for ordering the variables in your model, and then maintain that order in every function. In particular:
- set up the initial values in that order
- return the derivatives in that order
Be careful about how you nest the with(*, ...) calls.

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

  # do not change this
  foi <- pars$FoI[[i]]

  # attach the variables by name 
  with(list_Xvars(y, pars, i),{
    # compute H (if it isn't one of the variables) 
    H <- F_H(t, y, pars, i)

    # expose the parameters (see create_Xpar_Xname) 
    with(pars$Xpar[[i]], {
      # compute the derivatives 
      dX1 <- ... 
      dX2 <- ...
      ... 
      
      # concatenate the derivatives 
      derivs = c(dX1, dX2, ...) 
      
      # return the derivatives 
      return(derivs)
    })
  })
}

Variables

2. create_Xinits

create_Xinits_SIS

#' @title Make initial values for the SIS human model, with defaults
#' @param nStrata the number of strata in the model
#' @param Xopts a [list] to overwrite defaults
#' @param H0 the initial human population density
#' @param S0 the initial values of the parameter S
#' @param I0 the initial values of the parameter I
#' @return a [list]
#' @export
create_Xinits_SIS = function(nStrata, Xopts = list(), H0=NULL, S0=NULL, I0=1){with(Xopts,{
  if(is.null(S0)) S0 = H0 - I0
  stopifnot(is.numeric(S0))
  S = checkIt(S0, nStrata)
  I = checkIt(I0, nStrata)
  return(list(S=S, I=I))
})}

create_Xinits_newXname

#' @title Make initial values for the Xname human model, with defaults
#' @param nStrata the number of strata in the model
#' @param Xopts a [list] to overwrite defaults
#' @param X10 the initial value for X1 
#' @param X20 the initial value for X1 
#' @return a [list]
#' @export
create_Xinits_newXname = function(nStrata, Xopts = list(), X10=NULL, X20=1){with(Xopts,{
  stopifnot(is.numeric(X10))
  stopifnot(is.numeric(X20))
  X1 = checkIt(X10, nStrata)
  X2 = checkIt(X20, nStrata)
  return(list(X1=X1, X2=X2, ...))
})}

3. make_Xinits

In the xde_setup family of setup functions, initial values for a model fulfilling the dynamical component $\cal X,$ called newXname is set up by the function make_Xinits.

Generic

The make_Xinits function dispatches on pars$Xpar[[i]]

#' @title A function to set up Xpar
#' @description This method dispatches on `newXname`.
#' @param pars a [list]
#' @param i the host species index
#' @param Xopts a [list]
#' @return a [list]
#' @export
make_Xinits = function(pars, i, Xopts=list()){
  UseMethod("make_Xinits", pars$Xpar[[i]])
}

make_Xinits.SIS

The make_Xinits.SIS calls create_Xinits_SIS

#' @title Setup Xinits.SIS
#' @description Implements [make_Xinits] for the SIS model
#' @inheritParams make_Xinits
#' @return a [list] vector
#' @export
make_Xinits.SIS = function(pars, i, Xopts=list()){
  pars$Xinits[[i]] = with(pars, create_Xinits_SIS(pars$Hpar[[i]]$nStrata, Xopts, H0=Hpar[[i]]$H))
  return(pars)
}

make_Xinits.newXname

To modify, replace newXname

#' @title Setup Xinits.newXname
#' @description Implements [make_Xinits] for the newXname model
#' @inheritParams ramp.xds::make_Xinits
#' @return a [list] vector
#' @export
make_Xinits.newXname = function(pars, i, Xopts=list()){
  pars$Xinits[[i]] = with(pars, create_Xinits_newXname(pars$Hpar[[i]]$nStrata, Xopts, H0=Hpar[[i]]$H))
  return(pars)
}

4. make_indices_X

Generic

#' @title Add indices for human population to parameter list
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param pars a [list]
#' @param i the host species index
#' @return a [list]
#' @export
make_indices_X <- function(pars, i) {
  UseMethod("make_indices_X", pars$Xpar[[i]])
}

make_indices_X.SIS

#' @title Add indices for human population to parameter list
#' @description Implements [make_indices_X] for the SIS model.
#' @inheritParams ramp.xds::make_indices_X
#' @return none
#' @importFrom utils tail
#' @export
make_indices_X.SIS <- function(pars, i) {with(pars,{

  S_ix <- seq(from = max_ix+1, length.out=Hpar[[i]]$nStrata)
  max_ix <- tail(S_ix, 1)

  I_ix <- seq(from = max_ix+1, length.out=Hpar[[i]]$nStrata)
  max_ix <- tail(I_ix, 1)

  pars$max_ix = max_ix
  pars$ix$X[[i]] = list(S_ix=S_ix, I_ix=I_ix)
  return(pars)
})}

make_indices_X.newXname

#' @title Add indices for human population to parameter list
#' @description Implements [make_indices_X] for the newXname model.
#' @inheritParams make_indices_X
#' @return none
#' @importFrom utils tail
#' @export
make_indices_X.newXname <- function(pars, i) {with(pars,{

  X1_ix <- seq(from = max_ix+1, length.out=Hpar[[i]]$nStrata)
  max_ix <- tail(X1_ix, 1)

  X2_ix <- seq(from = max_ix+1, length.out=Hpar[[i]]$nStrata)
  max_ix <- tail(X2_ix, 1)

  ... 
  
  pars$max_ix = max_ix
  pars$ix$X[[i]] = list(X1_ix=X1_ix, X2_ix=X2_ix, ...)
  return(pars)
})}

5. list_Xvars

Generic

#' @title Return the variables as a list
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param y the variables
#' @param pars a [list]
#' @param i the host species index
#' @return a [list]
#' @export
list_Xvars <- function(y, pars, i) {
  UseMethod("list_Xvars", pars$Xpar[[i]])
}

list_Xvars.SIS

#' @title Return the variables as a list
#' @description This method dispatches on the type of `pars$Xpar`
#' @inheritParams list_Xvars
#' @return a [list]
#' @export
list_Xvars.SIS <- function(y, pars, i) {
  with(pars$ix$X[[i]],{
      S = y[S_ix]
      I = y[I_ix]
      H = S+I 
      return(list(S=S,I=I,H=H))})
}

list_Xvars.newXname

#' @title Return the variables as a list
#' @description This method dispatches on the type of `pars$Xpar`
#' @inheritParams ramp.xds::list_Xvars
#' @return a [list]
#' @export
list_Xvars.newXname <- function(y, pars, i) {
  with(pars$ix$X[[i]],{
      X1 = y[X1_ix]
      X2 = y[X2_ix]
      ...
      H = X1+X2+... 
      return(list(X1=X1,X2=X2,...,H=H))})
}

6. get_inits_X

By now the pattern should be obvious.

Generic

#' @title Return initial values as a vector
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param pars a [list]
#' @param i the host species index
#' @return none
#' @export
get_inits_X <- function(pars, i) {
  UseMethod("get_inits_X", pars$Xpar[[i]])
}

get_inits_X.SIS

Returns the initial values, in order, by name.

#' @title Return initial values as a vector
#' @description This method dispatches on the type of `pars$Xpar`.
#' @inheritParams get_inits_X
#' @return a [numeric] vector
#' @export
get_inits_X.SIS <- function(pars, i){
  with(pars$Xinits[[i]], return(c(S,I)))
}

get_inits_X.newXname

Write the function to return the initial values, in order, by name.

#' @title Return initial values as a vector
#' @description This method dispatches on the type of `pars$Xpar`.
#' @inheritParams ramp.xds::get_inits_X
#' @return a [numeric] vector
#' @export
get_inits_X.newXname <- function(pars, i){
  with(pars$Xinits[[i]], return(c(X1,X2)))
}

7. update_Xinits

Pass a vector y, extract the values of the variables, and modify the initial values.

Generic

#' @title Set the initial values from a vector of states
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param pars a [list]
#' @param y0 a vector of initial values
#' @param i the host species index
#' @return none
#' @export
update_Xinits <- function(pars, y0, i) {
  UseMethod("update_Xinits", pars$Xpar[[i]])
}

update_Xinits.SIS

#' @title Update inits for the SIS human model from a vector of states
#' @inheritParams update_Xinits
#' @return none
#' @export
update_Xinits.SIS <- function(pars, y0, i) {
  with(list_Xvars(y0, pars, i),{
    pars$Xinits[[i]] = create_Xinits_SIS(pars, list(), S0=S, I0=I)
    return(pars)
})}

update_Xinits.newXname

#' @title Update inits for the newXname human model from a vector of states
#' @inheritParams ramp.xds::update_Xinits 
#' @return none
#' @export
update_Xinits.newXname <- function(pars, y0, i) {
  with(list_Xvars(y0, pars, i),{
    pars = create_Xinits_newXname(pars, list(), X10=X1, X20=X2, ...)
    return(pars)
})}

Parameters

The second function sets up a list Xpar that dispatches dXdt and that has the parameter values.

8. create_Xpar_SIS

There is no S3 definition for the create_Xpar_* class of functions. Instead, we write a generic S3 function called make_Xpar that calls a function make_Xpar_* There is, nevertheless, a reasonably standard formula for make_Xpar_*

The first argument is nStrata the number of population strata for the species.
The first argument is Xopts(). Notice that Xopts() is an empty list by default, but if a non-empty list is passed, the values overwrite the defaults.
The parameters are attached by name. These parameter names must be the same ones used in dXdt
By default, they should have length(*)= nStrata
The last part of the function name is the string that dispatches the S3 functions. In create_Xpar_SIS, we write class(Xpar) <- "SIS"

create_Xpar_SIS

#' @title Make parameters for SIS human model, with defaults
#' @param nStrata is the number of population strata
#' @param Xopts a [list] that could overwrite defaults
#' @param b transmission probability (efficiency) from mosquito to human
#' @param c transmission probability (efficiency) from human to mosquito
#' @param r recovery rate
#' @return a [list]
#' @export
create_Xpar_SIS = function(nStrata, Xopts=list(),
                         b=0.55, r=1/180, c=0.15){
  with(Xopts,{
    Xpar = list()
    class(Xpar) <- "SIS"

    Xpar$b = checkIt(b, nStrata)
    Xpar$c = checkIt(c, nStrata)
    Xpar$r = checkIt(r, nStrata)

    return(Xpar)
  })}

create_Xpar_newXname

To format this for a new function, do a search and replace on newXname and replace the parmeter names and values $p1$ , $p2$ , $\ldots$

#' @title Make parameters for newXname human model, with defaults
#' @param nStrata is the number of population strata
#' @param Xopts 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
create_Xpar_newXname = function(nStrata, Xopts=list(),
                         p1=1, p2=2, p3=3){
  with(Xopts,{
    Xpar = list()
    class(Xpar) <- c("newXname")

    Xpar$p1 = checkIt(p1, nStrata)
    Xpar$p2 = checkIt(p2, nStrata)
    Xpar$p3 = checkIt(p3, nStrata)

    return(Xpar)
  })}

9. make_Xpar

In the xde_setup family of setup functions, setup for a model fulfilling the dynamical component $\cal X,$ called newXname is set up by the function make_Xpar.

Generic

This function must set up Xpar and it must also set up the initial values. By convention, the functions that get called are not S3 class functions.

The make_Xpar functions assigns to Xname the class(Xname) and dispatches on Xname

#' @title A function to set up Xpar
#' @description This method dispatches on `Xname`.
#' @param Xname a [character] string
#' @param pars a [list]
#' @param i the host species index
#' @param Xopts a [list]
#' @return a [list]
#' @export
make_Xpar = function(Xname, pars, i, Xopts=list()){
  class(Xname) <- Xname
  UseMethod("make_Xpar", Xname)
}

make_Xpar.SIS

The dispatched function SIS ignores xname, but it calls a function create_Xpar_SIS

#' @title Setup Xpar.SIS
#' @description Implements [make_Xpar] for the SIS model
#' @inheritParams make_Xpar
#' @return a [list] vector
#' @export
make_Xpar.SIS = function(Xname, pars, i, Xopts=list()){
  pars$Xpar[[i]] = create_Xpar_SIS(pars$Hpar[[i]]$nStrata, Xopts)
  return(pars)
}

make_Xpar.newXname

To use this, simply to a search and replace on newXname.

#' @title Setup Xpar.newXname
#' @description Implements [make_Xpar] for the newXname model
#' @inheritParams ramp.xds::make_Xpar
#' @return a [list] vector
#' @export
make_Xpar.newXname = function(newXname, pars, i, Xopts=list()){
  pars$Xpar[[i]] = create_Xpar_newXname(pars$Hpar[[i]]$nStrata, Xopts)
  return(pars)
}

Terms

All models in ramp.xds the modular design dictates that we must define a few functions that compute a few key quantities:

X is defined to be the effective density of infectious humans, retrieved by F_X
H is defined to be the size of each population stratum, retrieved by F_H
b is defined to be the fraction of infective bites that cause an infection, retrieved by F_b

10. F_X

The function F_X returns the effective density of infectious humans for each stratum. The notion is that $X$ is proportional to $H$ , where the constant is the fraction of bites on the strata that would infect a mosquito.

Generic

#' @title Size of effective infectious human population
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param y state vector
#' @param pars a list
#' @param i the host species index
#' @return a [numeric] vector of length `nStrata`
#' @export
F_X <- function(y, pars, i) {
  UseMethod("F_X", pars$Xpar[[i]])
}

F_X.SIS

As an example, we include the $SIS$ model from ramp.xds. In the SIS model, the returned value is $cI$ where $I$ is the density of infected individuals; the prevalence of infection is $X/H$ . The constant $c$ is the probability a mosquito would become infected after blood feeding on an infected human. Two things to note:

First, ramp.xds holds the indices by component under pars$ix, so pars$ix$X holds the indices for all the variables for all species, and pars$ix$X[[i]] holds the indices for the $i^{th}$ species. Putting them under X helps to avoid conflicts if, for example, another dynamical component happened to use the same variable name.

#' @title Size of effective infectious human population
#' @description Implements [F_X] for the SIS model.
#' @inheritParams F_X
#' @return a [numeric] vector of length `nStrata`
#' @export
F_X.SIS <- function(y, pars, i) {
  I = y[pars$ix$X[[i]]$I_ix]
  X = with(pars$Xpar[[i]], c*I)
  return(X)
}

The commented code relevant code that must be changed is commented out. The preferred style is to extract the relevant variables by name, then compute the effective infectious density,or individuals in this stratum, then return the value. In the manuscript $^1$ , this quantity was called X.

F_X.newXname

The following is a template for new model. The declaration of inheritParams from the SIS model in ramp.xds works naturally. The namespace from ramp.xds is needed, so our newXname function (above) has ramp.xds::F_X.

The steps are to extract the variables in

#' @title Size of effective infectious human population
#' @description Implements [F_X] for the SIS model.
#' @inheritParams ramp.xds::F_X
#' @return a [numeric] vector of length `nStrata`
#' @export
F_X.newXname <- function(y, pars, i) {
  ########################
  # extract: 
  # VAR <- y[pars$ix$X$...] 
  ########################
  with(pars$Xpar[[i]], 
       ########################
       # compute: 
       # X <- ... F(VAR) 
       ########################
  )
  return(X)
}

11. F_H

The total population size, $H$ , is often important, so we write a function that computes it.

Generic

#' @description This method dispatches on the type of `pars$Xpar`.
#' @param y state vector
#' @param pars a list
#' @param i the host species index
#' @return a [numeric] vector of length `nStrata`
#' @export
F_H <- function(y, pars, i) {
  UseMethod("F_H", pars$Xpar[[i]])
}

F_H.SIS

This is the function that computes $H = S+I$ for the SIS model.

#' @title Size of effective infectious human population
#' @description Implements [F_H] for the SIS model.
#' @inheritParams F_H
#' @return a [numeric] vector of length `nStrata`
#' @export
F_H.SIS <- function(y, pars, i){
  with(list_Xvars(y, pars, i), return(H))
}

F_H.newXname

The function F_X returns the effective density of infectious humans for each stratum. The following is a template for a model, assuming $H$ is computed in list_Xvars:

#' @title Size of effective infectious human population
#' @description Implements [F_H] for the newXname model.
#' @inheritParams ramp.xds::F_H
#' @return a [numeric] vector of length `nStrata`
#' @export
F_H.newXname <- function(y, pars, i){
  with(list_Xvars(y, pars, i), return(H))
}

12. F_b

We use the short parameter name b to denote the fraction of infective bites by mosquitoes that cause an infection. This is computed in another function, called Exposure under a model for environmental heterogeneity that computes a local FoI based on the local daily EIR. Then, the total FoI is computed from a travel model.

Generic

#' @title Infection blocking pre-erythrocytic immunity
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param y state vector
#' @param pars a list
#' @param i the host species index
#' @return a [numeric] vector of length `nStrata`
#' @export
F_b <- function(y, pars, i) {
  UseMethod("F_b", pars$Xpar[[i]])
}

`F_b.SIS`

In ramp.xds, we recognize the problem of translating the local EIR into the FoI. The local EIR is computed as one of the core model terms for each one of the model strata.

The function F_b returns the model-defined probability an infective bite by a mosquito would cause an infection. This could include the combined effects of inefficient transmission, pre-erytrocytic immunity, and perhaps some effects of blood stage immunity. In a model with no effect, the return value would be $1$ .

Once again, we provide the SIS model from ramp.xds as our example:

#' @title Infection blocking pre-erythrocytic immunity
#' @description Implements [F_b] for the SIS model.
#' @inheritParams F_b
#' @return a [numeric] vector of length `nStrata`
#' @export
F_b.SIS <- function(y, pars, i) {
  return(with(pars$Xpar[[i]], b))
}

F_b.newXname

Depending on the model, the value of b is either retrieved or it is computed from other variables in the model, and it is expected that length(b) = nStrata

#' @title Infection blocking pre-erythrocytic immunity
#' @description Implements [F_b] for the newXname model.
#' @inheritParams ramp.xds::F_b
#' @return a [numeric] vector of length `nStrata`
#' @export
F_b.newXname <- function(y, pars, i) {
  with(pars$Xpar[[i]],{ 
    ########################
    # retrieve or compute it 
    ########################
    b = pars$Xpar[[i]]$b
    ########################
    # return it 
    ########################
    return(b)
  })
}

Outputs

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

13. `parse_deout_X`

parse_deout_X takes the output of deSolve, parses the outputs using the variable indices, and returns the variables by name as a list.

Generic

#' @title Parse the output of deSolve and return the variables by name in a list
#' @description This method dispatches on the type of `pars$Xpar`. Adds the variables
#' from the X model to a list and returns it
#' @param deout a [matrix] of outputs from deSolve
#' @param pars a [list] that defines a model
#' @param i the host species index
#' @export
parse_deout_X <- function(deout, pars, i) {
  UseMethod("parse_deout_X", pars$Xpar[[i]])
}

parse_deout_X.SIS

#' @title Parse the output of deSolve and return variables for the SIS model
#' @description Implements [parse_deout_X] for the SIS model
#' @inheritParams parse_deout_X
#' @return none
#' @export
parse_deout_X.SIS <- function(deout, pars, i) {
  time = deout[,1]
  with(pars$ix$X[[i]],{
    S = deout[,S_ix+1]
    I = deout[,I_ix+1]
    H = S+I
    return(list(time=time, S=S, I=I, H=H))
})}

parse_deout_X.newXname

#' @title Parse the output of deSolve and return variables for the newXname model
#' @description Implements [parse_deout_X] for the newXname model
#' @inheritParams ramp.xds::parse_deout_X
#' @return none
#' @export
parse_deout_X.newXname <- function(deout, pars, i) {
  time = deout[,1]
  with(pars$ix$X[[i]],{
    X1 = deout[,X1_ix+1]
    X2 = deout[,X2_ix+1]
    ... 
    H = X1 + X2 + ...
    return(list(time=time, X1=X1, X2=X2, ..., H=H))
})}

14. F_pr

We require each model to compute the true PR.

Generic

#' @title Compute the *true* prevalence of infection / parasite rate
#' @description This method dispatches on the type of `pars$Xpar[[i]]`.
#' @param vars a list with the variables attached by name
#' @param Xpar a list defining a model for human
#' @return a [numeric] vector of length `nStrata`
#' @export
F_pr <- function(vars, Xpar) {
  UseMethod("F_pr", Xpar)
}

F_pr.SIS

#' @title Compute the "true" prevalence of infection / parasite rate
#' @description Implements F_pr for the SIS model.
#' @inheritParams F_pr
#' @return a [numeric] vector of length `nStrata`
#' @export
F_pr.SIS <- function(vars, Xpar) {
  pr = with(vars, I/H)
  return(pr)
}

The function F_pr computes the models true prevalence using the variable names. It is passed the varslist that is created by parse_deout (see above).

F_pr.newXname

#' @title Compute the true prevalence of infection / parasite rate
#' @description Implements F_pr for the newXname model.
#' @inheritParams F_pr
#' @return a [numeric] vector of length `nStrata`
#' @export
F_pr.newXname <- function(vars, Xpar) {
  pr = with(vars, with(Xpar,(return(X1/H))))
  return(pr)
}

15. HTC

Similarly, we require each function to return the human transmitting capacity, or HTC

Generic

#' @title Compute the human transmitting capacity
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param pars a [list]
#' @param i the host species index
#' @return none
#' @export
HTC <- function(pars, i) {
  UseMethod("get_inits_X", pars$Xpar[[i]])
}

HTC.SIS

#' @title Compute the HTC for the SIS model
#' @description Implements [HTC] for the SIS model with demography.
#' @inheritParams ramp.xds::HTC
#' @return a [numeric] vector
#' @export
HTC.SIS <- function(pars, i) {
  with(pars$Xpar[[i]],
    return(c/r)
  )
}

HTC.newXname

#' @title Compute the HTC for the newXname model
#' @description Implements [HTC] for the newXname model with demography.
#' @inheritParams ramp.xds::HTC
#' @return a [numeric] vector
#' @export
HTC.newXname <- function(pars, i) {
  with(pars$Xpar[[i]],
    #HTC <- 
    return(HTC)
  )
}

Plotting

16. xds_lines_X

xds_lines_X_SIS

#' Add lines for the density of infected individuals for the SIS model
#'
#' @param XH a list with the outputs of parse_deout_X_SIS
#' @param nStrata 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(XH, nStrata, clrs=c("darkblue","darkred"), llty=1){
  with(XH,{
    if(nStrata==1) {
      lines(time, S, col=clrs[1], lty = llty[1])
      lines(time, I, col=clrs[2], lty = llty[1])
    }
    if(nStrata>1){
      if (length(clrs)==2) clrs=matrix(clrs, 2, nStrata)
      if (length(llty)==1) llty=rep(llty, nStrata)

      for(i in 1:nStrata){
        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_newXname

#' Add lines for the density of infected individuals for the newXname model
#'
#' @param XH a list with the outputs of parse_deout_X_newXname
#' @param nStrata 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_newXname = function(XH, nStrata, clrs=c("darkblue","darkred"), llty=1){
  with(XH,{
    if(nStrata==1) {
      lines(time, S, col=clrs[1], lty = llty[1])
      lines(time, I, col=clrs[2], lty = llty[1])
    }
    if(nStrata>1){
      if (length(clrs)==2) clrs=matrix(clrs, 2, nStrata)
      if (length(llty)==1) llty=rep(llty, nStrata)

      for(i in 1:nStrata){
        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 pars 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(pars, i=1, clrs="black", llty=1, stable=FALSE, add_axes=TRUE){
  UseMethod("xds_plot_X", pars$Xpar[[i]])
}

xds_plot_X.SIS

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

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


  xds_lines_X_SIS(vars$XH[[i]], pars$Hpar[[i]]$nStrata, clrs, llty)
}

xds_plot_X.newXname

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

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

  xds_lines_X_newXname(vars$XH[[i]], pars$Hpar[[i]]$nStrata, clrs, llty)
}

Reiner RC Jr, Perkins TA, Barker CM, Niu T, Chaves LF, Ellis AM, et al. A systematic review of mathematical models of mosquito-borne pathogen transmission: 1970-2010. J R Soc Interface. 2013;10: 20120921.↩︎
Wu SL, Henry JM, Citron DT, Mbabazi 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 ↩︎

Adding new X models to ramp.xds

Human/Vertebrate Host Infection Dynamics

Reusable Code

Functions

Style

Derivatives

1. dXdt

dXdt.SIS

dXdt.newXname

Variables

2. create_Xinits

create_Xinits_SIS

create_Xinits_newXname

3. make_Xinits

Generic

make_Xinits.SIS

make_Xinits.newXname

4. make_indices_X

Generic

make_indices_X.SIS

make_indices_X.newXname

5. list_Xvars

Generic

list_Xvars.SIS

list_Xvars.newXname

6. get_inits_X

Generic

get_inits_X.SIS

get_inits_X.newXname

7. update_Xinits

Generic

update_Xinits.SIS

update_Xinits.newXname

Parameters

8. create_Xpar_SIS

create_Xpar_SIS

create_Xpar_newXname

9. make_Xpar

Generic

make_Xpar.SIS

make_Xpar.newXname

Terms

10. F_X

Generic

F_X.SIS

F_X.newXname

11. F_H

Generic

F_H.SIS

F_H.newXname

12. F_b

Generic

F_b.SIS

F_b.newXname

Outputs

13. parse_deout_X

Generic

parse_deout_X.SIS

parse_deout_X.newXname

14. F_pr

Generic

F_pr.SIS

F_pr.newXname

15. HTC

Generic

HTC.SIS

HTC.newXname

Plotting

16. xds_lines_X

xds_lines_X_SIS

xds_lines_X_newXname

17. xds_plot_X

Generic

xds_plot_X.SIS

xds_plot_X.newXname

`F_b.SIS`

13. `parse_deout_X`