Basic reproduction number of a two-group next-generation matrix as its dominant eigenvalue

Signature

T_K = k_11 + k_22; D_K = k_11 * k_22 - k_12 * k_21; R_0 = (T_K + sqrt(T_K^2 - 4 * D_K)) / 2

Function

Basic reproduction number function from transmission parameters or a next-generation matrix

Computational function

• Computational function for R_0 as the spectral radius of a next-generation matrix of any size

  Takes a next-generation matrix K with n groups, or the matrices F and V of a compartmental model from which K is F times the inverse of V, and returns R_0 as the spectral radius of K together with the dominant eigenvector scaled to sum to 1, the share of cases in each group once the generations settle. R and Python call an eigenvalue routine. Excel has no eigenvalue function, so it uses power iteration with helper columns, applying K to a vector of cases generation by generation, the growth factor per generation that Diekmann, Heesterbeek and Roberts take as the characterisation of R_0. The calculator line shows the two-group case over three generations started from one case in each group, beside the closed form of HE-FM-BRN-005, and the inputs therefore differ from the formula's variables, which cover two groups only.

  Inputs and outputs: K is the n by n next-generation matrix, entries zero or above, rows for the group of the new case and columns for the group of the infector, required unless F and V are given; F and V are the m by m new-infection and transition matrices at the disease-free equilibrium, passed as Fm and Vm in the code; k_11, k_12, k_21 and k_22 are the entries of a two-group K in the calculator line; a_1, a_2, b_1, b_2, c_1 and c_2 are the expected cases in each group in generations 1, 2 and 3; R_pow is the generation 3 total divided by the generation 2 total; R_two is the square root of the generation 3 total divided by the generation 1 total; T_K, D_K and R_0 are the trace, determinant and dominant eigenvalue of HE-FM-BRN-005; the R and Python versions return R0 and share, the dominant eigenvector scaled to sum to 1.

  Assumption: K is computed at the disease-free state and is non-negative, so its spectral radius is a real eigenvalue with a non-negative eigenvector. Power iteration on K settles when that eigenvalue is strictly larger in modulus than all others. A host and vector matrix with a zero diagonal has eigenvalues R_0 and minus R_0, so the one-generation ratio swings, and the two-generation ratio, or iteration on K plus the identity matrix as in the Excel version, is used instead. Three generations in the calculator line give an approximation, while the R and Python routines are exact to rounding.

  Worked example (Article's two-group matrix): Three generations from one case in each group give 2.8, 3.56 and 4.516 cases in total, so the one-generation ratio is about 1.2685 and the two-generation ratio about 1.270, against the exact R_0 of about 1.2685 (computed here for illustration). k_11 = 0.8; k_12 = 1.2; k_21 = 0.3; k_22 = 0.5; a_1 = 2; a_2 = 0.8; b_1 = 2.56; b_2 = 1; c_1 = 3.248; c_2 = 1.268; R_pow = 1.2685; R_two = 1.27; T_K = 1.3; D_K = 0.04; R_0 = 1.2685

  Worked example (Host and vector matrix with a zero diagonal): The one-generation ratio is 2.25 here and about 0.889 one generation later, while the two-generation ratio gives the exact 1.4142, the square root of 2 (computed here for illustration). k_11 = 0; k_12 = 4; k_21 = 0.5; k_22 = 0; a_1 = 4; a_2 = 0.5; b_1 = 2; b_2 = 2; c_1 = 8; c_2 = 1; R_pow = 2.25; R_two = 1.4142; T_K = 0; D_K = -2; R_0 = 1.4142

  Worked example (Three illustrative age groups, R and Python): For an illustrative matrix with rows (1.2, 0.4, 0.1), (0.6, 0.8, 0.2) and (0.1, 0.3, 0.4) the routines return R_0 of about 1.5726, with about 46.3%, 39.6% and 14.1% of cases in the three groups once the generations settle, and the Excel power iteration returns the same R_0 to six decimals after 100 helper columns. The calculator line does not apply to three groups. ngm_r0(matrix(c(1.2, 0.6, 0.1, 0.4, 0.8, 0.3, 0.1, 0.2, 0.4), 3)) returns R0 = 1.572561.

  Worked example (SEIR model built from F and V, R and Python): For van den Driessche's SEIR model with an illustrative beta S_0 of 0.5 per day, kappa of 0.2 per day, gamma of 0.25 per day and a natural death rate d of 1/25,550 per day, a 70-year lifespan, F has the single entry beta S_0 in row 1, column 2, and V has d + kappa and d + gamma on the diagonal with minus kappa in row 2, column 1. The routines return about 1.9993, the closed form kappa beta S_0/((d + kappa)(d + gamma)) of section 3.

  Excel: Put K in a range named NextGenMatrix and a column of n ones as the starting vector. Select the n cells of the next column and enter =(MMULT(NextGenMatrix,B2:B4)+B2:B4)/SUM(MMULT(NextGenMatrix,B2:B4)+B2:B4) as an array formula with Ctrl+Shift+Enter, shown for three groups with the start vector in B2:B4, then copy it across 100 columns and name the last column LastVector. =SUMPRODUCT(MMULT(NextGenMatrix,LastVector))/SUM(LastVector) returns R_0, and LastVector itself holds the case shares, since every column is scaled to sum to 1. Adding the vector to its image is iteration on K plus the identity, which settles even for a host and vector matrix, and only classic Excel functions are used.

  R: ngm_r0 <- function(K = NULL, Fm = NULL, Vm = NULL) { if (is.null(K)) K <- Fm %*% solve(Vm); ev <- eigen(K); i <- which.max(Mod(ev$values)); v <- Mod(ev$vectors[, i]); list(R0 = Mod(ev$values[i]), share = v/sum(v)) } Takes a matrix K, or Fm and Vm, and returns a list.

  Python: def ngm_r0(K=None, Fm=None, Vm=None): K = np.array(K, float) if Fm is None else np.array(Fm, float) @ np.linalg.inv(np.array(Vm, float)); w, vec = np.linalg.eig(K); i = int(np.argmax(np.abs(w))); v = np.abs(vec[:, i].real); return float(np.abs(w[i])), v/v.sum() Uses numpy imported as np and returns R_0 and the shares.

  Test (Excel power iteration agrees with the two-group closed form): For a two-group matrix the power-iteration result matches HE-FM-BRN-005 to six decimals. The expected result is TRUE. Excel check =ABS(SUMPRODUCT(MMULT(NextGenMatrix,LastVector))/SUM(LastVector)-(INDEX(NextGenMatrix,1,1)+INDEX(NextGenMatrix,2,2)+SQRT((INDEX(NextGenMatrix,1,1)-INDEX(NextGenMatrix,2,2))^2+4*INDEX(NextGenMatrix,1,2)*INDEX(NextGenMatrix,2,1)))/2)<1E-6

  Test (SEIR route matches the closed form): The R check abs(ngm_r0(Fm = matrix(c(0, 0, 0.5, 0), 2), Vm = matrix(c(0.2+1/25550, -0.2, 0, 0.25+1/25550), 2))$R0-0.2*0.5/((0.2+1/25550)*(0.25+1/25550))) < 1e-12 returns TRUE. A FALSE result shows the entry minus kappa left out of V, which returns 0 because no exposed case then moves on to the infectious stage.

  Common error (One-generation ratio on a host and vector matrix): Reading the ratio of successive generation totals as R_0 for the illustrative host and vector matrix gives 2.25 and then about 0.889, never settling, while R_0 is about 1.414; the two-generation ratio or iteration on K plus the identity removes the swing.

  Common error (Spectral bound of F minus V taken as R_0): The largest eigenvalue of F minus V, the Jacobian of the infected compartments, has the same sign as R_0 minus 1 by Theorem 1 of van den Driessche, but it is a growth rate. For the article's SIR infection it is 0.30 minus 0.10, or 0.2 per day, the early growth rate r, not R_0 of 3.

  Source: Diekmann O, Heesterbeek JA, Roberts MG. The construction of next-generation matrices for compartmental epidemic models. Journal of the Royal Society Interface. 2010;7(47):873-885. Section 1 on consecutive generations of infected individuals and the growth factor per generation as the characterisation of R_0, and section 3.2, R_0 = rho(K). van den Driessche P. Reproduction numbers of infectious disease models. Infectious Disease Modelling. 2017;2(3):288-303. Section 3, R_0 = rho(F V^-1), Theorem 1 and the SEIR example with F, V and R_0 = kappa beta S_0/((d + kappa)(d + gamma)).

  a_1 = k_11 + k_12; a_2 = k_21 + k_22; b_1 = k_11 * a_1 + k_12 * a_2; b_2 = k_21 * a_1 + k_22 * a_2; c_1 = k_11 * b_1 + k_12 * b_2; c_2 = k_21 * b_1 + k_22 * b_2; R_pow = (c_1 + c_2) / (b_1 + b_2); R_two = sqrt((c_1 + c_2) / (a_1 + a_2)); T_K = k_11 + k_22; D_K = k_11 * k_22 - k_12 * k_21; R_0 = (T_K + sqrt(T_K^2 - 4 * D_K)) / 2

Try this function

Implementations

• Excel

  Two-group dominant eigenvalue from a named matrix

  With the 2 by 2 next-generation matrix in a range named NextGenMatrix, columns for the infector's group and rows for the new case's group, the cell returns R_0, held in a cell named DominantEig. It uses (k_11 minus k_22)^2 plus 4 k_12 k_21 under the square root, which equals T_K^2 minus 4 D_K and avoids small negative values from rounding.

  =(INDEX(NextGenMatrix,1,1)+INDEX(NextGenMatrix,2,2)+SQRT((INDEX(NextGenMatrix,1,1)-INDEX(NextGenMatrix,2,2))^2+4*INDEX(NextGenMatrix,1,2)*INDEX(NextGenMatrix,2,1)))/2

Assumptions

• Next-generation matrix evaluated at the disease-free state with non-negative entries

  K is computed with everyone susceptible, so each entry is an expected count of new cases over the whole infectious period of one case and is zero or above; van den Driessche notes that F is entrywise non-negative and V a non-singular M-matrix, so the inverse of V and K are non-negative. Columns index the group of the infector and rows the group of the new case, as in the definition of Diekmann, Heesterbeek and Roberts. A transposed matrix has the same eigenvalues and so the same R_0, but its dominant eigenvector, the distribution of cases over groups, differs.

Worked examples

• Article's illustrative two-group next-generation matrix

  For the article's matrix with rows (0.8, 1.2) and (0.3, 0.5) the trace is 1.3 and the determinant 0.04, so R_0 is about 1.268. Neither diagonal entry reaches 1, yet transmission between the groups takes the system above the threshold.

  k_11 = 0.8; k_12 = 1.2; k_21 = 0.3; k_22 = 0.5; T_K = 1.3; D_K = 0.04; R_0 = 1.2685

• Rank-one two-group next-generation matrix whose R_0 equals the trace

  Rows (1.2, 0.6) and (0.8, 0.4), an illustrative matrix with proportional columns, have a determinant of zero, so R_0 equals the trace, 1.6, the case Diekmann, Heesterbeek and Roberts single out for 2 by 2 matrices.

  k_11 = 1.2; k_12 = 0.6; k_21 = 0.8; k_22 = 0.4; T_K = 1.6; D_K = 0; R_0 = 1.6

• Host and vector next-generation matrix with a zero diagonal

  With no direct transmission within either population, an illustrative 4 host infections per infected vector and 0.5 vector infections per infected host give R_0 of the square root of 2, about 1.414, the geometric mean form that van den Driessche describes for a vector-host model.

  k_11 = 0; k_12 = 4; k_21 = 0.5; k_22 = 0; T_K = 0; D_K = -2; R_0 = 1.4142

Common errors

• Largest diagonal entry of a next-generation matrix read as R_0

  For the article's matrix the larger diagonal entry, 0.8, suggests that spread dies out, while R_0 is about 1.268 because cases pass between the groups, so judging each group on its own transmission misses the threshold crossing.

• Square root dropped from a host and vector reproduction number or used for a host coverage target

  Van den Driessche notes that the literature often omits the square root for vector-host models, which keeps the threshold at 1 but changes the value, so the two forms cannot be compared across studies. In the illustrative host and vector matrix the product 2 would replace 1.414. Neither value gives a host coverage target through the threshold formula of HE-FM-BRN-003. Vaccinating a proportion p of hosts scales the host row of K by 1 minus p and gives R_0 of the square root of 2(1 minus p), which falls below 1 only when p exceeds 50%, while 1 minus 1/1.414 would suggest about 29.3% (computed here for illustration). Van den Driessche restricts herd immunity to control aimed at all individuals and uses type reproduction numbers for control aimed at one host type.

Sources

• Diekmann, Heesterbeek and Roberts on constructing the next-generation matrix

  Diekmann O, Heesterbeek JA, Roberts MG. The construction of next-generation matrices for compartmental epidemic models. Journal of the Royal Society Interface. 2010;7(47):873-885. Section 1, which defines element K_ij as the expected number of new cases with state-at-infection i generated by one individual just infected in state-at-infection j; section 2.1, equation 2.12, which obtains the dominant eigenvalue of a 2 by 2 matrix from its trace and determinant and notes that R_0 equals the trace when det K = 0; and section 3.2, which defines R_0 as the largest eigenvalue, R_0 = rho(K).

  View source →

• Diekmann, Heesterbeek and Metz on R_0 as a dominant eigenvalue

  Diekmann O, Heesterbeek JA, Metz JA. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology. 1990;28(4):365-382. The abstract, which defines the expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness in a completely susceptible population as the dominant eigenvalue of a positive linear operator.

  View source →

• van den Driessche on the next-generation matrix as F times the inverse of V

  van den Driessche P. Reproduction numbers of infectious disease models. Infectious Disease Modelling. 2017;2(3):288-303. Section 3, which splits the infected compartments into new infections F and other transitions V at the disease-free equilibrium, states that F is entrywise non-negative and V a non-singular M-matrix, that the (i, j) entry of F V^-1 is the expected number of secondary infections in compartment i produced by an individual introduced in compartment j and that R_0 = rho(F V^-1), with Theorem 1 linking R_0 below 1 to local stability; section 4.1, which states that herd immunity is applicable if a control strategy is aimed at all individuals and introduces type reproduction numbers for control aimed at one host population; and section 5, which gives the square root form for a vector-host model and notes that the square root is often omitted in the literature.

  View source →

• van den Driessche and Watmough on a precise definition of R_0 for compartmental models

  van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences. 2002;180:29-48. The abstract, which presents a precise definition of R_0 for a general compartmental model based on ordinary differential equations, shows that the disease-free equilibrium is locally asymptotically stable if R_0 is below 1 and unstable if it is above 1, and illustrates the definition with multigroup, staged progression and vector-host models.

  View source →

Canonical Identity