Approximating reproduction numbers: a general numerical method for age-structured models

1.   Introduction

Reproduction numbers are key quantities in epidemiology, as they are usually related to concepts such as the occurrence of epidemic outbreaks, the herd immunity level, the final size, and the endemic equilibrium ^[1]. They were originally introduced in the context of demography and ecology, where they typically characterize either population persistence or extinction ^[2]. The most known and used example of a reproduction number in epidemiology is the basic reproduction number (or ratio) $R_0$, which describes the average number of secondary cases produced by a typical infected individual during its entire infectious period, in a completely susceptible population ^[3]. For deterministic models, $R_0 = 1$ is a threshold that determines the stability of the disease-free equilibrium: it is stable when $R_0 < 1$ and unstable when $R_0 > 1$. In simple models, $R_0$ also typically characterizes the herd immunity threshold (i.e., the proportion of population that should be immunized to prevent the spread of the disease) via the expression $1-1/R_0$. Variations of the basic reproduction number when the population is partially immune or when transmission is affected by control measures are also known as effective and control reproduction numbers ^[4].

For more complicated models (e.g., with heterogeneity), looking at $R_0$ alone may not be satisfactory for epidemic control ^[5]. For instance, this is true when it is possible to apply control measures only to a certain group of individuals (e.g., mosquitoes rather than humans in vector-borne infections) or against certain transmission routes (e.g., vertical rather than horizontal transmission). Under these circumstances, a simple and explicit relation between $R_0$ and the herd immunity level may no longer exist ^[5]. In this case, the type reproduction number, usually denoted by $T$, is a more appropriate measure, as it describes the expected number of secondary cases in individuals of a certain type produced by one infected individual of the same type, either directly or through chains of infection passing through any sequence of the other types ^[6]. As such, $T$ can be directly linked to the amount of control measures to be applied to one specific group of individuals to stop the spread. An extension of the type reproduction number, the state reproduction number, was proposed in ^[7], which allows for focus types not only states describing new infections but potentially any epidemic state (e.g., asymptomatic phase and symptomatic phase) [8, Section 5.2]. Finally, ^[9,10] considered a further generalization of these quantities, the target reproduction number, which focuses on specific interactions between types, rather than all interactions that involve a given type of individuals. It is important to underline that although different reproduction numbers have different biological interpretations, they typically share the same threshold for epidemic extinction with $R_0$ ^[10,11].

For deterministic models formulated as ordinary differential equations, a well-established and widely-used framework to compute $R_0$ is that of linearizing the model around the disease-free equilibrium and then computing the spectral radius of a Next Generation Matrix (NGM) ^[3,12]. The method is based on the idea of interpreting infection transmission as a demographic process, where a new infection is considered as a birth in the demographic sense. Intuitively, the NGM maps the distribution of infected individuals in one generation to the distribution of newly infected individuals in the next generation ^[3]; this is mathematically derived from the model's coefficients by suitably splitting the Jacobian into two parts, one accounting for the birth/infection processes and one accounting for the transition processes (which include, for instance, changes in the epidemiological state, death, or acquisition of immunity) ^[13]. The simplicity of the NGM method has considerably increased its popularity in the last decades. At the same time, the potentially arbitrary splitting into birth and transition processes described above has given rise to many criticisms and misconceptions about $R_0$ and its generational interpretation; see for example ^[14]. In fact, the interpretation of birth and transition is typically left to the modeller ^[15] and, while different choices agree on the sign of $R_0-1$, they can lead to different values of $R_0$. We refer to ^[16] for a recent didactic note about this issue.

In the context of age-structured models, which are often formulated as integro-partial differential equations with nonlocal boundary conditions, reproduction numbers are again associated with the linearization of the model around an equilibrium and the splitting of the linearization into birth and transition; however, in this case, the operators act on infinite-dimensional spaces. For example, it is well known that for a susceptible-infected-removed (SIR) model structured by demographic age without vertical transmission, $R_0$ can be computed as the spectral radius of a Next Generation Operator (NGO) ^[3], which is obtained by splitting the linearized operator in two parts–one accounting for all the infection processes and one accounting for all the remaining processes–which are both linear and defined on a subspace of the state space $L^1$, with values in $L^1$ itself [8, Section 6.2]. Thus, the procedure is very similar to that of the NGM method, but considering a space of $L^1$ functions rather than $\mathbb{R}^d$, for $d$ a positive integer. Working in $L^1$ is quite straightforward when processes described by boundary conditions are considered as transition processes. However, things get more involved when boundary conditions involve birth/infection processes (for instance, think at vertical transmission, or simple infection in the case of models structured by infection age). In this context, different strategies have been proposed in the literature. In ^[18], the authors introduced sequences of "approximating problems", for which the relevant reproduction numbers can be computed via the NGO method in the $L^1$ framework and such that these "approximating" reproduction numbers converge to the one of the original problem, see ^[19] for its applications. Alternatively, another possible approach is to work in extended spaces of the form $\mathbb{R}^d\times L^1$ and to develop the spectral theory following the results of ^[11]; see ^[20] for details on the extended space method and ^[8,21] for applications in this context.

Since reproduction numbers for age-structured models are defined as the spectral radius of operators acting between infinite-dimensional vector spaces, their analytical computation is typically difficult, unless one makes additional simplifying assumptions on the model coefficients (e.g., separable mixing). To overcome this problem, several numerical methods have been proposed to approximate $R_0$ by discretizing the birth and transition operators first to derive a finite-dimensional approximation of the NGO. The NGO is positive and, typically, compact, hence its spectral radius is a dominant eigenvalue (in the sense of largest in magnitude) ^[22] and the spectral radius of the discrete operator gives an approximation of $R_0$. ^[23,24] proposed a Theta method and a backward Euler method, respectively, to discretize the birth and the transition operators relevant to an age-structured epidemic model with no vertical transmission and a finite age span. Being based on finite-order methods, these two approaches guarantee, under suitable smoothness assumptions on the coefficients, a finite order of convergence. An improvement of these methods was proposed in ^[25,26,27] by using pseudospectral collocation (thus potentially guaranteeing an infinite order of convergence for smooth coefficients ^[28]) in the case of models with nontrivial boundary conditions and a finite age span. However, all these methods rely on the definition of $R_0$ in the $L^1$ framework discussed above, thus including the boundary condition in the domain of the transition operator and suffering from a lack of flexibility in the choice of the birth and transition processes.

