NUMERICAL RECIPES FOR INVESTIGATING ENDEMIC EQUILIBRIA OF AGE-STRUCTURED SIR EPIDEMICS

. The subject of this paper is the analysis of the equibria of a SIR type epidemic model, which is taken as a case study among the wide family of dynamical systems of inﬁnite dimension. For this class of systems both the determination of the stationary solutions and the analysis of their local asymptotic stability are often unattainable theoretically, thus requiring the application of existing numerical tools and/or the development of new ones. Therefore, rather than devoting our attention to the SIR model’s features, its biological and physical interpretation or its theoretical mathematical analysis, the main purpose here is to discuss how to study its equilibria numerically, especially as far as their stability is concerned. To this end, we brieﬂy analyze the construction and solution of the system of nonlinear algebraic equations leading to the stationary solutions, and then concentrate on two numerical recipes for approximating the stability determining values known as the characteristic roots. An algorithm for the purpose is given in full detail. Two applications are presented and discussed in order to show the kind of results that can be obtained with these tools.


2676
DIMITRI BREDA, STEFANO MASET AND ROSSANA VERMIGLIO the many questions coming from the application side. The mutual interaction of these three aspects of mathematical modelling should always be fostered in order to counterbalance the inherent limits of each of them when considered individually. Part of the scope of this work is therefore to show how suitable numerical tools can help in advancing the theoretical analysis where it may appear to be stuck in light of considerable complexity of the model at hand.
As far as mathematics with time-dependent phenomena is concerned, it should be no surprise that the theory of dynamical systems today assumes a central position, a fact equally true for evolutions with discrete as well as continuous time. In the latter case, models driven by differential equations are dominant in many fields of sciences, and among the main concerns is the determination of their solution's behavior after a long time. We can expect to observe equilibrium solutions, but reality also rewards us many times with more complicated behaviors ranging from periodic solutions all the way to dynamics now commonly referred to as chaos. Which behavior will persist in a stable way is the main point: "what's the weather like tomorrow" then?
The subject of this work is the analysis of the equilibria of dynamical systems modeled by differential equations: their determination first and their stability second. The focus is on tools to take the path from the model to the key issue of asymptotic stability of the relevant equilibria. Of course, tools for this purpose are very well-established for Ordinary Differential Equations (ODEs) since the thesis of Lyapunov in 1892 (translated from Russian in [25]): one first computes the possible stationary solutions by solving a system of nonlinear algebraic equations, second linearizes the original model around the equilibrium whose stability has to be determined, and third extracts the eigenvalues of the matrix of coefficients of the resulting linear differential system. The stability is then determined from the real part of these eigenvalues.
In this paper we focus the attention on systems of infinite dimension, arising from Partial Differential Equations (PDEs), Functional Differential Equations (FDEs) or other abstract settings. Although the theoretical recipe for the analysis of the local stability remains unaltered in its substance, its application is accompanied by such an intrinsic amount of complexity that resorting to numerical tools is the rule rather than the exception. As far as the solution of the system of nonlinear algebraic equations for the equilibria is concerned, nonlinear solvers are the standard tools, whether or not employed in a continuation framework [1,18]. The crucial issue emerges when the stability of an equilibrium has to be inferred after the linearization of the model: this is truly a problem of infinite dimension and numerics is often mandatory.
Recently, a family of techniques based on approximating the spectrum of the infinitesimal generator of the semigroup of solution operators associated with the linearized model has been proposed and fully developed in the context of the stability analysis of equilibria of dynamical systems from several classes of differential equations; these include retarded, neutral-retarded and mixed-type FDEs [8,9], PDEs of reaction-diffusion type with delay in the reaction term [10], as well as diverse models in the context of population dynamics devoted to the analysis of age-structured populations and various types of diseases [3,6,7,12,13,17,24]. As it is our intention to concentrate on the numerical tools rather than on the analysis of the model features and their interpretation, we consider a prototypical but still very general case study. More specifically, we focus here on systems of nonlinear hyperbolic PDEs describing epidemic models of SIR type (see e.g. [2,26]), where the mechanisms of transmission of the disease depend both on the age of the individuals as well as on some weighted measure of the demographic population. In particular, devoting our attention to the numerical techniques, we also extend the analysis to the case of systems, rather than scalar equations as in [7], and moreover allow for the treatment of general discontinuities in the age-dependency of the model parameters. The algorithm for constructing the approximation matrix from whose eigenvalues the stability information can be recovered is also detailed, as to permit a direct implementation by potential users.
The prototype model is introduced in its generality in Section 2. The search for equilibria through the construction of a nonlinear algebraic system is addressed in Section 3. The linearized stability analysis is tackled in Section 4 through two common approaches, namely the characteristic equation in Section 4.1 and the spectrum of the infinitesimal generator in Section 4.2. As for the latter approach, which we choose to follow, the model is further generalized and the relevant numerical discretization is presented in Section 5. The numerical recipes given in the paper are finally applied in Section 6 to two different case studies taken from [13,17], where one may possibly appreciate the type of analysis these tools are able to support and the results that can be provided. Some concluding remarks and future perspectives are given in Section 7.
2. The model. Let us consider a population described by its age density N (a, t) at time t ≥ 0 with age a ∈ [0, a † ], where a † is the maximum age. It is assumed that the growth of the population obeys a model of Gurtin-MacCamy type [20], i.e.
is the demographic net reproduction rate and is the size of the population with weight kernel q. Furthermore, and µ and β are the intrinsic mortality and fertility rates, respectively. The function Φ describes the density-dependence of births; it is assumed to be Lipschitzcontinuous, non-increasing from Φ(0) = 1, decaying as lim x→+∞ Φ(x) = 0 and decreasing on [0, x 0 ] for some x 0 ∈ R such that R d 0 Φ(x 0 ) < 1. By assuming the normalizing condition a † 0 β(a)π(a)da = 1 (2) for the survival probability π(a) := exp(− a 0 µ(σ)dσ), it can be shown that the demographic equilibrium N * (a) = N * (0)π(a) exists with N * (0 loc (0, a † ) and q ∈ L ∞ (0, a † ) are assumed to be nonnegative in [0, a † ] together with a † 0 µ(a) da = +∞ since no individual may live past the maximum age, i.e. π(a † ) = 0.
Let us now assume that the population develops an SIR type epidemics [2]. Therefore, we set N (a, t) = S(a, t) + I(a, t) + R(a, t), where S(a, t), I(a, t), R(a, t) denote the age densities at time t of the susceptible, infective and removed individuals, respectively. Susceptibles may become infected at rate λ(a, t), the force of infection, and infected become removed at rate γ(a). Considering also natural deaths at rate µ(a) and an extra mortality α(a) due to the disease, the equations satisfied by S, I and R are where, neglecting transient maternal immunity, it is assumed that all newborns are susceptible and equally produced by all kinds of individuals. Observe that in the presence of the disease we have that ∂N (a,t) ∂t + ∂N (a,t) ∂a = −µ(a)N (a, t) − α(a)I(a, t), which shows how the disease-induced mortality α affects the population growth.
For the parameters regulating the infection mechanism we assume α(a) ≥ 0 and a constitutive rule for the force of infection λ(a, t) assigned for the moment under the standard assumption where F (a, σ) is the contact coefficient between susceptible individuals aged a and the infected aged σ. We will specify a more detailed behavior in the sequel. Existence and uniqueness of a solution of (3) can be proved by following the standard techniques presented e.g. in [21,27]. In this work, we focus on the determination of the nontrivial equilibria and on the analysis of their asymptotic stability entirely through numerical rather than analytical tools.
3. Stationary solutions. Searching for steady states of (3) leads to the system of ODEs where and λ * (a) = In order to gain some insight into the problem we need to specify some extra assumptions for the force of infection (4). To preserve generality as much as possible, and to encompass different epidemic models such as those treated in e.g. [13,17], let us set Namely, we consider especially contacts that take the form of WAIFW matrices ("Who Acquires Infection From Whom", [2]), i.e. susceptibles are divided into n age groups I l = (a l−1 , a l ), l = 1, . . . , n, with a 0 = 0 and a n = a † and, similarly, infected are divided into m possibly different age groups J r = (b r−1 , b r ), r = 1, . . . , m, with b 0 = 0 and b m = a † . In (8) above χ I denotes the characteristic function of the interval I, f lr are the contacts coefficients between susceptible and infected groups and, finally, g is a weight kernel for the latter. Substituting where W r (t) := Jr g(σ)I(σ, t)dσ, for which we set W (t) := (W 1 (t), . . . , W m (t)) T ∈ R m . Similarly, for the stationary solution(s) we set λ * (a) = n l=1 and W * := (W * 1 , . . . , W * m ) T ∈ R m . Direct integration of (5) leads to the stationary solution components where we set, for brevity, η * (a, Through (12) the search for equilibria is reduced to the solution of the system of nonlinear algebraic equations for the unknown vector (S * (0), W * T ) T ∈ R m+1 . The first equation comes from the initial condition for S * in (5), the other equations follow from (11).
It is clear that the trivial equilibrium q(a)S * (0)π(a)da = 1 and having observed that η * (a, 0) = 1 and ν * (a, 0) = 0. Excluding these immediate solutions and focusing on the nontrivial ones, i.e. the so-called endemic equilibria E e = (S * (0), W * T ) T for which W * ∈ R m \ {0}, system (13) is far from amenable to analytical solution even though very special and more specific behaviors were assumed for the model parameters and rates. Therefore, resorting to numerics represents the only way forward, and a number of standard and well-developed approaches can be adopted. A first is based on continuation, each step of which relies on Newton like iterations (see e.g. [23] and software packages such as AUTO, CONTENT, DSTOOL, MATCONT and XPPAUT). A second is based on selecting a nonlinear solver (either Newton or more sophisticated ones, e.g. fsolve in MATLAB) and trying to inspect the solution hyperplane by running the tool for a sufficiently fine mesh of nodes chosen as independent starting points. Both have their own pros and cons. For example, on the one hand, continuation needs a starting point for each equilibrium branch, on the other hand, the selection of initial guesses for the second procedure may soon result in difficulties unless the degree of complexity of the system is quite low. In Section 6 we will adopt both alternatives in two different case studies.
z(a, t) := R(a, t) − R * (a) to derive the linear system where the initial functions are given by and is the perturbation for the total population. The principle of linearized stability (see e.g. [21] in this context) tells us that the local asymptotic stability of the concerned equilibrium can be inferred from the stability of the zero solution of (14). The latter is given by the dominant exponential rate for the perturbations x, y and z, which in turn is deduced by analyzing the position in C of the so-called characteristic roots. Among others, two standard techniques are available to determine the characteristic roots theoretically, which we discuss below. The difference between them lies in the way their values are formulated and how, consequently, they can be determined. Nevertheless it ought to be stressed that there are infinitely many characteristic roots. This fact follows from the infinite-dimensional nature of the underlying dynamical system, which can be easily appreciated by considering the evolution of a general solution u(a, t) as the time evolution of the age profile u t in a suitable state space, where u t (a) := u(a, t) and typically is assumed to leave in L 1 for models of populations. Given the above, it is not surprising that (exceptional cases excluded) we need to resort again to numerical methods. The choice of one alternative or the other plays an essential role in this context.

The characteristic equation.
For population models such as those considered in this work it is possible to write a so-called characteristic equation. For this common technique (and the theory behind it) the reader is referred to e.g. [21]. As far as the derivation of a characteristic equation is concerned, a detailed procedure is developed in [13]. Here we recall the fundamental steps, but avoid a complete treatment for the considerations stated below.
Consider system (14), which is a set of hyperbolic linear PDEs with left inflow and bottom initial function for each component. Each equation can be integrated along the characteristics, which read a(t) = t + constant. The obtained solutions can then be inserted into the definition of x(0, t) in (14), as well as into each of the equations in (15). This can be verified to lead to a set of m + 1 integral equations of the form . The final step consists now of taking the Laplace transform, or less formally to substitute exponential functions for the unknown solutions, i.e. x(0, t) = e zt v 0 and w l (t) = e zt v l , l = 1, . . . , m, for z ∈ C. The procedure leads to nontrivial solutions where the entries of the matrixĈ areĉ lr (z) = a † 0 e −za c lr (a)da, l, r = 0, . . . , m.
Equation (16) is the characteristic equation and its solutions z ∈ C are the characteristic roots. There are infinitely many due to its transcendental nature. Moreover, since the entries ofĈ are entire functions, they form a discrete set in the complex plane, and any right half-plane contains only a finite number of them. From a theoretical point of view, the determination of the characteristic roots as solutions of the characteristic equation is rather limiting, in spite of the fact that it is one of the preferred tools adopted by mathematical analysts and applied scientists. Its use is, in fact, restricted to very simple cases or lucky instances demanding too often strong assumptions on the model parameters. A natural idea could then be that of exploiting the characteristic equation in combination with numerical solvers for nonlinear equations. This requires the above derivation of the equation coefficients anyhow and furthermore, this derivation is an unavoidable source of roundoff perturbations even when it is assumed to start from exact model data. The point is that (16) is a quasi-polynomial equation, e.g. a combination of exponentials and polynomials. As proved in [4] for the closely related case of delay differential equations, computing the roots of this class of equations might, in general, give raise to ill-conditioned problems, whose solutions are therefore not very trustworthy. We will refine this comment at the end of the next section.

4.2.
Generalization and the semigroup approach. Given the considerations at the end of the previous section, we would like to recommend the following alternative yet common approach suggested by the results from the theory of semigroups; see e.g. [16] for a comprehensive treatment. This approach is closer to the dynamical systems point of view and, moreover, allows for devising a suitable numerical discretization process in a quite intuitive manner, as presented in Section 5.
To underline the potential of this approach and of the relevant discretization we consider the general linear system of d equations of Gurtin-MacCamy type for an unknown vector population p(a, t) ∈ R d where • N : [0, a † ] → R d×d is a possibly piecewise defined matrix-valued function; • H : X → X is a linear bounded operator of finite rank s acting on the Banach where H(a) = (H 1 (a), . . . , H s (a)) and Kϕ = ( with matrix kernel K 0 ∈ L ∞ [0, a † ], R d×d . Observe that (17) includes those models coming from the linearization around a stationary solution P * (a) of the general nonlinear system of Gurtin-MacCamy type ∂t + ∂P(a, t) ∂a = −M(a, S(t))P(a, t), where the vital rates are given through matrix-valued mortality M : [0, a † ] × R s → R d×d and fertility B : [0, a † ] × R s → R d×d , which depend on the vector of sizes with matrix kernel G ∈ L ∞ [0, a † ], R s×d . In particular p(a, t) = P(a, t) − P * (a),  (14) fitting into the format (21) and (17), respectively. Finally, what follows has been fully described in [7] as far as a single population is concerned, i.e. for the case d = 1. Here we present the theoretical steps for d ≥ 1 which basically remain unaltered.
Problem (17) can be restated as the abstract evolution equation on the state space X above, defined for where A : dom (A) ⊆ X → X is the unbounded linear operator given by Observe that A in (24) is the infinitesimal generator of the strongly continuous is the solution of (17). It can be shown that the solution U (t) of (23) coincides with the solution p (·, t) of (17) whenever p 0 ∈ dom (A). Moreover, the semigroup T is eventually compact and so the spectrum of the infinitesimal generator A consists of eigenvalues only. Therefore, the asymptotic stability of the zero solution of (17) is determined by analyzing the position in C of the eigenvalues of A. In particular, asymptotic stability holds if and only if all these eigenvalues lie in the left half-plane (z) < 0; see e.g. [15]. The main point now is the result that the spectrum of A coincides with the set of characteristic roots, i.e. the solutions of (16); see e.g. [15] for a rigorous theorem or [7] for an explicit proof in the scalar case of (17). This double characterization generalizes the well-known situation for the eigenvalues of a matrix as roots of its characteristic polynomial to the quasi-polynomial case (hence to infinite dimension). As stated in the celebrated work [29] and as mentioned in Section 4.1, computing characteristic roots directly is not numerically advisable in general. This is why we show how to approximate (a finite number of) the eigenvalues of A in Section 5.

4.3.
Fitting of the original model. In this section we show which choices lead the original epidemic model (3) and its linearized version (14) to fit into the format (21) and (17), respectively. First of all, clearly d = 3 and P(a, t) = (S(a, t), I(a, t), R(a, t)) T ∈ R 3 .

The matrix vital rates read
with Q(t) defined in (1). The vector of sizes has dimension s = m + 1, namely with W (t) defined after (10), and the relevant matrix kernel in (22) reads The force of infection λ in M above is unchanged with respect to (9), even though we stress now the dependence on the sizes: As far as the linearization is concerned, the mortality matrix becomes The finite rank operator H introduced in (17) with Q * and B * given by (6) and (7), respectively. 5. Spectrum approximation. We now present how to approximate (a finite number of) the eigenvalues of the infinitesimal generator A defined in (24) as the eigenvalues of a matrix that constitutes a discrete approximation of A. In particular, we adopt the pseudospectral technique guaranteeing high accuracy; it has been proposed to investigate the behavior of solutions for different classes of mathematical models, e.g. delay differential equations in [8], Gurtin-MacCamy equations with nonlinearity restricted to the fertility rate in [6], and general nonlinear Gurtin-MacCamy equations in [7]. Here we extend to systems of equations the scheme developed in [7] for the scalar case.
The method is unaltered in its substance, and, e.g., the proof of convergence follows the same lines as in [7]; the final result on spectral accuracy is stated as Theorem 5.1 below. We revisit the main steps of the algorithm to deal explicitly with the technicalities induced, on the one hand, by a system dimension of d > 1 and, on the other hand, by the fact that we distinguish here the possible agediscontinuities in the model coefficients, i.e. the intervals I l , l = 1, . . . , n, in (8), from those determining (part of) the sizes, i.e. the intervals J r , r = 1, . . . , m.
Assume there exist 0 = a 0 < a 1 < · · · < a n = a † such that the matrix function N in (17), the functions H k and the kernels K k , k = 1, . . . , s, in (18) Observe that θ Such a polynomial is piecewise represented as in [0, a † ] (the same notation is used for the other model functions and kernels). Let us remark that F −1 → I d since a1 0 |l N1,0 (a)|da → 0 as N 1 → ∞. Now we replace the infinite-dimensional operator A on X by the finite-dimensional linear operator A N on X N given by Tedious calculations with careful attention to the indices yield explicit expressions for the entries of the canonical matrix A N ∈ R dN ×dN representing the operator A N . We omit the details here and rather directly give a pseudocode to build A N explicitly in Algorithm 1. There we use for brevity the compact index notation i l := l−1 k=1 N k + i, i = 1, . . . , N l , l = 1, . . . , n, where, according to the sequential ordering of the components in X N , i l refers to the i-th mesh node inside the l-th interval. Recall also the superpositions at the age discontinuities a l = θ (l) As a final observation let us remark that, in general, the integral functionals K 0 in (20) and K k , k = 1, . . . , s, in (19) may not be calculated exactly. If this is the case, an efficient quadrature must be adopted, and we advise a piecewise Clenshaw-Curtis formula with weights ω (r) i and the same Cheyshev nodes θ (r) i used for the main collocation on each interval [a r−1 , a r ], r = 1, . . . , n; see e.g. [28]. Proceeding in this way, in fact, one always obtains for the general integral for k = 0, 1, . . . , m + 1, j = 1, . . . , N r and r = 1, . . . , n. Remark 1. In Algorithm 1 the computation of F −1 is a safe step since the dimension d is sufficiently low in general. If that is not the case, then as a common practice it is better to avoid the computation of the inverse and rather solve the corresponding linear system whenever the product with the inverse is required; for this just a single factorization of F is necessary at the beginning.
We now state for completeness the spectral accuracy result, the proof of which can be obtained by following the convergence analysis in [7]. Let us underline that the convergence of infinite order translates into very high accuracy with rather low matrix dimension, which in turn allows for extremely fast eigenvalues computation by standard algebraic eigensolvers (e.g. eig in MATLAB). This means that a single run for a fixed choice of model parameters and of the relevant equilibria can be efficiently repeated for robust analysis purposes. This is very helpful for producing bifurcation diagrams and detecting stability boundaries in combination with the computation of equilibria. for j = 1 to N r and r = 1 to n do 4: k l Nr,j 5: end for 6: end for 7: for i, j = 1 to N l and l = 1 to n do 8: i )δ ij (δ ij Kronecker's delta) 9: end for 10: for i = 1 to N 1 do 11: for j = 1 to N r and r = 1 to n do 12: for r = 1 to n − 1 do 17:  for r = 1 to n − 1 do 25: end for 27: end for Remark 2. The fact that the error constant C is proportional to |λ| reflects on the use of high degrees for the approximating polynomial when the eigenvalues have, e.g., high frequency. This aspect is confirmed by the numerical experiments reported at the end of Section 6.1, similar to those in [7] for low-frequency eigenvalues.
6. Case studies. In this section we perform the analysis of the stationary solutions of the general model (3) for two specific cases that were introduced and studied in [13] and [17], respectively; we refer the reader interested in the qualitative and mathematical analysis to these two references for biological discussions and conclusions.
6.1. SI epidemics with extra disease mortality [13]. Assume γ = 0 and R 0 ≡ 0 in (3), i.e. susceptible individuals that become infected never recover and either die naturally at age-dependent rate µ(a) or due to the disease at constant rate α. The model then reads ∂I(a, t) ∂t + ∂I(a, t) ∂a = λ(a, t)S(a, t) − [µ(a) + α] I(a, t), with age-specific contamination rate F and age-specific infectiousness g.
In [13] sufficient conditions for the uniqueness of the endemic equilibria have been derived. Due to the limits imposed by the complexity of the model, specific choices of the model parameters were then assumed in order to show numerically the possibility of the existence of multiple endemic equilibria; here we assume different choices, namely a † = 1, α = 10, β ≡ g ≡ q ≡ 1, µ(a) = tan(πa/2), R d 0 = X and Φ(x) = max 1 − x X , 0 for a given X and 3 4 , hence n = 3, m = 1, f 11 = f 31 = 1 and f 21 = 0 in (8). The definition of F gives rise to three age subintervals to be considered for the piecewise numerical approach as discussed in Section 5. The model has the two sizes S * (0) and W * according to (11). The system of two nonlinear equations can be reduced to just a single equation for W * , which is solved by standard Newton iterations; see [13] for details. The result is a function of varying α and, for several choices of X in Φ, it is plotted in Figure  1; this figure shows that, for sufficiently large values of the extra mortality α and of the parameter X, multiple solutions are possible.
Within this framework we consider the above special choices for the parameters and X = 30, letting α vary along the corresponding bold equilibrium curve in Figure 1. For each point on the curve (including those corresponding to the disease-free equilibrium for W * = 0) we are able to compute numerically the right-most characteristic root with the procedure in Section 5. The overall situation is illustrated in Figure 2 (a), where solid and dashed segments denote stable and unstable behavior, respectively. Some corresponding spectra are plotted in the remaining panels (b)-(g) of Figure 2; they confirm the presence of steady-state bifurcations at the points labelled b and e in panel (a) as can be seen in the corresponding panels (b) and (e).
Destabilization due to the emergence of periodic solutions may also occur. An example, for β(a) = q(a) = g(a) = sin (πa) with X = 40 (and all the other previous choices unchanged), is shown in Figure 3, where stability of the equilibrium branch ceases at a Hopf bifurcation, indicated by the dot labelled b in panel (a), and whose spectrum is represented in panel (b).  To give an idea of the computational effort, some data are reported in Table 1. In particular, each spectrum was approximated by using a constant number N of nodes in each of the three age subintervals, leading to a piecewise polynomial of degree N and to a final matrix of dimension dN × dN with d = 2 and N = 3N . The construction of the matrix and the solution of the corresponding eigenvalue problem took about t sp seconds on average on a MacBook Pro 2.53 GHz Intel Core 2 Duo with 4 GB at 1067 MHz. This means that roughly t eq minutes are sufficient to inspect 200 equilibrium points on the selected branch. All spectra of Figures Figure 4 the behavior of the error for some high-frequency eigenvalues of Figure  2 (b), specifically the 4th, 8th, 12th and 16th ones with positive imaginary part, counted starting from the real one. The error is measured by taking as reference the eigenvalues computed with N = 1000, thus following the same procedure used in [7]. In particular, panel (a) illustrates the theoretical spectral convergence stated in Theorem 5.1, showing also the role that |λ| plays in the error constant C. The latter aspect is shown clearly in panel (b), which represents the error for a fixed choice of the piecewise approximating polynomial, in this case N = 500, as a function of |λ|. Let us underline that the result of Theorem 5.1 holds asymptotically with N , therefore it is not a surprise that eigenvalues with large absolute value (or high frequency in general) require rather large values of N . Nevertheless, the choice of very high degree polynomials is not improper from the numerical point of view when dealing with interpolation of analytic functions (as the eigenfunctions of the infinitesimal generator; see also the convergence analysis in [7]) on Chebyshev nodes; see e.g. [28]. 6.2. SIR epidemics with age groups infection transmission [17]. Assume α = 0 in model (3), i.e. infected individuals may die only naturally if not removed by immunity at constant rate γ. Also set fertility to depend only on the age, i.e. R d 0 Φ ≡ 1. Moreover, due to the absence of α, the demographic population N can be assumed to be in its stationary state N * (a) = N * (0)π(a) so that the change of variables x(a, t) = S(a,t) N * (a) , y(a, t) = I(a,t) N * (a) can be performed. Furthermore, since the normalized removed can be obtained as z(a, t) = 1 − x(a, t) − y(a, t), the third equation can be neglected from (3), which can then be verified to read ∂y(a, t) ∂t + ∂y(a, t) ∂a = λ(a, t)x(a, t) − γy(a, t), x(0, t) = 1, y(0, t) = 0, x(a, 0) = x 0 (a), y(a, 0) = y 0 (a), with x 0 (a), y 0 (a) ≥ 0 and x 0 (a)+y 0 (a) ≤ 1; compare e.g. with [22] . As for the force of infection, we assume a standard WAIFW form with g ≡ 1, n = m = 3 with the same age group subdivisions and constant contact rates f lr for all l, r = 1, 2, 3. In if a ∈ [a 1 , a 2 ], where the parameters φ * l , l = 1, 2, 3, satisfy The nonlinear system (26) has been solved for φ * l , l = 1, 2, 3, independently with the package MATCONT and with the function fsolve of MATLAB. The (identical) results are depicted in Figure 5 (a) w.r.t. the parameter φ * 1 of (26) for varying a 2 , which is taken as the bifurcation parameter; other parameter values are f 13 = 2 × 10 4 , f 22 = 2, f 31 = 1, a 1 = 2.5, a † = 50 and γ = 1. The overall resulting structure confirms that, over a central range of a 2 , about (12.991, 13.692), the model actually has three endemic equilibria.
This equilibrium diagram can be enriched with the relevant stability information obtained again through the spectrum approximation of Section 5. The overall bifurcation picture is represented in Figure 5 (a) for φ * 1 as a function of a 2 . The spectra corresponding to (some of the) bifurcation points are represented in the remaining panels (b)-(e) of Figure 5. For instance, looking at panel (a) from the left for increasing a 2 , one sees that at a 2 9.012 a Hopf bifurcation occurs with the emergence of a periodic solution that, in the supercrititcal case, would inherit the asymptotically stability; see panel (e). A deeper analysis would be required now for the study of the asymptotic stability of the periodic solutions, since multiple complex-conjugate pairs of roots appear in the right half-plane as can be seen in panel (d). Again, to give an idea of the computational effort, some data are reported in Table 2, whose entries have the same meaning of those in Table 1. In particular, d = 2, N = 3N and all spectra of Figure 5 were computed with N = 200.  Table 2. Computational data for the SIR epidemic model. 7. Conclusions and further developments. By considering a prototype model of SIR type epidemics with age-structure, we have shown how effective numerical methods for the approximation of the characteristic roots of equilibria in an infinite dimensional setting can be devised. The resulting numerical tools can be rather helpful from different points of view. First, they can be used for obtaining numerical confirmation in those cases where analytical statements can be proved; these are very limited situations in practice, often restricted to the local analysis of the trivial equilibrium. Second, these tools can be used conversely for inspecting the system behavior for wide ranges of the parameters, in order to get at least an approximative idea of what is going on; in this way the theoretical analysis can be directed towards specific aspects of the dynamics. An illustrative example of this is described in [3], where the rigorous proof of the local asymptotic stability for the positive equilibrium for any choice of the system parameters of an SEIR model was attempted and achieved only after an extensive numerical check for unstable regions in the parameters plane. Third, the numerical tools may be useful also in disproving conjectures and closing open problems by providing valid (counter-) examples; an example is the analysis in [17], where the presence of multiple endemic equilibria for a SIR transmission disease was shown, contrary to what was suggested in [19]. Furthermore, the possibility of destabilization of equilibria to limit cycles through Hopf bifurcation was established, hence answering an open question raised in [22].
Of course it must be said that, even though the methods adopted are always accompanied by a rigorous analysis of their convergence properties, still the infinitedimensional nature of the problem under consideration is responsible for the lack of exact analytical results. In our opinion this fact alone is enough to require the cross-checking of the numerical outcomes by independent techniques. As far as this is concerned, the presented approach based on the discretization of the associated infinitesimal generator (a direct method) should be accompanied by simulations via standard time-integrators (an indirect method); clearly the efficiency of direct methods (due to their spectral accuracy; see Theorem 5.1) is unbeatable when compared to very long simulations, especially when parameters are varied as well.
A comment on the work presented here may be that the analysis of equilibria is only the first step when one wants to give a general bifurcation picture of the overall system behavior. We certainly agree with the statement. For this reason, and motivated by many complicated situations we have experienced in several classes of dynamical systems arising from PDEs and FDEs, we are directing our current and future efforts towards the analysis of the stability of periodic solutions (through the computation of Floquet multipliers; see e.g. [11]) as well as to the detection of chaotic dynamics (through the computation of Lyapunov exponents; see e.g. [5]). Both techniques concern non-autnomous linearized problems and, therefore, the attention should be moved from the discretization of the infinitesimal generator to that of the evolution family (a two-parameters dependent semigroup; see e.g. [14]).
As a final remark, we would like to stress that the development of these methods and the associated algorithms should be followed by the necessary effort to convert them to suitable user-friendly software packages. This will allow the practitioners to use these tools. This is a relevant point, because, after all, we started in the introduction from the observation of natural phenomena, which is rarely driven by numerical analysts. With this in mind, we invite the interested reader to contact the authors for the Matlab codes used in this paper.