facilityepimath

an R package to analyze differential equation models of infectious disease in healthcare facilities

Damon Toth

University of Utah

Healthcare Facility Transmission

• Patient-to-patient transmission is important but difficult to document
  • Inpatient with infection onset may have acquired organism before or after admission
  • Extent to which transmission-interrupting interventions will reduce infections is unclear
• Transmission threshold effects may be especially important
  • Can a single high-transmission facility sustain/amplify transmission by itself (within-faciliy \(R_0>1\))?
  • Interventions that reduce \(R_0<1\) can be highly impactful and cost-beneficial
• Mathematical models can help address these topics

facilityepimath R package

• Provide functions for simulation-free calculations of useful quantities for a differential equation model of healthcare facility transmission

  • Facility basic reproduction number (\(R_0\))

  • Model equilibrium at a given distribution of admission states

• Allow flexible user input of heterogenous patient states and length of stay

  • Heterogeneous susceptibility and transmissibility of states

  • State-dependent discharge / death rate option

  • Time-of-stay-dependent discharge / death hazard option

Differential equation model

General model can have any number of state compartments

Generic example model with two “colonized” compartments \(C_1\) and \(C_2\) and two “susceptible” compartments \(S_1\) and \(S_2\):

\(\frac{dS_1}{dt} = -(s_{21}+(a_{11}+a_{21})\alpha+\omega_1+h(t))S_1 + s_{12}S_2 + r_{11}C_1 + r_{12}C_2\)

\(\frac{dS_2}{dt} = s_{21}S_1 - (s_{12}+(a_{12}+a_{22})\alpha+\omega_2+h(t))S_2 + r_{21}C_1 + r_{22}C_2\)

\(\frac{dC_1}{dt} = a_{11}\alpha S_1 + a_{12}\alpha S_2 - (c_{21}+r_{11}+r_{21}+\omega_3+h(t))C_1 + c_{12}C_2\)

\(\frac{dC_2}{dt} = a_{21}\alpha S_1 + a_{22}\alpha S_2 + c_{21}C_1 - (c_{12}+r_{12}+r_{22}+\omega_4+h(t))C_2\)

The acquisition rate \(\alpha\) in each equation: \(\alpha = \beta_1 C_1 + \beta_2 C_2\)

State-dependent removal (death or discharge) rates: \(\omega_i\)

State-independent, time-of-stay-dependent removal hazard: \(h(t)\)

Model components for \(R_0\) calculation

Transitions between, and out of, the \(S\) states in the absence of colonized patients: \[ S = \left( \begin{matrix} -s_{21}-\omega_1 & s_{12} \\ s_{21} & -s_{12}-\omega_2 \end{matrix}\right) \] Transitions between, and out of, the \(C\) states: \[ C = \left( \begin{matrix} -c_{21}-\sum r_{i1}-\omega_3 & c_{12} \\ c_{21} & -c_{12}-\sum r_{i2}-\omega_4 \end{matrix}\right) \] \(S\)-to-\(C\) state transitions when an acquisition occurs: \[ A = \left( \begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}\right) \]

Transmission rates from each colonized compartment: \[ \beta = \left( \begin{matrix} \beta_1 \\ \beta_2 \end{matrix}\right) \] Admission state distribution for susceptible compartments: \[ \theta = \left( \begin{matrix} \theta_1 \\ \theta_2 \end{matrix}\right) \] Moment-generating function associated with \(h(t)\): \[ M(t) \]

\(R_0\) formula

\[ R_0=\frac{\beta^T\int_0^\infty(1-F(t))e^{Ct}\int_0^te^{-C\tau}Ae^{S\tau}\theta\,d\tau\,dt}{1^T\int_0^\infty(1-F(t))e^{St}\theta\,dt} \]

• \(e^{St}\) and \(e^{Ct}\): matrix exponentials

• \(F(t)\): cdf of the distribution with \(h(t)\) as the hazard function: \(h(t)=\frac{F^\prime(t)}{1-F(t)}\)

  • If no state-dependent removal (\(\omega_i=0\)), this is the length of stay distribution

• Integrals of product of matrix exponential and \((1-F(t))\) can be expressed using the moment-generating function \(M(t)\) evaluated at the eigenvalues of the matrix

  • Formula can generate closed-form solutions when matrix exponentials can be expressed symbolically

  • We also derived numerical methods to calculate \(R_0\) by this formula - encoded in the facilityR0() function of our facilityepimath R package…

Example

\[ \frac{dS}{dt}=-(\alpha+h(t))S+\gamma C \] \[ \frac{dC}{dt}=\alpha S-(\delta_c+\gamma+h(t))C \]\[ \frac{dC_{cd}}{dt}=\delta_cC-h(t)C_{cd} \]

• \(S\): single susceptible state
• \(C\): colonized, not detected
• \(C_{cd}\): colonized, clinically detected
• \(\alpha=\beta(C+(1-\epsilon)C_{cd})\)
• \(\epsilon\): contact precaution effectiveness
• \(\gamma\): clearance of colonization rate
• \(\delta_c\): progression to clinical detection/infection rate

Matrix representation: \[ S = 0, \quad C = \left(\begin{matrix} -\delta_c-\gamma & 0 \\ \delta_c & 0 \end{matrix}\right), \quad A = \left(\begin{matrix} 1 \\ 0 \end{matrix}\right), \quad \theta=1, \quad \beta=\left(\begin{matrix} \beta \\ \beta(1-\epsilon)\end{matrix}\right) \]