In this paper, we introduce a general numerical method to approximate the reproduction numbers of a large class of multi-group age-structured models. We follow the idea of ^{[23,24,25,26,27]} by identifying a birth and a transition operator and discretizing them via pseudospectral collocation. To work with complete flexibility in the definition of the birth/infection and transition processes, we build our approach on the extended space framework by ^[8,21]. We define a suitable integral mapping from the extended space to the space of absolutely continuous functions, so that point evaluation is well defined and the birth processes included in the boundary conditions become part of the action of a new operator with a trivial boundary condition. The idea of going to the integrated framework has been previously successfully applied in ^{[29,30,31,32]}.

We assume that the maximum age is finite, which is equivalent to require that the survival probability (i.e., the probability of still being alive or infectious depending on the context) is zero after the maximum age. We focus on the applications of the method and we postpone the proof of convergence to a manuscript in preparation by the authors ^[33].

The paper is organized as follows. In Section 2, we consider a prototype linear multi-group age-structured model and, with reference to it, we illustrate the theoretical framework and the reformulation via integration of the state in the age variable. The numerical method is described in Section 3, alongside additional details about its implementation. In Section 4, we present some applications of the method to epidemic models taken from the literature. To illustrate the flexibility of the approach, we compute different types of reproduction numbers, depending on different interpretations of the age variable and the transmission term. To facilitate the reading, the modeling details and specific computations are collected in Appendix A. Finally, in Section 5 we discuss some concluding remarks.

2.   A general theoretical framework

In this section, we introduce a general, linear, multi-group, age-structured, population model, which encompasses many models of the literature typically obtained from the linearization of a nonlinear model around an equilibrium. With reference to this prototype model and following ^[8,21,34], we first recall the definition of the basic reproduction number and other relevant reproduction numbers useful to address some control strategies of infectious diseases within the extended space framework. Then, we introduce an integral mapping to the space of absolutely continuous functions and provide the equivalent definitions within the $AC$-framework, which is more advantageous for the development of the discretization technique presented in Section 3.

Let $a^{\dagger}\in (0, +\infty)$ denote the maximum age. We consider the following linear, multi-group, age-structured, population model:

where $x(t, \cdot) \in X: = L^1([0, a^\dagger], \mathbb{R}^d)$ for $t\ge 0$, $\beta \in L^\infty([0, a^\dagger]^2, \mathbb{R}^{d\times d})$, $b, \delta \in L^\infty([0, a^\dagger], \mathbb{R}^{d\times d})$, and $d$ is a positive integer. The $d \times d$ matrices $\beta(a, \xi)$ and $b(a)$ are assumed to be non-negative. The $d \times d$ matrix $\delta(a)$ has non-positive diagonal elements and all the off-diagonal elements are non-negative. Hence, $\delta(a)$ is an essentially non-negative matrix and the associated fundamental solution matrix is a non-negative, non-singular matrix [8, pag. 77].

In the context of infectious disease modeling, (2.1)–(2.2) enables one to describe several types of transmission routes and biological processes: the boundary term (2.2) can account for either vertical transmission (when $a$ represents demographic age) or for natural infection (when $a$ represents infection age), while the right-hand side of (2.1) includes horizontal transmission terms, as well as the removal of individuals via death or recovery.

To allow for flexibility in the definition of reproduction numbers associated to different ways of inflow into the infected compartments, we assume that $\beta$ and $b$ can be divided into the following:

where the non-negative matrices $\beta^{+}, \beta^{-}, b^{+}, b^{-}$ are chosen to conveniently split the inflow processes into two parts where $\beta^+, b^+$ collect the birth processes and $\beta^-, b^-$ collect the transition processes. The $+/-$ notation was inspired by ^[7].

We work in the extended state space of the density function and its boundary value, hence we consider the space $Z: = \mathbb{R}^d \times X$, equipped with the following norm:

where $|\cdot|$ is a norm in $\mathbb{R}^d$ and $\| \cdot \|_X$ is the standard $L^1$ norm, see ^[21] and [8, Chapter 6.4.2]. Furthermore, we introduce the subspace $Z_0 = \{ 0 \} \times X \subset Z$, where $0$ is now used to denote the null vector of $\mathbb{R}^d$.

We define the baseline transition operator $\mathcal{M}^Z\colon D(\mathcal{M}^Z)(\subset Z_0) \to Z$ as follows:

and the bounded linear birth operators $\mathcal{B}^Z_{\pm} \colon Z_0 \to Z$ such that

From the assumption on $\delta,$ it follows that $\mathcal{M}^Z$ is invertible, and then we can define the following:

To compute the reproduction number for the birth process described by $\mathcal{B}^Z_+$, we define the transition operator $\mathcal{M}_{-}^Z: = \mathcal{M}^Z-\mathcal{B}^Z_{-},$ whose domain coincides with $D(\mathcal{M}^Z).$ From $\mathcal{M}_{-}^Z = (\mathcal I-\mathcal{K}^Z_{-})\mathcal{M}^Z$, we have that if $\rho(\mathcal K_{-}^Z) < 1$, then $\mathcal{M}_{-}^Z$ is invertible [34, Section 4].

Then, the reproduction number $R$ for the birth process $\mathcal{B}^Z_{+}$ and the transition operator $\mathcal{M}_{-}^Z$ is the spectral radius of the positive operator

If (2.3) is a compact operator with positive spectral radius, then the reproduction number is a dominant real eigenvalue with an associated non-negative eigenfunction ^[22]. Note that, in line with the interpretation, the operator $(\mathcal I-\mathcal K^Z_-)^{-1} = \sum_{i = 0}^\infty (\mathcal K^Z_-)^i$ captures the chains of transmission through any sequence of any other type not included in $\mathcal{B}_+^Z$.

