3 Flow Between States

Here we describe the definition of the basic model of flows among states. Later chapters describe extensions to this basic model, both implemented (TODO) and unimplemented (TODO).

The per-capita rate of flow between any two states can be defined using the add_rate function.

sir = (sir
 %>% add_rate("S", "I", ~ (I) * (beta) * (1/N))
 %>% add_rate("I", "R", ~ (gamma))
)

The first add_rate function says that on each simulated day, individuals flow from the S box to the I box at a rate equal to the product of I, beta, and 1/N. This rate is a per-capita rate, and so the flow of individuals from S to I is S times the rate. Similarly, the second add_rate function says that I times gamma individuals flow from I to R every simulated day.

Here are the rules of these rate formulas (although these can be relaxed to some degree using techniques covered in later chapters).

  • both state variables and parameters can be referred to by name as variables in the formulas
  • variables must be encapsulated with parentheses
  • variables can be converted to their inverse (1/x) or complement (1-x)
  • variables, their inverses, and their complements can be combined using + and * operators

Examples of valid rate formulas for this model include:

  • (1-gamma) * (S) + (gamma)
  • (1/I)
  • (beta) * (gamma) + (beta) * (1/N) * (1-gamma)

Examples of invalid rate formulas include:

  • 1-gamma – no parentheses around this complement
  • beta * gamma – no parentheses around variables (but this can be addressed using struc objects – TODO: link to these)
  • ((S) + (R)) * (beta) – general grouping parentheses and common factors are not allowed (but this can be addressed using struc objects or intermediate computations – TODO: link to these)
  • (E) * (beta) – the variable E is in neither the parameter nor state variable

By specifying these rates, we now have a difference between the initial and final state vector.

state_init(sir)
##  S  I  R 
## 99  1  0
state_final(sir)
##           S           I           R 
##  0.01240091 12.18596876 87.80163033

The remainder of this chapter describes the theory behind these rates, and clarifies what they mean.

To start using this simple SIR model, please move on the the next chapter where you will learn how to do Simulation and explore epidemic trajectories.

3.1 State Flows

The numbers of individuals in each compartment is stored in the state vector, \(\mathbf{s}\), which contains one element for each compartment. At each time step, \(t\), the state vector is given by \(\mathbf{s}(t)\). The relationship between the state vector at successive time steps is given by the following.

\[ \mathbf{s}(t+1) = \mathbf{s}(t) + \mathbf{f}^{\text{in}}(t) - \mathbf{f}^{\text{out}}(t) \]

where \(\mathbf{f}^{\text{in}}\) and \(\mathbf{f}^{\text{out}}\) are the inflow and outflow vectors.

Continuing our concrete SIR model example, individuals flow from a box of susceptible individuals to infected to recovered individuals.

\[ S(t+1) = S(t) - \frac{\beta I(t)}{N}S(t) \] \[ I(t+1) = I(t) + \frac{\beta I(t)}{N}S(t) - \gamma I(t) \] \[ R(t+1) = R(t) + \gamma I(t) \]

where \(\beta\), \(\gamma\), and \(N\) are the transmission rate and population size parameters.

In this model the inflow to the \(S\) box is zero and the outflow from \(R\) is also zero. Therefore we can express this model in general terms as the following.

\[ \mathbf{s} = \begin{bmatrix} S \\ I \\ R \end{bmatrix} \] \[ \mathbf{f}^{\text{in}} = \begin{bmatrix} 0 \\ \frac{\beta I S}{N} \\ \gamma I \end{bmatrix} \] \[ \mathbf{f}^{\text{out}} = \begin{bmatrix} \frac{\beta I S}{N} \\ \gamma I \\ 0 \end{bmatrix} \]

3.2 Flow Matrix

Note that the outflow from \(S\) and inflow to \(I\) in the previous example has an identical magnitude. This is a common pattern in compartmental models. Many outflows are balanced perfectly by an associated inflow, to model individuals flowing from one compartment to another. McMasterPandemic assumes this balancing of flows to be the default situation, and therefore expresses both inflows and outflows in terms of an \(n\) by \(n\) flow matrix, \(\mathbf{F}\), that only requires specifying a single expression for each inflow-outflow pair. The element in the \(i\)th row and \(j\)th column of the flow matrix gives the flow from state \(i\) to state \(j\). The default inflow and outflow vectors can therefore be computed as the column sums and row sums respectively. Continuing the SI model example, we have the following flow matrix.

\[ \mathbf{F} = \begin{bmatrix} 0 & \frac{\beta I}{N}S & 0 \\ 0 & 0 & \gamma I \\ 0 & 0 & 0 \\ \end{bmatrix} \]

There are times however where one wants a particular flow from one state to another to include only the inflow component and not the outflow. For example in cases where death removes individuals from the population. We consider this and other examples of asymmetric flow in later chapters (TODO).

3.3 Rate Matrix

In modelling it is often more convenient to define per-capita rates of flow between compartments rather than total flow. The rate matrix, \(\mathbf{M}\), contains these per-capita rates. The elements, \(F_{ij}\), of the flow matrix can be computed from the rate matrix and the state vector as follows.

\[ F_{ij} = M_{ij}s_i \]

