How to add a new model to the library: Human Infection

Required Functions

The package motrap extends the package exDE, which implements a modular, flexible, and extensible framework for building ordinary and delay differential equations models for malaria and other mosquito-transmitted pathogens. The underlying mathematics has been explained in Spatial Dynamics of Malaria Transmission ¹ One of the most important functions of motrap is to host the model families in the model libraries for the main dynamical components. This vignette discusses how to add a new model family describing human infection dynamics and immunity.

Adding a new model to the dynamical component of class $\cal X$ – human infection dynamics, including immunity – follows a standard template. The S3 class functions are defined in human-interface.R, and several examples are posted in the exDE github repository, including the SIS model human-SIS.R.

This code can serve as a template for a new model. Simply execute the following command, and start editing Xmod-NEW.R

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

Terms, etcetera

`F_X`

The function F_X returns the effective density of infectious humans for each stratum. The following is a template for a model:

#' @title Size of effective infectious human population
#' @description Implements [F_X] for the NEW model.
#' @inheritParams exDE::F_X
#' @return a [numeric] vector of length `nStrata`
#' @export
F_X.NEW <- function(t, y, pars) {
  ########################
  # extract: 
  # VAR <- y[pars$ix$X$...] 
  ########################
  with(pars$Xpar, 
       ########################
       # compute: 
       # X <- ... F(VAR) 
       ########################
  )
  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.

Example

As an example, we include the $SIS$ model from exDE. In the SIS model, the returned value is $cX$ where prevalence 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, exDE holds the indices by component under pars$ix, so pars$ix$X holds the index names for all the variable names. Putting them under X helps to avoid conflicts if, for example, another dynamical component happened to use the same variable name.
The declaration of inheritParams from the SIS model in exDE works naturally. The namespace from exDE is needed, so our NEW function (above) has exDE::F_X

#' @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(t, y, pars) {
  X = y[pars$ix$X$X_ix]
  with(pars$Xpar, c*X)
}

`F_X`

The function F_X returns the effective density of infectious humans for each stratum. The following is a template for a model:

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

`F_b`

In exDE, 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$.

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

Once again, we provide the SIS model from exDE 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) {
  with(pars$Xpar, b)
}

Dynamics

The style guidelines are to:

Extract the variables by name using the indices attached to Xpar
Use with(pars$Xpar, to attach the parameter values by name
Compute the derivatives using human readable formulas.
Return the derivatives as a concatenated vector

This should be straightforward, but there’s one very important wrinkle. If you want to be able to use the exDE capabilities for cohort dynamics and vital dynamics, then the dynamics must be extended carefully to call the demographic matrix function dHdt.

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

  ###############################
  # get variables by name from y 
  # ...  
  ###############################

  with(pars$Xpar, {
    ###############################
    # Compute the derivatives 
    # 
    # dX <- FoI*(H - X) - r*X + dHdt(t, X, pars)
    # dH <- Births(t, H, pars) + dHdt(t, H, pars) 
    ###############################
    
    ###############################
    # Return the derivatives 
    # ...  
    ###############################
    return(c(dX1, dX2, dX3, dH))
  })
}

`dXdt`

In exDE, we write two versions of the derivatives function, to make it easy for the average user to scale up.

dXdt.SISdx <- function(t, y, pars, FoI) {
  with(pars$Xpar, {
    X <- y[X_ix]
    H <- F_H(t, y, pars)

    dX <- FoI*(H - X) - r*X

    return(c(dX))
  })
}

`dXHdt`

If the model will have human demography, 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 the derivative describing changes in human population density, dH.

dXdt.SISdH <- function(t, y, pars, FoI) {
  with(pars$Xpar, {

    H <- F_H(t, y, pars)
    X <- y[X_ix]

    dX <- FoI*(H - X) - r*X + dHdt(t, X, pars)
    dH <- Births(t, H, pars) + dHdt(t, H, pars)

    return(c(dX, dH))
  })
}

Setup

`setup_X`

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.

#' @title Setup Xpar.NEW
#' @description Implements [setup_X] for the NEW model
#' @inheritParams exDE::setup_X
#' @return a [list] vector
#' @export
setup_X.NEW = function(pars, Xname, Xopts=list()){

  pars$Xname = "NEW"
  pars = make_Xpar_NEW(pars, Xopts)
  pars = make_Xinits_NEW(pars, Xopts)

  return(pars)
}

setup_X.SIS = function(pars, Xname, Xopts=list()){
  pars$Xname = "SIS"
  pars = make_Xpar_SIS(pars, Xopts)
  pars = make_Xinits_SIS(pars, Xopts)

  return(pars)
}

`make_Xpar_...`

#' @title Make parameters for NEW human model, with defaults
#' @param pars a [list]
#' @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
make_Xpar_NEW = function(pars, Xopts=list(),
                         p1=1, p2=2, p3=3){
  with(Xopts,{
    Xpar = list()
    class(Xpar) <- c("NEW")

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

    pars$Xpar = Xpar
    return(pars)
  })}

The function that makes parameters should assign

make_Xpar_SIS = function(pars, Xopts=list(),
                         b=0.55, r=1/180, c=0.15){
  with(Xopts,{
    Xpar = list()
    class(Xpar) <- c("SISdX", "SISdXdH")

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

    pars$Xpar = Xpar
    return(pars)
})}

`make_Xinits_...`