We generically refer to "reproduction number" to include several different interpretations as specific cases, including the basic reproduction number $R_0$ and the type reproduction number $T$, as well as more general definitions ^[6,8,35].

● If $\mathcal{B}^Z_+$ contains all processes leading to new infections, then (2.3) is the standard NGO and its spectral radius is precisely $R_0$.

● If $\mathcal{B}^Z_+$ only contains a subset of the new infections (e.g., horizontal vs vertical, vector vs host) and $\mathcal{B}^Z_-$ contains all other processes, then (2.3) is the type reproduction operator and its spectral radius is the type reproduction number, and $\mathcal{K}^Z_+$ and $\mathcal{K}^Z_-$ are the type-specific NGOs [8, pag. 477].

● If $\mathcal{B}^Z_+$ contains the inflow in a generic compartment (possibly not a state-at-infection), then the spectral radius of (2.3) is the state reproduction number according to the terminology in ^[7].

Moreover, depending on the assumptions on population immunity and interventions, (2.3) can capture the concepts of effective and control reproduction numbers. Additionally, (2.3) provides a unifying abstract framework to introduce the numerical method in Section 3.

Now, let us consider the Banach space $Y: = AC([0, a^\dagger], \mathbb{R}^d)$ equipped with the norm

and its closed subspace $Y_0 = \{\psi\in Y\ |\ \psi(0) = 0\}$. The integral operator $\mathcal V\colon Z\to Y$ given by

defines an isomorphism between $Z$ and $Y$ and between $Z_0$ and $Y_0$.

By defining $y(t, \cdot): = \mathcal V(0;x(t, \cdot))$, the model (2.1)–(2.2) for $x(t, \cdot) \in X$ is equivalent to the following multi-group model:

for $y(t, \cdot) \in Y_0$. Note that, in (2.4), we consider an absolutely continuous function as a measure. Now, we introduce the birth operators $\mathcal{B}^Y_{\pm} \colon Y_0\to Y$

and the transition operators $\mathcal{M}^Y, \mathcal{M}_{-}^Y:D(\mathcal M^Y) (\subset Y_0) \to Y$

where $D(\mathcal M^Y) : = \mathcal{V}(D(\mathcal M^Z)).$ Explicitly, they read as

and

It is clear that the spectral radius of $\mathcal{H}^Y: = \mathcal{V}\mathcal{H}^Z\mathcal{V}^{-1}$ coincides with $R$. In fact, the relations $\sigma(\mathcal H^Z) = \sigma(\mathcal H^Y)$ [36, Section 2.1] and $\sigma_p(\mathcal H^Z) = \sigma_p(\mathcal H^Y)$ [37, Proposition 4.1] hold, and there is a one-to-one correspondence of the eigenfunctions via the operator $\mathcal{V}$. Moreover, the compactness results for $\mathcal{H}^Z$ in $Z$ can be easily extended to the corresponding operators in $Y$ via the isomorphism $\mathcal{V}.$

3.   Numerical approximation via pseudospectral method

In this section, we illustrate how to approximate the reproduction number of a class of models that can be recast in the framework of Section 2. The idea is to derive a finite-dimensional approximation $\mathcal H_N^Y$ of $\mathcal H^Y$, and to approximate the eigenvalues of the latter through those of the former. This is achieved by separately discretizing the operators $\mathcal B_+^Y$ and $\mathcal M_-^Y$ in Section 2 via pseudospectral collocation ^[28,38].

In this section, we only work in the space $Y$ and, to simplify the notation, we drop the superscript $Y$ from the operators acting on $Y$. We adopt the MATLAB-like notation according to which elements of a column vector are separated by "; ", while elements of a row vector are separated by ", ".

Let us consider the space $Y_N$ of algebraic polynomials of degree at most equal to $N$, for $N$ a positive integer, on $[0, a^\dagger]$, taking values in $\mathbb{R}^d$, and its subspace $Y_{0, N}: = \{\psi_N\in Y_N\ |\ \psi_N(0) = 0\}$. For the Chebyshev zeros $\Theta_N: = \{a_1 < \dots < a_N\}$ in $(0, a^\dagger)$ ^[39], we introduce the restriction operator $\mathcal R_N\colon Y\to \mathbb{R}^{dN}$ defined as

and the prolongation operator $\mathcal P_{0, N}\colon \mathbb{R}^{dN}\to Y_{0, N}\subseteq Y_0$ defined as^*

^*Observe that $\Psi_i\in \mathbb{R}^d$ for every $i = 1, \dots, N$.

for $\{\ell_{0, i}\}_{i = 0}^N$ the Lagrange basis relevant to $\Theta_{0, N}: = \{a_0 = 0\}\cup \Theta_N$, i.e.,

Observe that $\mathcal R_N\mathcal P_{0, N} = \mathcal I_{ \mathbb{R}^{dN}}$ and that the composition

defines the Lagrange interpolation operator $\mathcal L_{0, N}\colon Y_0\to Y_{0, N}$ relevant to $\Theta_{0, N}$.

We derive the finite-dimensional approximations $\mathcal B_{N}\colon \mathbb{R}^{dN}\to \mathbb{R}^{dN}$ and $\mathcal M_{N}\colon \mathbb{R}^{dN}\to \mathbb{R}^{dN}$ of $\mathcal B_+$ and $\mathcal M_-$, respectively, as follows:

Then, the finite-dimensional approximation $\mathcal H_N\colon \mathbb{R}^{dN}\to \mathbb{R}^{dN}$ of $\mathcal H$ is obtained as $\mathcal H_N: = \mathcal B_N\mathcal M_N^{-1}$. Finally, we can use the eigenvalues of $\mathcal H_N$ to approximate those of $\mathcal H$. The discrete eigenvalue problem can be solved either by using the standard MATLAB function $\texttt{eig}$ or by solving the generalized eigenvalue problem $\mathcal B_{N} = \lambda \mathcal M_{N}$ ^[26]. Correspondingly, observe that the eigenvectors of $\mathcal H_N$ give an approximation of the values of the eigenfunctions of $\mathcal H$ at the Chebyshev zeros. Thus, an approximation of the eigenfunctions of $\mathcal H$ can be obtained by interpolating the eigenvectors of $\mathcal H_N$ at the nodes in $\Theta_{N}$ and, subsequently, an approximation of the eigenfunctions of $\mathcal H^Z$ can be obtained by applying $\mathcal V^{-1}$ to those polynomials.

3.1.   Implementation issues

Here, for the sake of simplicity, we restrict to the case $d = 1$, and we give an explicit description of the entries of the matrices representing the discretized operators. These follow from the following cardinal property of the Lagrange polynomials:

from which it is easy to see that the entries of the matrices $\mathcal B_{N}$ and $\mathcal M_{N}$ are explicitly given by

and

If the integrals in (3.1) and (3.2) cannot be analytically computed, then we need to approximate them via a quadrature formula. In this regard, for a function $\psi\colon [0, a^\dagger]\to \mathbb{R}^d$, we make the following approximation:

where $w_1, \dots, w_N$ are the Fejer's first rule-quadrature weights relevant to the Chebyshev zeros ^[40]. As for the integrals in $[0, a_i]$, $i = 1, \dots, N$, inspired by ^[41], we use the $i$-th row of the inverse of the following (reduced) differentiation matrix:

When dealing with models where the structuring variable $a$ lives in a "large" domain, it can be convenient to consider the change of variable $\alpha = a/a^\dagger$, where $\alpha\in[0, 1]$. Then, by defining $u(t, \alpha): = x(t, a^\dagger \alpha)$, (2.1)-(2.2) can be rewritten as follows:

Via integration, one can derive the corresponding equation for (2.4):

Then, the numerical approach can be easily adapted. Moreover, observe that one could also be interested in using different interpretations of "age" in the same model; for example see (4.3) which considers both infection age and disease age, or [42, Section 2.1], where an SIR model was proposed with susceptible individuals structured by the demographic age, infected individuals structured by the infection age, and removed individuals structured by the recovery age. In this context, different interpretations of the age variable can require one to work with different age intervals (or, equivalently, different maximum ages). The method can be easily extended to account for this case by resorting to the scaling described above (with appropriate modifications).

Finally, in the presence of breaking points (i.e., discontinuities in the model coefficients or in their derivatives), it is preferable to resort to a piecewise approach. More in the detail, given $0 = \bar a_0 < \bar a_1 < \dots < \bar a_M = a^\dagger$ the breaking points of the coefficients, for $M$ a positive integer, we approximate a function $\psi\in Y$ via a continuous function $\psi_N$ such that ${\psi_{N}}_{|_{[\bar a_{i-1}, \bar a_i]}}$ is a polynomial of degree at most $N$ on $[\bar a_{i-1}, \bar a_i]$ for every $i = 1, \dots, M$.

In this case, in order to simplify the implementation, a possible choice is to extend the mesh (of Chebyshev zeros plus the left endpoint) within each interval by adding the right endpoint. This still ensures convergence of interpolation under the choice made at the beginning of the section [43, Theorem 4.2.4]. Alternatively, it can be convenient to choose the Chebyshev extremal nodes as discretization points and the Clenshaw–Curtis quadrature formula to approximate the integrals in $[0, a^\dagger]$ ^[44,45]. The latter choice has been widely used in the literature for pseudospectral methods, and experimentally shows convergence properties comparable to those of Chebyshev zeros ^[45]. In the codes available at https://cdlab.uniud.it/software, we implement this latter choice.

4.   Age-structured epidemic models and reproduction numbers in applications

In this section, we introduce some examples of linear age-structured models in the context of infectious disease dynamics which are obtained from the linearization of a nonlinear model around an equilibrium. For each of them, we compute different types of reproduction numbers depending on different interpretations of the age variable and the transmission term. The examples are chosen to cover a range of cases: the first two models are somewhat simplified, which allows us to work with continuous rates and scalar equations and to have explicit expressions for the reproduction numbers; the third example is a system of equations that is useful to reflect on the interpretation of the birth processes; and finally, the last example involves a system of equations and application to real data, which requires piecewise constant parameters. All the reproduction numbers are computed using the method of Section 3 with $N = 100$ and the piecewise version of the numerical approach in the presence of breaking points. The modeling details and the linearization around the equilibria are described in the appendix, while Table 1 shows how the test examples fit into the general framework of Section 2. As we assume to work with a finite maximum age, we implicitly assume that the death/removal rates are infinite after the maximum age.

4.1.   An epidemic model structured by infection age

Let us consider the spread of an infectious disease in a closed population, with the infected individuals structured by infection age, in the presence of isolation measures upon detection of symptoms. We refer to [8, Section 5.3] and Appendix A.1 for further details about the nonlinear model, and to ^[46] for an application of an extended version of this model to study the impact of contact tracing on the containment of COVID-19.

Let $i(t, \tau)$ denote the density of infected individuals at time $t\geq 0$ and infection age $\tau\in [0, \tau^\dagger]$. The linearization around the disease-free equilibrium reads as follows:

The non-negative functions $b, \, \gamma \colon [0, \tau^\dagger] \to \mathbb{R}$ describe the per capita infection rate and recovery rate, respectively. In particular, we assume that $\gamma$ accounts for both natural recovery and isolation upon symptom onset, and we take the following:

where $\gamma_0, \epsilon, D$ are non-negative parameters whose interpretation is specified in Table 2.

In the absence of isolation from symptoms ($\epsilon = 0$), the basic reproduction number is $R_0 = {\int}_0^{\tau^\dagger} b(\tau) \mathrm{e}^{-\gamma_0 \tau} \mathop{}\!\mathrm{d} \tau$. In the presence of isolation ($\epsilon > 0$), the control reproduction number is as follows:

Note that, even though an explicit expression for $R_c$ is available in this case, its computation from the analytical formula requires numerical approximations.

In Figure 1, we investigated the impact of a delay from symptom onset to diagnosis ($D$) and the fraction of symptomatic infections ($\epsilon$) using parameter values inspired by the COVID-19 literature, collected in Table 2. $R_c$ increases linearly with the baseline transmission parameter $R_0$, and isolation of symptomatic individuals effectively reduces $R_c$ and promotes controllability (left panel). Moreover, the isolation is more effective if the proportion of symptomatic individuals is larger or the delay from symptoms to isolation is shorter (right).

Finally, in Figure 2, we illustrate the convergence behavior of the approximation error with respect to $R_0 = 1$ for increasing values of $N$ with $\epsilon = 0$.

4.2.   A multi-strain epidemic model with host age structure

We consider a model with two classes of infected individuals structured by demographic age, which describes the dynamics of two competing strains in a host population ^[49]. For applications of similar models to real infectious diseases, we refer to ^[50], where the case of influenza was discussed. In the model, susceptible individuals can be infected either with strain $1$ or with strain $2$, and enter the class of individuals infected with each strain. Cross-immunity is assumed, so that individuals recovered with any strain are immune towards both strains, and immunity is assumed to be permanent. In ^[49], it is shown that the coexistence of both strains in an endemic equilibrium is not possible when the parameters do not depend on age, but is possible for age-dependent parameters. The model derivation is described in detail in Appendix A.2, and the parameters are summarized in Table 3. We assume that the total population is at the (nontrivial) demographic steady state, and we neglect disease-induced mortality. The system can have four equilibria, of which three are boundary equilibria (one disease-free and two with only one strain present in the population) and one is an endemic coexistence equilibrium.

Let $i_j(t, a)$ denote the density of individuals infected with strain $j$, for $j = 1, 2$, at time $t\geq 0$ and demographic age $a \in [0, a^\dagger]$. The linearization around the disease-free equilibrium reads as follows:

with $\beta_j(a, \xi) = K_j(a) q_j(\xi) P^*(\xi)$, see also Table 3, and

The basic reproduction number $R_{0, j}$ for strain $j$ can be computed by individually considering each scalar equation in the absence of the other strain, and explicitly reads as follows:

Figure 3 shows the values of $R_{0, 1}$ and $R_{0, 2}$ varying the parameters $c_1$, $m_1$, and $c_2$, $m_2$, respectively, which are introduced for mathematical convenience to characterize the shape of the age-specific susceptibility of individuals (see also Table 3). The large values of the basic reproduction numbers for these parameter choices ensure the existence of the boundary equilibria where one strain is endemic in the population, and allow us to study the invasion reproduction numbers, as explained below.

Consider the boundary equilibria $E^*_1 = (s^*_1, i^*_1, 0)$ and $E^*_2 = (s^*_2, 0, i^*_2)$, where only strain 1 or strain 2, respectively, is present in the population (where $s^*_k(a)$ and $i^*_k(a)$ are the densities of susceptible and infected individuals at equilibrium divided by the stable age distribution $P^*(a)$). The equilibria $E^*_k$, $k = 1, 2$, do not admit an analytic expression, but their values can be solved numerically, as explained in Appendix A.2. The linearization at $E^*_k$ for $k = 1, 2$ and $j\neq k$ reads as follows:

The invasion reproduction number $R_k^j$, which describes whether strain $j$ can invade the equilibrium set by strain $k$, admits the following explicit expression:

Note that, in this case, computing the reproduction numbers from the analytical formula requires one to numerically approximate not only the integrals, but also the equilibria.

When both invasion reproduction numbers are larger than 1, the coexistence equilibrium is stable and the two strains can persist in the population ^[49].

Figure 4 shows the invasion reproduction numbers $R_2^1$ and $R_1^2$ when varying the parameters $m_1, m_2, c_1$, and $c_2$. As expected, each of these parameters has an opposite impact (either decreasing or increasing) on $R_2^1$ and $R_1^2$. The parameter $m_1$ is the only one that positively impacts the invasion of strain 1 (increasing $R_2^1$) and negatively impacts strain 2 (decreasing $R_1^2$).

In Figure 5, we plot the absolute errors of approximation for $R_{0, 1}$ and $R_{0, 2}$ for increasing $N$ with respect to the reference values computed from the analytical formulas. The numerical convergence is of infinite order, which is consistent with the fact that the parameters are of class $C^\infty$.

4.3.   An epidemic model with symptomatic and asymptomatic transmission

We consider the asymptomatic transmission model described in ^[7]. Upon infection, individuals enter an asymptomatic phase characterized by an infection age $\tau_1 \in [0, \tau_1^\dagger]$. Then, individuals can either recover without developing symptoms at a rate $\gamma_1(\tau_1)$, or they can develop symptoms at a rate $b_{21}(\tau_1)$, upon which they enter the symptomatic phase, which is characterized by a disease age (time since the onset of symptoms) $\tau_2 \in [0, \tau_2^\dagger]$, and from which they recover with a rate $\gamma_2(\tau_2)$.

Let $i_1(t, \tau_1)$ denote the density of asymptomatic individuals at time $t\ge0$ and infection age $\tau_1$, and let $i_2(t, \tau_2)$ denote the density of symptomatic individuals at time $t\ge0$ and disease age $\tau_2$. Assuming that the total susceptible population is normalized to one, the linearization around the disease-free equilibrium reads

where $b_{11}(\tau_1)$ is the per capita infection rate at the infection age $\tau_1$ in the asymptomatic phase, and $b_{12}(\tau_2)$ is the infection rate at the disease age $\tau_2$ in the symptomatic phase.

The parameter values used in the numerical simulations are listed in Table 4. Here, inspired by [15, Section 4] and ^[7], we can consider different definitions of "birth". According to the standard interpretation of $R_0$ as the number of new infections generated by one average infectious individual, we consider all new infections coming from asymptomatic and symptomatic individuals as birth processes, and consider the development of symptoms as transition process. On the other hand, since asymptomatic individuals are invisible from the point of view of the public health system (in the absence of other interventions like test-and-trace and asymptomatic testing), one could be interested in studying the effectiveness of control measures to the class of symptomatic individuals only (e.g., isolation upon symptoms) ^[7]. In this case, one could interpret the entrance in the class of symptomatic individuals as birth, while the asymptomatic phase is included in the transition operator, see also ^[51]. Therefore, we denote the state reproduction number of symptomatic individuals by $T_S$. As expected, the two reproduction numbers, $R_0$ and $T_S$, have different values in general, but the same threshold at 1, as seen in Figure 6 when varying $r: = b_{11}/b_{21}$. Additionally, Figure 6 illustrates another important theoretical property: the state reproduction number $T_S$ is finite (and well defined) only when the spectral radius of the NGO relevant to asymptomatic transmission is smaller than one. When the latter becomes larger than one, then asymptomatic individuals alone can sustain the epidemic, hence interventions that are targeted to only symptomatic individuals are not sufficient to control its spread. This feature is reflected in the behavior of $T_S$, which tends to infinity and becomes not well defined.

4.4.   A model for the spread of Rubella with vertical transmission

Rubella, also known as German measles or three-day measles, is an acute and contagious viral infection that can be vertically transmitted ^[52]. It is not particularly severe in children and adults, but infection during pregnancy can result in the so called congenital rubella syndrome, which can result in fetal death or congenital malformations in infants ^[53]. For women infected during early pregnancy (first trimester), the probability of passing the virus to the fetus is reported to be 90% ^[53] and, since there is no treatment for Rubella, the design of vaccination policies plays a fundamental role. In the last few decades, this has triggered a series of works by Anderson and colleagues, see for example ^[54,55,56].

Here, we consider a model inspired by ^[54] for the spread of congenital rubella syndrome in the United Kingdom. The model definition and the derivation of the disease-free equilibrium and the corresponding linearization are described in detail in Appendix A.3. The parameter definitions and values used in the numerical tests are collected in Table 5. Note that, in the literature, the term control reproduction number is typically used in the presence of interventions such as vaccinations. To simplify our terminology, here we refer to $R_0$ even in the presence of vaccinations.

Let $e(t, a)$ and $i(t, a)$ denote the density of exposed (not infectious) and infectious individuals, respectively, at time $t\geq 0$ and demographic age $a \in [0, a^\dagger]$, and let $s^*(a)$ denote the density of susceptible individuals at equilibrium, which is determined by the vaccination rate $\nu$ (see Appendix A.3 for the details). The linearization around the disease-free equilibrium reads as follows:

where $\beta(a, \xi) = s^*(a) \Pi(\xi) k(a, \xi)$ and $b(a) = qf(a)\Pi(a)$, see Table 5 for more details.

Following ^[54], we assume that the transmission rate $k$ is piecewise constant among six age groups, i.e., given $0 = \bar a_0 < \bar a_1 < \dots \bar a_6 = a^\dagger$,

