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 C1 and C2 representing patients who are colonized with an infectious organism, and two compartments S1 and S2 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 α 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:

α = β1C1 + β2C2

Facility Basic Reproduction Number R0

We will demonstrate how to calculate the basic reproduction number R0 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 S1 and S2 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 sij 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 ω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 C1 and C2:

$$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 cij 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 rij rates govern clearance of colonization (thus, transitions back to one of the susceptible states), while the ω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 aij 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 S1 patients move only to state C1 and state S2 patients move only to state C2. The aij values also describe the susceptibility of the two S states; in the above example, state S2 patients are twice as susceptible as state S1 patients.

The transm vector: transmission rates

The next required input is a vector transm containing the β coefficients (transmission rates from each colonized compartment) appearing in the α equation above: (β1, β2)

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

The βj values describe the transmissibility of each colonized state; in this example, state C2 patients are 50% more transmissible than state C1 patients. When the levels of colonization in the facility are C1 and C2, the acquisition rate of an Si patient is aii(β1C1 + β2C2)

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): (θ1, 1 − θ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 R0 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 ωi values, which must be set such that patients in any state must be guaranteed to eventually reach a state for which ω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) ≠ 0

When the time-of-stay-dependent removal rate h(t) ≠ 0, 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 ω 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 R0 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) ≠ 0

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 − pa)θ1, (1 − pa)(1 − θ1), paκ1, pa(1 − κ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 C1.

The init vector should sum to 1.

facilityeq function call with h(t) ≠ 0

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 S1, S2, C1, and C2, 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 ω 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