equilibrium-and-R0

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

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