equilibrium-and-R0

library(facilityepimath)

The functions in this package analyze compartmental differential equation models, where the compartments are states of individuals who are co-located in a healthcare facility, for example, patients in a hospital or patients/residents of long-term care facility.

Our example here is a hospital model with four compartments: two compartments \(C_1\) and \(C_2\) representing patients who are colonized with an infectious organism, and two compartments \(S_1\) and \(S_2\) representing patients who are susceptible to acquiring colonization with that same organism. The two colonized and susceptible states are distinguished by having potentially different infectivity to other patients and vulnerability to acquisition from other patients, respectively.

A system of differential equations with those 4 compartments may take the following general form:

\[ \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\) appearing in each equation, and governing the transition rates between the S compartments and the C compartments, is assumed to depend on the number of colonized patients in the facility, as follows:

\[ \alpha = \beta_1 C_1 + \beta_2 C_2 \]

Facility Basic Reproduction Number \(R_0\)

We will demonstrate how to calculate the basic reproduction number \(R_0\) of this system using the facilityR0 function. The following components of the system are required as inputs to the function call below.

The S matrix; pre-invasion susceptible state transitions

A matrix S governing the transitions between, and out of, the states \(S_1\) and \(S_2\) in the absence of any colonized patients “invading” the facility (i.e., all state transitions except for acquiring colonization):

\[S = \left( \begin{matrix} -s_{21}-\omega_1 & s_{12} \\ s_{21} & -s_{12}-\omega_2 \end{matrix}\right) \]

s21 <- 1; s12 <- 2
omega1 <- omega2 <- 0.1
S <- rbind(c(-s21-omega1, s12), c(s21, -s12-omega2))

The \(s_{ij}\) rates can be used to model transitions between different patient states that might alter susceptibility to acquiring the modeled organism, for example risky treatment procedures, drug exposures, or protective measures. The \(\omega_i\) rates model removal from the facility via discharge or death.

The C matrix; colonized state transitions

Next, we require a matrix C governing the transitions between, and out of, the states \(C_1\) and \(C_2\):

\[C = \left( \begin{matrix} -c_{21}-r_{11}-r_{21}-\omega_3 & c_{12} \\ c_{21} & -c_{12}-r_{12}-r_{22}-\omega_4 \end{matrix}\right) \]

c21 <- 0.1; c12 <- 0
r11 <- r22 <- 0.1; r12 <- r21 <- 0
omega3 <- omega4 <- 0.1
C <- rbind(c(-c21-r11-r21-omega3, c12), c(c21, -c12-r12-r22-omega4))

The \(c_{ij}\) rates can be used to model transitions between different colonized patient states that might be observable in data and/or alter the patient’s infectivity to other patients, for example clinical infection, detection status, or placement under protective measures. The \(r_{ij}\) rates govern clearance of colonization (thus, transitions back to one of the susceptible states), while the \(\omega_i\) rates model removal from the facility via discharge or death.

The A matrix; susceptible-to-colonized state transitions and relative susceptibility

The next input requirement is a matrix A describing the S-to-C state transitions when an acquisition occurs:

\[A = \left( \begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}\right) \]

a11 <- 1; a22 <- 2; a12 <- a21 <- 0
A <- rbind(c(a11,a12),c(a21,a22))

The \(a_ij\) rates describe which of the two colonized states can be entered from each of the two susceptible states at the moment of acquiring colonization; in the above example, state \(S_1\) patients move only to state \(C_1\) and state \(S_2\) patients move only to state \(C_2\). The \(a_ij\) values also describe the susceptibility of the two \(S\) states; in the above example, state \(S_2\) patients are twice as susceptible as state \(S_1\) patients.

The transm vector: transmission rates

The next required input is a vector transm containing the \(\beta\) coefficients (transmission rates from each colonized compartment) appearing in the \(\alpha\) equation above: \((\beta_1,\beta_2)\)

beta1 <- 0.2; beta2 <- 0.3
transm <- c(beta1, beta2)

The \(\beta_j\) values describe the transmissibility of each colonized state; in this example, state \(C_2\) patients are 50% more transmissible than state \(C_1\) patients. When the levels of colonization in the facility are \(C_1\) and \(C_2\), the acquisition rate of an \(S_i\) patient is \(a_{ii}(\beta_1C_1+\beta_2C_2)\)

The initS vector: susceptible state admission distribution

A vector initS containing the admission state probabilities for the susceptible compartments only (i.e., the pre-invasion system before a colonized patient is introduced): \((\theta_1,1-\theta_1)\)