#' @title Make initial values for the NEW human model, with defaults
#' @param pars a [list]
#' @param Xopts a [list] to overwrite defaults
#' @param X10 the initial values of the parameter X1
#' @param X20 the initial values of the parameter X2
#' @param X30 the initial values of the parameter X3
#' @return a [list]
#' @export
make_Xinits_NEW = function(pars, Xopts = list(), X10=1, X20=2, X30=3){with(Xopts,{
  inits = list()
  inits$X10 = checkIt(X10, pars$nStrata)
  inits$X20 = checkIt(X20, pars$nStrata)
  inits$X30 = checkIt(X30, pars$nStrata)
  pars$Xinits = inits
  return(pars)
})}

make_Xinits_SIS = function(pars, Xopts = list(), X0=1){with(Xopts,{
  inits = list()
  inits$X0 = checkIt(X0, pars$nStrata)
  pars$Xinits = inits
  return(pars)
})}

`make_indices_X`

#' @title Add indices for human population to parameter list
#' @description Implements [make_indices_X] for the NEW model.
#' @inheritParams exDE::make_indices_X
#' @return none
#' @importFrom utils tail
#' @export
make_indices_X.NEW <- function(pars) {
  
  pars$Xpar$X1_ix <- seq(from = pars$max_ix+1, length.out = pars$nStrata)
  pars$max_ix <- tail(pars$Xpar$X1_ix, 1)
  
  pars$Xpar$X2_ix <- seq(from = pars$max_ix+1, length.out = pars$nStrata)
  pars$max_ix <- tail(pars$Xpar$X2_ix, 1)
  
  pars$Xpar$X3_ix <- seq(from = pars$max_ix+1, length.out = pars$nStrata)
  pars$max_ix <- tail(pars$Xpar$X3_ix, 1)
  
  return(pars)
}

#' @title Add indices for human population to parameter list
#' @description Implements [make_indices_X] for the SIS model.
#' @inheritParams make_indices_X
#' @return none
#' @importFrom utils tail
#' @export
make_indices_X.SIS <- function(pars) {
  pars$Xpar$X_ix <- seq(from = pars$max_ix+1, length.out = pars$nStrata)
  pars$max_ix <- tail(pars$Xpar$X_ix, 1)
  return(pars)
}

Initial Values

`update_inits_X`

#' @title Update inits for the SIS human model from a vector of states
#' @inheritParams exDE::update_inits_X 
#' @return none
#' @export
update_inits_X.NEW <- function(pars, y0) {
  with(pars$ix$X,{
    X10 = y0[X1_ix]
    X20 = y0[X2_ix]
    X30 = y0[X3_ix]
    pars = make_Xinits_SIS(pars, list(), X10, X20, X30)
    return(pars)
})}

update_inits_X.SIS <- function(pars, y0) {
  X0 = y0[pars$ix$X$X_ix]
  make_Xinits_SIS(pars, list(), X0)
  return(pars)
}

`get_inits_X`

#' @title Return initial values as a vector
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param pars a [list]
#' @return none
#' @export
get_inits_X.NEW <- function(pars){
  with(pars$Xinits,{
    return(c(X1, X2, X3)) 
  })
}

#' @title Return initial values as a vector
#' @description This method dispatches on the type of `pars$Xpar`.
#' @param pars a [list]
#' @return none
#' @export
get_inits_X.SIS <- function(pars){
  pars$Xinits$X0
}

Outputs

`parse_deout_X`

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

#' @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) {
  time = deout[,1]
  Hlist <- parse_deout_H(deout, pars)
  with(Hlist,{
    X = deout[,pars$ix$X$X_ix+1]
    return(list(time=time, X=X, H=H))
})}

`F_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).

#' @title Compute the "true" prevalence of infection / parasite rate
#' @description Implements [F_pr] for the NEW model.
#' @inheritParams exDE::F_pr
#' @return a [numeric] vector of length `nStrata`
#' @export
F_pr.NEW<- function(varslist, pars) {
  pr = with(varslist$XH, X1/H)
  return(pr)
}

`HTC`

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

Plotting

`xde_plot_X`

#' Plot the density of infected individuals for the NEW model
#'
#' @inheritParams xde_plot_X
#' @export
xde_plot_X.NEW = function(pars, clrs="black", 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,
         plot(time, 0*time, type = "n", ylim = c(0, max(H)),
              ylab = "# Infected", xlab = "Time"))

  xde_lines_X(vars$XH, pars, clrs, llty)
}

`xde_lines_X`

#' Add lines for the density of infected individuals for the NEW model
#'
#' @inheritParams xde_lines_X
#'
#' @export
xde_lines_X.NEW = function(XH, pars, clrs="black", llty=1){
  with(XH,{
    if(pars$nStrata==1) lines(time, X1, col=clrs[1], lty = llty[1])
    if(pars$nStrata>1){
      if (length(clrs)==1) clrs=rep(clrs, pars$nStrata)
      if (length(llty)==1) llty=rep(llty, pars$nStrata)
      for(i in 1:pars$nStrata){
        lines(time, X1[,i], col=clrs[i], lty = llty[i])
      }
    }
  })}

February 12, 2024

Required Functions

Terms, etcetera

F_X

F_X

F_b

Dynamics

dXdt

dXHdt

Setup

setup_X

make_Xpar_...

make_Xinits_...

make_indices_X

Initial Values

update_inits_X

get_inits_X

Outputs

parse_deout_X

F_pr

HTC

Plotting

xde_plot_X

xde_lines_X

`F_X`

`F_X`

`F_b`

`dXdt`

`dXHdt`

`setup_X`

`make_Xpar_...`

`make_Xinits_...`

`make_indices_X`

`update_inits_X`

`get_inits_X`

`parse_deout_X`

`F_pr`

`HTC`

`xde_plot_X`

`xde_lines_X`