Example \(R_0\) calculation in R

library(facilityepimath)
#Assumed values for carbapenem-resistant Enterobacterales (CRE):
gam <- 1 / 387       #clearance rate (per day)
eps <- 0.5         #contact precaution effectiveness

#Values fit to CRE data from long-term acute care hospitals (LTACH):
bet <- 0.0510      #transmission rate (per day)
deltaC <- 0.00845  #progression to clinical detection rate (per day)

S <- 0
C <- rbind(c(-deltaC - gam, 0), c(deltaC, 0))
A <- rbind(1, 0)
transm <- bet*c(1, 1 - eps)
initS <- 1

#Length of stay as a mixed gamma distribution (fit to LTACH):
px <- 0.580    #probability of following exponential distribution
rx <- 0.0285   #exponential distribution rate parameter
rg <- 0.179    #gamma distribution rate parameter
k <- 5.74      #gamma distribution shape parameter

mgf <- function(x, deriv = 0) MGFmixedgamma(x, prob = c(px, 1 - px), rate = c(rx, rg), shape = c(1, k), deriv)

facilityR0(S, C, A, transm, initS, mgf)

[1] 1.240727

Equilibrium calculation: facilityeq()

Same inputs as facilityR0(), plus recovery matrix \(R=(\gamma,0)\), and init vector includes \(C\) states (importation rate)

R <- cbind(gam, 0)
getEquilibriumFacilityPrevalence <- function(importationRate){
  init <- c(1 - importationRate, importationRate, 0)
  equilib <- facilityeq(S, C, A, R, transm, init, mgf)
  sum(equilib[2:3])
}

getEquilibriumFacilityPrevalence(0.2) #observed ~ 20% CRE importation at Chicago LTACHS

[1] 0.5002234

getEquilibriumFacilityPrevalence(0.1) #10% importation

[1] 0.3900065

getEquilibriumFacilityPrevalence(0.01) #1% importation

[1] 0.2266591

getEquilibriumFacilityPrevalence(0.001) #0.1% importation

[1] 0.1967067

Consequence of \(R_0>1\): high facility prevalence even when importation is infrequent

Intervention model: seek and decolonize

#Intervention parameters:
deltaS <- 1 / 28  #surveillance detect rate / day
gamD <- 1 / 10    #clearance rate / day from 
                  #decolonizing drug treatment

S <- rbind(c(0, 0), c(0, 0))
C <- rbind(c(-deltaS - deltaC - gam, 0, 0),
           c(deltaS, -deltaC - gamD, 0),
           c(deltaC, deltaC, 0))
A <- rbind(c(1, 0), c(0, 1 - eps), c(0, 0))
transm <- bet * c(1, 1 - eps, 1 - eps)
initS <- c(1, 0)

facilityR0(S, C, A, transm, initS, mgf)

[1] 0.8002351

\[ S = \left(\begin{matrix} 0 & 0 \\ 0 & 0\end{matrix}\right), \quad C = \left(\begin{matrix} -\delta_s-\delta_c-\gamma & 0 & 0 \\ \delta_s & -\delta_c-\gamma_d & 0 \\ \delta_c & \delta_c & 0\end{matrix}\right), \quad A = \left(\begin{matrix} 1 & 0 \\ 0 & 1-\epsilon \\ 0 & 0 \end{matrix}\right), \quad \theta=\left(\begin{matrix} 1 \\ 0 \end{matrix}\right), \quad \beta=\left(\begin{matrix} \beta \\ \beta(1-\epsilon) \\ \beta(1-\epsilon) \end{matrix}\right) \]

Intervention model: reduce length of stay

• LTACH length of stay (LOS) is a 58/42 mixture of patients following

  • Exponential distribution: mean 35 days, standard deviation 35 days
  • Gamma distribution: mean 32 days, standard deviation 13 days

• Mean LOS of ~34 days could be reduced while retaining LTACH status (>25 d)

• Does it matter (for \(R_0\)) which group(s) are targeted for LOS reduction?

  • Method 1: increase discharge rate of all patients (rx and rg)
  • Method 2: increase discharge rate of high-variance patients (rx)
  • Method 3: increase discharge rate of low-variance patients (rg)
  • Method 4: decrease fraction of high-variance patients (px)

  mgf <- function(x, deriv=0) MGFmixedgamma(x, prob = c(px, 1 - px), rate = c(rx, rg), shape = c(1, k), deriv)

Intervention model: reduce length of stay

Which LOS statistics correlate with \(R_0\)?

Summary

• facilityepimath calculates \(R_0\) and model equilibrium for facility models

  • Represent heterogeneous, dynamic susceptibility and transmissibility

    • Including intervention-induced states (surveillance, decolonization, etc.)

  • Represent realistic, complicated length of stay distributions

    • Details of distribution matter in unexpected ways

• Conclusions for CRE in long-term acute care hospitals

  • \(R_0>1\) for a single stay is a realistic possibility

  • Mid-stay surveillance for carriers can reduce \(R_0<1\)

  • Decolonization could greatly lower the surveillance resources needed

  • Length of stay efforts should emphasize both mean and variance reduction