theta1 <- 0.9; theta2 <- 1 - theta1
initS <- c(theta1, theta2)

The initS vector should sum to 1.

facilityR0 function call with \(h(t)=0\)

The time-of-stay-dependent removal rate \(h(t)\) is the remaining component in the master equations above. When \(h(t)=0\), the above components are sufficient to calculate the facility \(R_0\) as in the following example:

facilityR0(S,C,A,transm,initS)
#> [1] 1.763978

Note that when \(h(t)=0\), the removal rates (discharge and death) are entirely governed by the \(\omega_i\) values, which must be set such that patients in any state must be guaranteed to eventually reach a state for which \(\omega_i>0\); otherwise, patients can have an infinitely long length of stay in the facility.

The mgf function: moment generating function associated with \(h(t)\neq0\)

When the time-of-stay-dependent removal rate \(h(t)\neq0\), there is one additional argument required as input to the facilityR0 function: a function mgf(x,deriv) that is the moment-generating function (and its derivatives) of the distribution for which \(h(t)\) is the hazard function. This is the length of stay distribution when the state-dependent removal rates \(\omega\) are all zero, as in the following example.

omega1 <- omega2 <- omega3 <- omega4 <- 0
S <- rbind(c(-s21-omega1, s12), c(s21, -s12-omega2))
C <- rbind(c(-c21-r11-r21-omega3, c12), c(c21, -c12-r12-r22-omega4))

The length of stay distribution can be any statistical distribution with non-negative range, as long as the moment generating function (mgf) and its derivatives can be evaluated, as this calculation is required within the \(R_0\) formula. We currently provide three mgf functions in this package: one for the exponential distribution, MGFexponential, one for the gamma distribution MGFgamma, and one for a mixed gamma distribution MGFmixedgamma; the latter distribution can employ a weighted mixture of any number of different gamma distributions.

The mgf(x, deriv) function to be passed to facilityR0() must be defined as in the following example, which uses a gamma distribution.

mgf <- function(x, deriv=0) MGFgamma(x, rate=0.1, shape=3.1, deriv)

facilityR0 function call with \(h(t)\neq0\)

The call to facilityR0() is the same as above with the additional argument mgf (which defaulted to NULL when not provided):

facilityR0(S,C,A,transm,initS,mgf)
#> [1] 2.570984

Facility model equilibrium

We will demonstrate how to calculate the equilibrium of the full system of equations, with a given set of initial conditions, i.e. distribution of patient states at admission, using the facilityeq function. The following components of the system are required as inputs to the function call below.

The matrix S, the matrix C, the matrix A, the vector transm, and the (optional) function mgf are the same as those required by the facilityR0 function as described above. The following are new components:

The R matrix: recovery rates

The matrix R describing C-to-S state transition rates: \[R = \left( \begin{matrix} r_{11} & r_{12} \\ r_{21} & r_{22} \end{matrix}\right) \]

R <- rbind(c(r11,r12),c(r21,r22))

The init vector: admission state distribution

A vector init containing the admission state probabilities for all four compartments: \(((1-p_a)\theta_1,(1-p_a)(1-\theta_1),p_a\kappa_1,p_a(1-\kappa_1))\)

pa <- 0.05; kappa1 <- 1; kappa2 <- 1-kappa1
init <- c((1-pa)*theta1, (1-pa)*theta2, pa*kappa1, pa*kappa2)

Here we have modeled an importation probability (probability of an admitted patient being colonized at admission) of 5%, and assumed that all admitted colonized patients are in state \(C_1\).

The init vector should sum to 1.

facilityeq function call with \(h(t)\neq0\)

We can now calculate the equilibrium of the facility model:

facilityeq(S, C, A, R, transm, init, mgf)
#> [1] 0.2447589 0.1233193 0.1880693 0.4438525

The result is the portion of patients in the facility who are in each of the state \(S_1\), \(S_2\), \(C_1\), and \(C_2\), respectively, at equilibrium.

facilityeq function call with \(h(t)=0\)

We can also leave out the mgf argument when \(h(t)=0\), but should first reintroduce positive \(\omega\) values in the diagonal of the \(S\) and \(C\) matrices to represent discharge of patients:

omega1 <- omega2 <- omega3 <- omega4 <- 0.1
S <- rbind(c(-s21-omega1, s12), c(s21, -s12-omega2))
C <- rbind(c(-c21-r11-r21-omega3, c12), c(c21, -c12-r12-r22-omega4))
facilityeq(S, C, A, R, transm, init)
#> [1] 0.3685012 0.1739270 0.1650132 0.2925586