The rate matrix for the SI model is given by the following expression.

\[ \mathbf{M} = \begin{bmatrix} 0 & \frac{\beta I}{N} & 0 \\ 0 & 0 & \gamma \\ 0 & 0 & 0 \end{bmatrix} \]

The elements of the rate matrix are determined by expressions involving elements of the state vector and parameters, which we collect into a paramter vector, \(\mathbf{\theta}\). For the SI model \(\mathbf{\theta}\) contains two parameters.

\[ \mathbf{\theta} = \begin{bmatrix} \beta \\ \gamma \\ N \end{bmatrix} \]

3.4 Rate Matrix Dependence on State Variables and Parameters

Currently it is not possible to specify elements of the rate matrix as arbitrary arithmetic expressions involving state variables and parameters – although we plan to add this functionality. However, there is a reasonable degree of flexibility.

Each element of the rate matrix can be any expression that obeys the following rules. Any element, \(x\), of either the parameter or state vector can be used to define a factor in one of the following three forms.

  • Identity: \(x\)
  • Complement: \(1-x\)
  • Inverse: \(1/x\)

We collect these user-defined factors into a factor vector, \(\mathbf{y}\). Factors can be repeated in \(\mathbf{y}\) if required. Any number of factors can be multiplied together using * to produce a product. Any number of factors and products can be added together using +.

There is the following higher level nested structure associated with the factor vector, \(\mathbf{y}\).

  • All factors associated with the, \(i\)th, non-zero rate matrix element, \(M_{(i)}\), are grouped together in a contiguous block within \(\mathbf{y}\)
  • Within the \(i\)th block, all factors associated with the \(j\)th product (\(j = 1 ... n_i\)) in that block are grouped together in a contiguous sub-block
  • Within the \(i,j\)th sub-block, all factors are given an index, \(k = 1 ... m_{ij}\)

With these definitions, the dependence of any non-zero rate matrix element on the parameters and state variables is given by the following expression.

\[ M_{(i)} = \sum_{j=1}^{n_i} \prod_{k=1}^{m_{ij}} y_{ijk} \]

where \(y_{ijk}\) is the \(k\)th factor associated with the \(j\)th product associated with the \(i\)th non-zero rate matrix element.

from to n-fctrs n-prdcts n-vrbls state-dependent time-varying sum-dependent
S I 3 1 3 TRUE FALSE FALSE
I R 1 1 1 FALSE FALSE FALSE

McMasterPandemic is designed to be modular and allow multiple definitions of valid expressions for the rate matrix. In later chapters we describe this possibility (TODO).

3.5 Connections with Classic McMasterPandemic

params <- read_params("ICU1.csv")
classic_macpan = make_base_model(
  params = params,  state = NULL,
  start_date = "2021-09-10",
  end_date = "2021-10-10"
)
knitr::kable(rate_summary(classic_macpan, include_parse_info = FALSE))
from to
E_to_Ia E Ia
E_to_Ip E Ip
Ia_to_R Ia R
Ip_to_Im Ip Im
Ip_to_Is Ip Is
Im_to_R Im R
Is_to_H Is H
Is_to_ICUs Is ICUs
Is_to_ICUd Is ICUd
Is_to_D Is D
ICUs_to_H2 ICUs H2
ICUd_to_D ICUd D
H2_to_R H2 R
H_to_R H R
Is_to_X Is X
S_to_E S E

3.6 Topological Sort

The ordering of compartments with which individuals flow can be determined using the topological_sort function, provided that the graph is a DAG (i.e. there are no loops such as waning immunity). So for example the topologically sorted ordering of the classic McMasterPandemic model is given by the following.

topological_sort(classic_macpan)
##  [1] "S"    "V"    "E"    "Ia"   "Ip"   "Im"   "Is"   "H"    "ICUs" "ICUd"
## [11] "X"    "H2"   "D"    "R"

This can be used to order categorical variables for sorting tables and arranging plot facets.

knitr::kable(classic_macpan
 %>% rate_summary(include_parse_info = FALSE)
 %>% mutate(from = factor(from, topological_sort(classic_macpan)))
 %>% arrange(from)
)
from to
S_to_E S E
E_to_Ia E Ia
E_to_Ip E Ip
Ia_to_R Ia R
Ip_to_Im Ip Im
Ip_to_Is Ip Is
Im_to_R Im R
Is_to_H Is H
Is_to_ICUs Is ICUs
Is_to_ICUd Is ICUd
Is_to_D Is D
Is_to_X Is X
H_to_R H R
ICUs_to_H2 ICUs H2
ICUd_to_D ICUd D
H2_to_R H2 R 
model = classic_macpan
state_nms = base::setdiff(topological_sort(model), c("X", "V"))
nodes = data.frame(
  id = state_nms,
  label = state_nms,
  title = state_nms,
  physics = FALSE)
visNetwork(
  nodes,
  select(rate_summary(model), from, to) %>%
    mutate(arrows = "to")) %>%
  #visLegend() %>%
  #visLayout() %>% 
  visInteraction(zoomSpeed = 0.1) %>% 
  visPhysics(stabilization = TRUE) %>% 
  visHierarchicalLayout(direction = "LR")