Next generation matrix reproduction number
Motivation
- General approach to calculating an \(R_0\) for a system of ODEs
- For example, given a system describing contact rates between population groups, what is the resulting \(R_0\)
- Need to think about which states contribute towards new infections, going back to the definition what does a typical infection look like when there are lots of infection classes?
- Want to understand how vaccination through reducing the susceptible population impacts the reproduction number (Mburu et al. 2021)
Derivation
See Van den Driessche and Watmough (2002) and Diekmann, Heesterbeek, and Roberts (2009).
\[ \mathbf{\dot{x}} = F(\mathbf{x}). \]
Linearise around DFE,
\[ \mathbf{\dot{x}} = (T - V)\mathbf{x}. \]
Derivation
\[ \mathbf{\dot{x}} = (T - V)\mathbf{x}. \]
Given \(\psi(0)\) as a typical infected individual at time \(0\), \(\psi(t)\) is the subsequent infected individuals at time \(t\) and,
\[ \dot{\mathbf{\psi}}(t) = -V\mathbf{\psi}(t) \]
The total number of infected individuals are \(\psi(\infty)\) and are,
\[ \psi(\infty) = \int_{t=0}^\infty \exp(-V t) \textrm{d}t\psi(0). \]
Derivation
\[ \psi(\infty) = \int_{t=0}^\infty \exp(-V t) \textrm{d}t\psi(0). \]
Implies
\[ \mathbf{\psi}(\infty) = V^{-1}\mathbf{\psi}(0). \]
Disease transition example 1
Consider the following system where \(I_1\) and \(I_2\) are infection classes,
\[ \begin{aligned} \dot{I_1} &= f_1 I_1 - \gamma_1 I_1, \\ \dot{I_2} &= f_2 I_2 - \gamma_2 I_2. \end{aligned} \]
\(f_1\) and \(f_2\) are forces of infection for the two states
\(\gamma_1\) and \(\gamma_2\) are the recovery rates
Disease transition example 1
\[ \begin{aligned} \dot{I_1} &= f_1 I_1 - \gamma_1 I_1, \\ \dot{I_2} &= f_2 I_2 - \gamma_2 I_2. \end{aligned} \]
\[ V = \begin{pmatrix} \gamma_1 & 0 \\ 0 & \gamma_2 \end{pmatrix} \]
The inverse of this matrix is,
\[ V^{-1} = \begin{pmatrix} 1/\gamma_1 & 0 \\ 0 & 1/\gamma_2 \end{pmatrix} \]
Disease transition example 1
\[ V^{-1} = \begin{pmatrix} 1/\gamma_1 & 0 \\ 0 & 1/\gamma_2 \end{pmatrix} \]
The time spent in each disease state is
\[\mathbf{\psi}(\infty) = V^{-1}\mathbf{\psi}(0).\]
\[\mathbf{\psi}(\infty) = V^{-1}\psi(t) = \begin{pmatrix} 1/\gamma_1 & 0 \\ 0 & 1/\gamma_2 \end{pmatrix}\begin{pmatrix} \psi_1(0) \\ \psi_2(0) \end{pmatrix}.\]
Disease transition example 1
\[\mathbf{\psi}(\infty) = V^{-1}\psi(t) = \begin{pmatrix} 1/\gamma_1 & 0 \\ 0 & 1/\gamma_2 \end{pmatrix}\begin{pmatrix} \psi_1(0) \\ \psi_2(0) \end{pmatrix}.\]
\[\mathbf{\psi}(\infty) = \begin{pmatrix} \psi_1(0)/\gamma_1 \\ \psi_2(0)/\gamma_2 \end{pmatrix}.\]
And so an individual in state \(1\) spends \(1/\gamma_1\) time in that state and an individual in state \(2\) spends \(1/\gamma_2\) time in that state
Disease transition example 2
Consider a different system,
\[ \begin{aligned} \dot{I_1} &= f_1 I_1 - \gamma_1 I_1, \\ \dot{I_2} &= \gamma_1 I_1 - \gamma_2 I_2. \end{aligned} \]
Here an individual once infected is in state \(1\), then transitions to state \(2\) before recovering.
\[ V = \begin{pmatrix} \gamma_1 & 0 \\ -\gamma_1 & \gamma_2 \end{pmatrix} \]
Disease transition example 2
\[ V = \begin{pmatrix} \gamma_1 & 0 \\ -\gamma_1 & \gamma_2 \end{pmatrix} \]
\[ V^{-1} = \frac{1}{\gamma_1\gamma_2}\begin{pmatrix} \gamma_2 & 0 \\ \gamma_1 & \gamma_1 \end{pmatrix} \]
Disease transition example 2
\[ V^{-1} = \frac{1}{\gamma_1\gamma_2}\begin{pmatrix} \gamma_2 & 0 \\ \gamma_1 & \gamma_1 \end{pmatrix} \]
\[ \begin{aligned} \psi(\infty) &= V^{-1}\psi(0) \\ & = \frac{1}{\gamma_1\gamma_2}\begin{pmatrix} \gamma_2 & 0 \\ \gamma_1 & \gamma_1 \end{pmatrix} \begin{pmatrix} \psi_1(0) \\ \psi_2(0) \end{pmatrix} \\ &= \begin{pmatrix} \psi_1(0)/\gamma_1 \\ \psi_1(0)/\gamma_2 + \psi_2(0)/\gamma_2 \end{pmatrix}. \end{aligned} \]
Disease transition for age-structured model
For age group indexed by \(a\) the infection class for an SIR-like model looks like,
\[ \dot{I}_a = f_aS_a - \gamma I_a. \]
\(f_a\) is the force of infection term which accounts for the interaction between other age groups.
Assumption that the recovery rate is the same for all age-groups. This greatly simplifies the analysis…
Disease transition for age-structured model
\[ \dot{I}_a = f_aS_a - \gamma I_a. \]
\[ V = \begin{pmatrix} \gamma & 0 & \ldots & 0 \\ 0 & \gamma & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \gamma \\ \end{pmatrix} = \gamma I \]
Disease transition for age-structured model
\[ V = \begin{pmatrix} \gamma & 0 & \ldots & 0 \\ 0 & \gamma & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \gamma \\ \end{pmatrix} = \gamma I \]
\[ \begin{aligned} \psi(\infty) &= V^{-1}\psi(0) = \gamma^{-1}\psi(0). \end{aligned} \]
i.e. regardless of which age-group you start in, the total recovery time is \(\gamma^{-1}\).
Back to the derivation
\[ \dot{x} = (F - V)x \]
\(F\) represents the transitions from a non-infected to a new infected class i.e. the new infections generated.
\[\int_{0}^\infty F\psi(t) dt\]
is the total number of new infections from the initial infections \(\psi(0)\).
Expected new infections
\[\int_{0}^\infty F\psi(t) dt\]
We can integrate this since \(F\) is constant in time,
\[F\int_{0}^\infty\psi(t) dt = FV^{-1}\psi(0).\]
\(FV^{-1}\) is the number of new infections in each compartment for an initial distribution of infections
Interpreting the entries of the matrix \(FV^{-1}\)
The \((j,k)\) entry of \(V^{-1}\) is the expected time an individual initially in compartment \(k\) spends in compartment \(j\)
The \((i,j)\) entry of \(F\) is the rate at which infected individuals in compartment \(j\) produce new infections in compartment \(i\)
So the \((i,k)\) entry of \(FV^{-1}\) is the expected number of new infections in compartment \(i\) produced by an infection in compartment \(k\).
Next generation matrix
\(FV^{-1}\) is therefore called the next generation matrix. It’s largest eigenvalue represents the largest number of new infections expected for an arrangement of initial infections and therefore,
\[R_0 = \rho(FV^{-1})\]
i.e. the R0 is the spectral radius of the next generation matrix (the largest eigenvalue)
Next generation matrix for an age-structured model
From Funk et al. (2019)
\[ \dot{I}_a = S_a \sum_b \beta_{ab} \frac{I_b}{N_b} - \gamma I_a, \]
where \(\beta_{ab} = p_{inf} \phi_{ab}\) is the infectivity rate, which is the probability of infection per contact multiplied by the contact rate \(\phi_{ab}\).
Next generation matrix for an age-structured model
\[ \dot{I}_a = S_a \sum_b \beta_{ab} \frac{I_b}{N_b} - \gamma I_a, \]
Linearise around the disease-free equilibrium,
\[ \dot{x}_a = (\sum_b \beta_{ab} \frac{N_a}{N_b} - \gamma) x_a, \]
Next generation matrix for an age-structured model
\[ \dot{x}_a = (\sum_b \beta_{ab} \frac{N_a}{N_b} - \gamma) x_a, \]
\(F_{ab} = \beta_{ab}\frac{N_a}{N_b} = p_{inf}\phi_{ab}\frac{N_a}{N_b}\) and \(V_{ab} = \gamma\delta_{ab}\).
Therefore the next-generation matrix (\(K\)) is,
\[K = (FV^{-1})_{ab} = \sum_{c} p_{inf}\phi_{ac}\gamma^{-1}\frac{N_a}{N_c}\delta_{cb} = p_{inf}\gamma^{-1}\phi_{ab}\frac{N_a}{N_b}\]
Adjusting for vaccination
Given an effective vaccinated proportion \(r_a\) in age-group \(a\), the modified next-generation matrix is,
\[K' = p_{inf}\gamma^{-1}\phi_{ab}\frac{N_a(1-r_a)}{N_b}\]
Defining the contact-adjusted immunity
In a well-mixed population the relationship between the effective \(R\) in a vaccinated population and \(R_0\) is,
\[R = (1 - r)R_0.\]
Where \(r\) is the proportion of effectively vaccinated individuals
Defining the contact-adjusted immunity
Where \(r\) is the proportion of effectively vaccinated individuals
we can re-arrange this to get,
\[r' = (1 - R/R_0) = \left(1 - \frac{\rho(K')}{\rho(K)}\right)\]
Note that \(p_{inf}\) and \(\gamma\) cancel out so this quantity is dependent only on the contact rate, population size, and proportion effectively vaccinated.
Contact matrices
Synthetic age-structured contact matrices can be constructed using the methods described in Prem, Cook, and Jit (2017) and Prem et al. (2021).
The original contact matrices were derived from the POLYMOD study which sampled a number of European countries (Mossong et al. 2008).
Method adjusts for population size and can be decomposed into settings: school, work, home, and all.
This can be easily generated using the conmat package
Contact matrices
This can be easily generated using the conmat package
Synthetic contact matrix for all
Synthetic contact matrix for work
Synthetic contact matrix for home
Synthetic contact matrix for school
Example vaccination distribution
Contact rate adjusted immunity example
library(epimixr)
# proportion immune in each age group
immunity <- sero_data$positive/100
# contact adjusted immunity
ca_immunity <- epimixr::adjust_immunity(synthetic_bc_pop_polymod[["all"]],
immunity = immunity)
# population immunity comparison
unadj_immunity <- sum(immunity * bc_pop$population)/sum(bc_pop$population)
print(scales::percent(ca_immunity))[1] "64%"[1] "76%"