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
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.
S matrix; pre-invasion susceptible state
transitionsA 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) $$
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.
C matrix; colonized state transitionsNext, 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.
A matrix; susceptible-to-colonized state
transitions and relative susceptibilityThe 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) $$
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.
transm vector: transmission ratesThe next required input is a vector transm containing
the β coefficients
(transmission rates from each colonized compartment) appearing in the
α equation above: (β1, β2)
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)
initS vector: susceptible state admission
distributionA 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)
The initS vector should sum to 1.
facilityR0 function call with h(t) = 0The 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:
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.
mgf function: moment generating function associated
with h(t) ≠ 0When 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.
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:
R matrix: recovery ratesThe matrix R describing C-to-S state transition rates:
$$R = \left(
\begin{matrix}
r_{11} & r_{12} \\
r_{21} & r_{22}
\end{matrix}\right)
$$
init vector: admission state distributionA 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) ≠ 0We can now calculate the equilibrium of the facility model:
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) = 0We 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: