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₁ and C₂ representing patients who are colonized with an infectious organism, and two compartments S₁ and S₂ 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{d S_{1}}{d t} = - (s_{21} + (a_{11} + a_{21}) α + ω_{1} + h (t)) S_{1} + s_{12} S_{2} + r_{11} C_{1} + r_{12} C_{2}$

$\frac{d S_{2}}{d t} = s_{21} S_{1} - (s_{12} + (a_{12} + a_{22}) α + ω_{2} + h (t)) S_{2} + r_{21} C_{1} + r_{22} C_{2}$

$\frac{d C_{1}}{d t} = a_{11} α S_{1} + a_{12} α S_{2} - (c_{21} + r_{11} + r_{21} + ω_{3} + h (t)) C_{1} + c_{12} C_{2}$

$\frac{d C_{2}}{d t} = a_{21} α S_{1} + a_{22} α S_{2} + c_{21} C_{1} - (c_{12} + r_{12} + r_{22} + ω_{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:

α = β₁C₁ + β₂C₂

Facility Basic Reproduction Number R₀

We will demonstrate how to calculate the basic reproduction number R₀ 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₁ and S₂ in the absence of any colonized patients “invading” the facility (i.e., all state transitions except for acquiring colonization):

$S = (\begin{matrix} - s_{21} - ω_{1} & s_{12} \\ s_{21} & - s_{12} - ω_{2} \end{matrix})$

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 ω_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₁ and C₂:

$C = (\begin{matrix} - c_{21} - r_{11} - r_{21} - ω_{3} & c_{12} \\ c_{21} & - c_{12} - r_{12} - r_{22} - ω_{4} \end{matrix})$

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 ω_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 = (\begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix})$

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₁ patients move only to state C₁ and state S₂ patients move only to state C₂. The a_ij values also describe the susceptibility of the two S states; in the above example, state S₂ patients are twice as susceptible as state S₁ 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: (β₁, β₂)

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

The β_j values describe the transmissibility of each colonized state; in this example, state C₂ patients are 50% more transmissible than state C₁ patients. When the levels of colonization in the facility are C₁ and C₂, the acquisition rate of an S_i patient is a_ii(β₁C₁ + β₂C₂)

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 − θ₁)

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₀ 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 R₀ 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 = (\begin{matrix} r_{11} & r_{12} \\ r_{21} & r_{22} \end{matrix})$

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

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 S₁, S₂, C₁, and C₂, 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

- Facility Basic Reproduction Number R₀
- Facility model equilibrium

equilibrium-and-R0

Facility Basic Reproduction Number R0

The S matrix; pre-invasion susceptible state transitions

The C matrix; colonized state transitions

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

The transm vector: transmission rates

The initS vector: susceptible state admission distribution

facilityR0 function call with h(t) = 0

The mgf function: moment generating function associated with h(t) ≠ 0

facilityR0 function call with h(t) ≠ 0