Next Generation Matrix Method in Epidemiology

• Next Generation Matrix Method is a matrix-based framework that quantifies key epidemiological thresholds and reproduction numbers in compartmental models.
• It employs an (F, V) decomposition to clearly separate infection transmission dynamics from recovery or removal, streamlining epidemic analysis.
• The method applies to both deterministic ODE and discrete Petri net models, extending to stochastic extinction analysis and final epidemic size estimation.

The next generation matrix (NGM) method is a rigorous, matrix-based approach for quantifying key epidemiological thresholds in compartmental epidemic and population dynamics models. Its primary function is to delineate the reproductive potential of a pathogen by mathematically encoding infection processes and transitions between states, providing the basic or gross reproduction number ($R_0$ or $R_g$), extinction probabilities, and, with suitable constructions, the final size of an epidemic. The NGM is central to both deterministic ordinary differential equation (ODE) models and, via extensions, to discrete-event frameworks such as Petri nets. Its scope encompasses the characterization of beginning and ending epidemic phases, supports a spectrum of model complexities, and connects closely to multitype branching process theory.

1. Mathematical Formulation of the Next Generation Matrix

The NGM method begins by expressing the system's dynamics near the disease-free equilibrium (DFE) in terms of a minimal set of "disease" or "infected" compartments. The governing equations are recast in the form: $$\frac{d\mathbf{i}}{dt} = \mathcal{F}(\mathbf{x}) - \mathcal{V}(\mathbf{x})$$ where $\mathcal{F}$ collects new infection inflow terms and $\mathcal{V}$ compiles all other transitions (natural removals, recoveries, inter-compartment flows).

Linearizing at the DFE, one computes the Jacobians $F = D_{\mathbf{i}} \mathcal{F}$ and $V = D_{\mathbf{i}} \mathcal{V}$, yielding matrices of size equal to the number of infected compartments. The next generation matrix is then: $K = F V^{-1}$ The dominant eigenvalue (spectral radius) $\rho(K)$ provides the basic reproduction number $R_0$, which quantifies the expected number of secondary cases produced by a typical infectious individual in a fully susceptible population (Avram et al., 2023, Reckell et al., 3 Jul 2025).

2. Construction Algorithms and Decomposition Uniqueness

Modern developments present algorithmic frameworks for the (F, V) decomposition (Avram et al., 2023). The user must specify an admissible set of disease variables, typically those that are zero at the DFE (the "zeroable set"). The procedure is as follows:

• Identify infectious compartments $\{\mathbf{i}\}$.
• For each equation governing $\mathbf{i}$, assign all terms corresponding to new infection inflows to $\mathcal{F}$, and all other terms to $\mathcal{V}$.
• Compute the Jacobians $F$ and $V$ with respect to $\mathbf{i}$, evaluated at the DFE.
• Adjust signs or positivity in $F$ to ensure $F$ contains only non-negative terms.
• $R_0$ is then $\rho(F V^{-1})$.

This process can be implemented algorithmically, for instance via Mathematica scripts that automate the decomposition and matrix construction (Avram et al., 2023). However, the decomposition is not unique; different splits between $F$ and $V$ are possible, leading to different but epidemiologically equivalent matrices under correct admissibility conditions. Notably, choosing a minimal sufficient set of infectious compartments often results in analytically tractable forms for $R_0$, while maximal choices may introduce unnecessary complexity.

3. Interpretation: Epidemic Phases and Thresholds

The NGM encapsulates two distinct roles, depending on its construction:

• Beginning Phase (Outbreak Threshold): In classical applications, the spectral radius of the NGM, often called the gross reproduction number $R_g$, defines the DFE stability. If $R_g < 1$, the DFE is locally asymptotically stable; if $R_g > 1$, sustained transmission (an epidemic) is possible. For multi-strain or structured models (e.g., drug-sensitive and drug-resistant tuberculosis), $R_g$ can be formulated as: $$R_g = \max\left\{R_1,\, R_2,\, \frac{R_3}{(1-R_1)(1-R_2)}\right\}$$ where $R_1$, $R_2$, and $R_3$ are strain- and pathway-specific reproduction numbers computed from the characteristic equation of $F V^{-1}$ constructed with a diagonal (or nearly diagonal) $V$ (Yang, 2020).
• Ending Phase (Final Size): A different construction, typically with a "full" (off-diagonal) $V$, yields an NGM whose spectral radius $p$ relates inversely to the asymptotic fraction of the population that remains susceptible at endemic equilibrium, $s^* = 1/p$ (Yang, 2020). The quadratic polynomial satisfied by $s^*$ is derived directly from the equilibrium equations: $$Pol(s) = R_{10}R_{20}\, s^2 - \left[R_{10}(1-R_{21}) + R_{20}(1-R_{11}) + R_{31}\right] s + (1-R_{11})(1-R_{21}) - R_{32} = 0$$ The smallest positive root corresponds biologically to $s^*$. For simple SIR/SEIR models, $s^* = 1/R_0$ emerges as a classical result.

4. Applications to Stochastic Extinction and Model Types