so that we can estimate it from existing data using the well-known procedure of [56, Appendix A] (that we recall in Appendix A.3 for the reader's convenience). In this regard, we consider force of infection data from two different datasets: one for the South East of England in 1980 (case a) and one for Leeds in 1978 (case b), as summarized in Table 8. These datasets fix the age class division at 5, 10, 15, 20, and 30 years of age. The piecewise form in (4.5) gives us a Who Acquires Infection From Whom (WAIFW) matrix $(k_{ij})_{i, j = 1, \dots 6}$, which collects the contact rates between different age groups. We consider three different forms of the WAIFW, which capture different features in the transmission patterns. The estimated values of $k_{ij}$ are illustrated in Figure 7. More details on the parameters and data used for the estimation are available in Appendix A.3.

We compute the basic reproduction number, $R_0$, and the type reproduction number for horizontal transmission, $T_H$, for different choices of the vaccination rate $v$ and the WAIFW matrix. The results are collected in Table 6 (case a) and Table 7 (case b), rounded to four decimal digits.

The computed values of $R_0$ and $T_H$ are never substantially different (we can appreciate some differences only at the third decimal digit), thus suggesting that, for these data-informed parameter values, vertical transmission has a substantially smaller effect on the epidemic spread compared to horizontal transmission. Note that, for case b, the results obtained with WAIFW2 and WAIFW3 are identical: this is actually a consequence of the particular force of infection data in Table 8, which is the same for the two age groups 20–29 and 30–75 years old, see also [54, pag. 324]. Additionally, the values in Table 6 and Table 7 numerically illustrate the known relations between $T_H$ and $R_0$, namely $0 < T_H < R_0 < 1$ and $1 < R_0 < T_H$ ^[34].

Both $R_0$ and $T_H$ decrease for increasing $v$ for all choices of WAIFW matrix, see Figure 8. This is expected since a larger vaccination rate reduces the density of the susceptible population at equilibrium. On the other hand, different choices of the WAIFW matrix may result in slightly different quantitative values of $R_0$ and $T_H$, which reflects how different assumptions on the mixing patterns between individuals of different ages can affect transmission. This difference is even more evident looking at the eigenfunction of the type reproduction operator relevant to $T_H$ (normalized in the $L^1$-norm), see Figure 9.

Finally, we compute the type reproduction number for vertical transmission $T_V$. Figure 10 shows how $T_V$ depends on $\nu$ (left) and $q$ (right) for case b with WAIFW1. $T_V$ is a decreasing function of $\nu$ and linearly increases with $q$ (the ad hoc values of $\nu$ are chosen for illustrative purposes).

5.   Discussion and conclusions

In this paper, we have proposed a general numerical method to approximate the reproduction numbers of a class of age-structured population models with a finite age span. We presented applications to epidemic models that show how the method can compute different reproduction numbers, including the basic and the type reproduction numbers. Additionally, these examples show that, even when analytical expressions for the reproduction numbers are available, their computation may still require numerical approximations. Hence, our approach may represent an efficient and general alternative.

To our knowledge, this is the first numerical method that can approximate the many types of reproduction numbers for any splitting of the processes into birth and transition. This flexibility is made possible by working with the equivalent formulation for the age-integrated state, which has several advantages. First, it allows us to interpret processes described by the boundary condition as perturbations of an operator with trivial domain. Second, since the state is continuous, we can work with polynomial interpolation. Moreover, the additional regularity provided by the integral mapping permits us to investigate the convergence of the approximated eigenvalues without the need to look for a characteristic equation as in ^[27]. In fact, we can prove under mild (and biologically meaningful) assumptions that if $\mathcal M_-$ is invertible with bounded inverse, then there exists a positive integer $\bar N$ such that $\mathcal M_{N}$ is also invertible with bounded inverse for all $N\ge \bar N$, and that the eigenvalues of $\mathcal H_N$ converge to those of $\mathcal H$ with rate that depends on the regularity of the relevant eigenfunctions. In this paper, we experimentally investigated the convergence in some cases where the reproduction numbers were known, and we refer to the work in preparation ^[33] for the theoretical details.

Here, we focused on models with one structuring variable, namely age. However, to obtain a more realistic portrait of the dynamics of a population, one could also be interested in considering models with two (or more) structures (e.g., demographic age and infection age), see for example ^[57,58] and references therein. Unfortunately, considering an additional structuring variable brings in many difficulties in the theoretical and numerical study of these models. In fact, the hyperbolic nature of the infinitesimal generator of the semigroup makes the stability analysis more involved in the presence of discontinuities that propagate along the characteristic lines; for instance, these could be generated by corner singularities in the domain of the structuring variables, which is typically a rectangle in $\mathbb{R}^2$. In this context, working with the integrated state permits us to work with more regular spaces and to have an explicit expression for the infinitesimal generator, without the need for additional smoothness assumptions on the model coefficients, or compatibility assumptions on the boundary conditions as in ^[25]. A first work in this direction is ^[29], where the integrated state framework was used to approximate the spectrum of the infinitesimal generator relevant to the semigroup of a linear, scalar model with two structures. Following this idea, we plan to extend the method presented in this paper to models with more than one structure.

A limitation of this work is the assumption of finite age span. Hence, another interesting extension of this method regards models with unbounded structuring variables, which are common in the literature when considering probability distributions with unbounded support. For handling this problem numerically, one can either truncate the domain or resort to interpolation on exponentially weighted spaces and Laguerre-type nodes. For recent applications of these techniques to delay and renewal equations, which can also be used to model structured populations ^[59], see ^[60,61].

A.   Modeling details

In this appendix, we collect some further modeling and computational details regarding the models considered in Section 4, including the original formulation of the models in their nonlinear form. We use capital letters ($S$, $I$, $R \dots$) to denote the age-densities (or numbers depending on the context), and small letters to denote the densities divided by either the total size or the age distribution of the host population (e.g., $s(t, a) = S(t, a)/P(t, a)$).

A.1.   An epidemic model structured by infection age

Model (4.1) is obtained by partitioning the population into three classes (susceptibles, infected and removed), where only the infected class is structured by the infection age, see for example [62, Chapter 7] and [8, Section 5.3]. Let $S(t)$ and $R(t)$ denote the number of susceptible and removed individuals, respectively, at time $t\ge0$, and let $I(t, \tau)$ denote the density of infected individuals at time $t\ge0$ and infection age $\tau\in[0, \tau^\dagger]$. We assume that the total population $P: = S(t)+R(t)+{\int}_0^{\tau^\dagger}I(t, \tau) \mathop{}\!\mathrm{d} \tau$ is constant (no demographic turnover). The model reads as follows:

and can be rewritten in terms of the new variables

as

A.2.   A multi-strain epidemic model with host age structure

We consider a simplified version of the model proposed in ^[49], with two classes of infected individuals structured by demographic age, describing the dynamics of two competing strains in a host population. Let $S(t, a)$ denote the density of susceptible individuals at time $t\ge0$ and demographic age $a\in [0, a^\dagger]$, and let $I_1(t, a)$ and $I_2(t, a)$ denote the density of infectious individuals with strain $1$ and $2$, respectively, at time $t\ge0$ and demographic age $a\in [0, a^\dagger]$. We neglect additional disease-induced mortality. The full nonlinear model reads as follows:

where $Q(t) : = {\int}_0^{a^\dagger} P(t, \xi) \mathop{}\!\mathrm{d} \xi$, $P(t, a) : = S(t, a)+I_1(t, a)+I_2(t, a)$ is the age distribution of the host population, and the force of infection $\lambda_j$ satisfies

where $q_j$ is the age-specific infection rate for strain $j$, and $K_j$ is the age-specific susceptibility of susceptible individuals to strain $j$. For the demographic process, $P(t, a)$ is assumed to be at demographic equilibrium, i.e., the death rate $\mu$ and the per capita fertility rate $f$ are such that ${\int}_0^{a^\dagger} f(a) \Pi(a) \mathop{}\!\mathrm{d} a = 1$, where

is the survival probability. If $R_0^d > 1$, then there exists an endemic equilibrium such that the population profile at equilibrium, $P^*(a)$, satisfies

for $P_0 = {\int}_0^{a^\dagger} P^*(a) \mathop{}\!\mathrm{d} a = \Phi^{-1}(1/R_0^d)$.

To simplify the model, we define the variables

and obtain the following nonlinear system of equations:

where the force of infection can be written (equivalently) as follows:

The linearization around the disease-free equilibrium is reported in (4.2). Regarding the boundary equilibria $E_1^* = (s^*_1, i^*_1, 0)$ and $E_2^* = (s^*_2, 0, i^*_2)$, where only one strain is present in the population, they satisfy, for $k = 1, 2$,

with $s^*_k(0) = 1$, $i^*_k(0) = 0$, and $\lambda^*_j(a) = K_j(a) {\int}_0^{a^\dagger} q_j(\xi) P^*(\xi) i^*_j(\xi) \mathop{}\!\mathrm{d} \xi$, $j = 1, \, 2.$

Numerical solution of the endemic equilibrium. Note that system (A.3) for $s_k^*$ and $i_k^*$ cannot be solved analytically. To compute the equilibrium values to perform the numerical tests in the main text, the solutions were approximated numerically. Consider the equilibrium $E_1^*$ for illustrative purposes. For $n\in \mathbb{N}$, we discretize the interval $[0, a^\dagger]$ using Chebyshev extremal nodes $\{0 = a_0 < a_1 < \dots < a_n = a^\dagger \}$. Then, the derivative at each node is approximated using spectral differentiation, and the integral is approximated using Clenshaw–Curtis quadrature formulas with weights $w_{n, j}$. Let $S_n, I_n \in \mathbb{R}^{n+1}$ be two vectors such that, for $j = 0, \dots, n$, each entry $S_{n, j}$ approximates $s^*_1(a_j)$, and each entry $I_{n, j}$ approximates $i^*_1(a_j)$. Let $D \in \mathbb{R}^{(n+1)\times (n+1)}$ be the differentiation matrix associated with the nodes. Finally, let $K_n, Q_n \in \mathbb{R}^{n+1}$ such that $K_{n, j} = K_1(a_j)$, and $Q_{n, j} = w_{n, j} \, q_1(a_j) P^*(a_j)$. Then, we can write the following approximating system for the unknowns $S_n, I_n$:

where $*$ denotes the element-wise product. In practice, to facilitate the convergence to the nontrivial solution, we divide both terms by $(Q_n \cdot I_n)$ and solve the corresponding system.

A.3.   A model for the spread of Rubella with vertical transmission

We consider a model inspired by ^[54]. Let $M(t, a), S(t, a), E(t, a), I(t, a)$ and $Z(t, a)$ denote the density of individuals who are protected by maternal antibodies, susceptible, infected but not infectious, infectious, and immune (acquired naturally or via vaccination), respectively, at time $t\ge 0$ and demographic age $a\in [0, a^\dagger]$. The model reads as follows:

with the following boundary conditions:

where

is the force of infection, for $\hat k(a, \xi)$ the transmission rate between one individual of age $a$ and one individual of age $\xi$. We refer to Table 5 for the interpretation of the parameters. Note that an individual is assumed to have permanent immunity once infected.

Following [8, Chapter 6], we assume that ${\int}_{0}^{a^\dagger}f(a)\Pi(a) \mathop{}\!\mathrm{d} a = 1$, where $\Pi$ is the survival probability defined in (6.1), and that the age density of the host population $P(t, a): = M(t, a)+S(t, a)+E(t, a)+I(t, a)+Z(t, a)$ has already attained the stable age distribution $P(t, a) = P^*(a)$ defined in (6.2), for some $P_0 > 0$, see for example [8, Chapter 6]. For convenience, we define the standardized transmission rate as follows:

Then, if we consider new variables

we can reduce to the following model:

with the following boundary conditions:

where the force of infection can be expressed as follows:

The disease-free equilibrium $E^*: = (m^*(a), s^*(a), e^*(a), i^*(a), z^*(a))$ explicitly reads as follows:

where

Assuming a constant vaccination rate $v(a)\equiv v$ and $\eta > v$, the density of susceptible population reads as follows:

where

Observe that we have $s^*(a)\equiv 1$ for $v\equiv 0$. The linearization around the disease-free equilibrium is given in (4.4).

We estimate $k$ using real-world prevalence data taken from [54, Table 2] and their comments at page 324, which are collected in Table 8. To do this, we assume that $k$ is piecewise constant in the six age groups defined in Table 8, i.e.,

We assume three different structures for the WAIFW matrix $(k_{ij})_{i, j = 1, \dots 6}$, which are listed below:

In Table 9, we list the values of $k_i$, $i = 1, \dots, 6$, obtained from Table 8 using the procedure described below.

Estimation of age-dependent transmission rates. Here, we recall the procedure described in [56, Appendix A] to estimate the age-dependent transmission rate $k$ under the hypothesis that it is piecewise constant among different age groups, i.e.,

for given $0 = \bar a_0 < \bar a_1 < \dots < \bar a_n = a^\dagger$. We assume that the age-specific mortality rate $\mu$ has the following form:

so that ${\int}_0^{a^\dagger} \Pi(a) \mathop{}\!\mathrm{d} a = a^\dagger$, and that the age-specific force of infection

is known.

Then, the following algorithm can be applied:

(i) define $\psi_0: = 0$ and $\psi_i: = \sum_{j = 1}^i \lambda_j(\bar a_j-\bar a_{j-1})$ for $i = 1, \dots, n$;

(ii) define $\Psi_i: = \exp\left(-\psi_{i-1}\right)-\exp\left(-\psi_{i}\right)$ for $i = 1, \dots, n$;

(iii) solve the linear problem $\lambda_i = \frac{1}{\gamma} \sum_{j = 1}^n k_{ij}\Psi_j$ for $i = 1, \dots, n$.

Observe that the linear system in (iii) is over-determined; thus some hypotheses on the structure of the WAIFW matrix $(k_{ij})_{i, j = 1, \dots, n}$ are needed. We refer the reader to ^[56] for some possible choices.

Acknowledgments

We thank the reviewers for their valuable comments, which helped us to improve the presentation of the paper.

The authors are members of the INdAM Research group GNCS and of the UMI Research group "Modellistica socio-epidemiologica". The work of Simone De Reggi and Rossana Vermiglio was supported by the Italian Ministry of University and Research (MUR) through the PRIN 2020 project (No. 2020JLWP23) "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics", Unit of Udine (CUP: G25F22000430006). Francesca Scarabel acknowledges a Heilbronn Small Grant Call award by HIMR that partly funded a research visit during which some of this research was conducted.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Code availability

MATLAB demos are available at GitLab via https://cdlab.uniud.it/software.

Conflict of interest

Rossana Vermiglio is a special issue editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.