The NGM framework supports extension to stochastic settings via multitype branching process approximations. Extinction probabilities for an epidemic can be formulated as solutions to a fixed point equation (the "Bacaer equation") written in terms of $(F, V)$. For the Markovian class of models, the extinction probability vector $\mathbf{q}$ solves: $\Phi(\mathbf{q}) = 0$ where $\Phi$ reflects the branching process offspring distribution derived directly from the (F, V) decomposition (Avram et al., 2023). In SIR-type cases, this simplifies to a scalar polynomial equation whose solutions yield the probability of early extinction (fade-out) versus an epidemic.

For Petri net models, a class of directed bipartite graphs representing compartments and transitions, the NGM computation is entirely analogous. Places correspond to epidemiological compartments, and transitions (with weights) define token flows representing infection, progression, recovery, etc. The matrices $F$ and $V$ are constructed by evaluating the partial derivatives of transition functions at the DFE. The spectral radius of $F V^{-1}$ again yields $R_0$. This compatibility allows rigorous calculation of $R_0$ in both continuous ODE and discrete-event PN settings (Reckell et al., 3 Jul 2025).

Model framework	$F$ construction	$V$ construction	$R_0$ formula
ODE (deterministic)	Jacobian of new infection terms	Jacobian of transition/removal terms	$\rho(F V^{-1})$
Petri net (PN)	Derivatives of infection transitions	Derivatives of loss/transfer transitions	$\rho(F V^{-1})$

5. Algorithmic Implementations and Practical Considerations

Algorithmic approaches permit the automated generation of $(F, V)$ for a broad class of models. Notable points:

• The NGM method is implemented in available Mathematica scripts (e.g., NGM[mod, inf]), which split the model's Jacobian into positive "input" and complementary terms (Avram et al., 2023).
• For high-dimensional, multi-strain, or vector-host models, rational univariate representations (RUR) can reduce fixed-point systems for extinction probabilities or threshold calculations to univariate polynomials, facilitating analysis.
• While the (F, V) split can be performed in alternative ways for analytic simplification (by grouping more terms in $F$ or $V$), each approach may require parameter admissibility constraints to preserve mathematical structure (e.g., non-negativity).
• Numerical tools such as GPenSim can verify equivalence of ODE and PN NGM computations, confirming analytic results in simulation environments (Reckell et al., 3 Jul 2025).

6. Model Generality, Flexibility, and Theoretical Extensions

The NGM method is relevant for an extensive hierarchy of epidemic and population models:

• Single and Multi-strain Compartmental Models: Standard SIR/SEIR models, TB models with drug-sensitive and resistant strains, and multi-strain/multi-host configurations.
• Vector-Host Systems: The NGM captures square-root and maximum relationships for $R_0$ across transmission pathways (e.g., $R_0 = \max \{R_1, R_2\}$ in certain multi-strain scenarios).
• Nonlinear and Nonstandard Incidence Functions: The procedure accommodates polynomial, saturated, or otherwise nonlinear incidence, as demonstrated by direct computation of Jacobians at the DFE (Reckell et al., 3 Jul 2025).
• Discrete-event Systems (Petri Nets): Modular and visual modeling via Petri nets extends the reach of the NGM to systems with hybrid, discrete, or event-driven dynamics, maintaining equivalence with continuous ODE models.

The NGM's flexibility is illustrated by its ability to generalize classical results (e.g., $s^* = 1/R_0$) to models with treatment failure, superinfection, and complex transmission networks (Yang, 2020). The minimal disease set concept further enables tractable analysis in otherwise intractable high-dimensional models (Avram et al., 2023).

7. Limitations, Selection Criteria, and Theoretical Implications

While the method furnishes robust criteria for epidemic threshold and final size, several nuances are noted:

• The (F, V) decomposition is not unique; different admissible choices may yield equivalently valid threshold parameters, but analytic convenience and interpretability vary across decompositions (Avram et al., 2023).
• Maximal compartment selections for the infectious set always work but may result in less practical expressions; minimal sufficient sets yield simpler $R_0$ formulas that are advantageous for theoretical analysis and parameter inference.
• In certain models, extra conditions on parameters must be imposed to ensure non-negativity of the matrices and the epidemiological validity of $R_0$ (Avram et al., 2023).
• Comparison with Jacobian factorization methods shows that, in many cases, both methodologies provide identical threshold results, but the (F, V) approach offers a more general and systematic framework, particularly when automated.
• The link between the NGM and branching process theory suggests further theoretical insights into extinction and invasion phenomena, especially in the context of stochastic model generalizations.

---

The next generation matrix method is a foundational tool in modern quantitative epidemiology and mathematical biology, providing principled calculation of epidemic thresholds, extinction probabilities, and final-size relations across deterministic and stochastic frameworks. Its algorithmic and generalizable structure, validated by theoretical and computational examples, supports a broad spectrum of modeling applications, from classical ODE-based compartment models to discrete Petri net representations, enabling robust analysis even in the presence of model complexity, multiple strains, or intricate infection pathways (Yang, 2020, Avram et al., 2023, Reckell et al., 3 Jul 2025).