key: cord-0957592-700dxcb7
authors: Breda, Dimitri; Florian, Francesco; Ripoll, Jordi; Vermiglio, Rossana
title: Efficient numerical computation of the basic reproduction number for structured populations
date: 2020-08-27
journal: J Comput Appl Math
DOI: 10.1016/j.cam.2020.113165
sha: 7fcc0e5e2ba4c494759931bb1b298842c0a3680a
doc_id: 957592
cord_uid: 700dxcb7

As widely known, the basic reproduction number plays a key role in weighing birth/infection and death/recovery processes in several models of population dynamics. In this general setting, its characterization as the spectral radius of next generation operators is rather elegant, but simultaneously poses serious obstacles to its practical determination. In this work we address the problem numerically by reducing the relevant operators to matrices through a pseudospectral collocation, eventually computing the sought quantity by solving finite-dimensional eigenvalue problems. The approach is illustrated for two classes of models, respectively from ecology and epidemiology. Several numerical tests demonstrate experimentally important features of the method, like fast convergence and influence of the smoothness of the models’ coefficients. Examples of robust analysis of instances of specific models are also presented to show potentialities and ease of application.

Aim of this paper is to give a concise, efficient and effective method for the numerical computation of the basic reproduction number, a fundamental tool in the analysis of ecological and epidemiological models of population dynamics, as well as in tackling evolutionary aspects in those populations, see e.g. [4, Section 2.1] and the references therein. We refer to [17] as a starting reference on this quantity, commonly denoted R 0 .

The need for a numerical approach is soon justified since we are interested in population models posed on infinite-dimensional spaces, where R 0 is usually characterized as the spectral radius of suitable operators, mainly known as next generation operators, see, e.g., [11] . As in general there is no explicit expression for R 0 readily practicable for computation (and in some cases nor the relevant operators are explicitly known), the main idea is to reduce to finite dimension through some discretization techniques. Then approximations to R 0 are obtained by computing the (dominant) eigenvalues of the resulting matrices.

As far as this general strategy is concerned, to the best of our knowledge [25] is the first scheme appearing in the literature. Therein an elementary discretization based on the celebrated Euler's method is proposed to approximate the basic reproduction number of a class of age-structured epidemics. Very recently a slight improvement has been presented in [16] , obtained by employing θ-methods, always in the context of age-structured epidemics. Made exception for [14] , which is anyway related to the approach proposed here (see below), [16, 25] are the only two works concerning the proper numerical computation of R 0 in continuously structured populations.

Here, on the one hand, we intend to largely improve the seminal idea of [25] by applying state-of-the-art discretizations based on pseudospectral polynomial collocation (see, e.g., [1, 2] ). The general outcome is a more reliable tool, with faster convergence ideally of infinite order (known as spectral accuracy, see, e.g., [29] ). This represents a great advantage compared to the finite-order convergence of [16, 25] , as it translates into much more accurate approximations obtained with much smaller matrices. As a consequence the required computational load reduces, in terms of both time and memory usage. The latter is a favorable feature when stability and bifurcation analyses are the final target in presence of varying or uncertain model parameters, and of course this is frequently the case in realistic contexts.

On the other hand, we aim at increasing the flexibility of the numerical approach, widening the applicability also to other models, e.g., from ecology. In particular, we illustrate two paradigmatic classes of models, namely spatially-structured cell populations (whose prototype representative is referred to as model A in the sequel) and infectious diseases in age-structured populations (whose prototype representative is referred to as model B in the sequel). This choice potentially serves as a basis for extending the proposed method to more involved systems of population dynamics (e.g., with several structuring variables). If a structured epidemic model describing the novel pandemia of coronavirus was available, such as Model B or any of its extensions, this method could be applied to compute R 0 for different values of the parameters.

The contents of this research are organized as follows. In Section 2 we present the two classes of models mentioned above, introducing all the necessary ingredients and the operators leading to the definition of R 0 . The numerical treatment is proposed in Section 3, where explicit expressions for the discretizing matrices are constructively given, also for the sake of implementation for those interested (Matlab codes are freely available at http://cdlab.uniud.it/software). In Section 4 we first illustrate the numerical properties of the proposed techniques, analyzing experimentally error and convergence, and then we finally perform some robust analysis on specific instances of the models of interest. In Section 5 we draw conclusions and we present ongoing investigations and future developments of the approach.

Let us remark that here we deliberately quit to tackle a theoretical investigation of the convergence features of the proposed methods. This would indeed require suitable tools, falling out of the scope of the current experimental study, as well as enough room for a proper description. The formal proofs as well as possible improvements are given in the related work [7] . As far as other developments are concerned, we also plan to extend this approach in order to treat periodic problems, that is, to population models where the environment is time-periodically varying. In this respect we suggest the reading of the very recent work [22] .

Finally, we wish to mention that the present research originated from [13] , where an embryonic study of collocation techniques for the numerical computation of R 0 is proposed (also the work [14] mentioned above is inspired from [13] ). Moreover, let us also note that pseudospectral techniques have become a reference tool in the numerical treatment of infinite-dimensional dynamical systems in view of either stability and bifurcation analyses in the context of population dynamics, see, e.g., [6] or [8] and the references therein.

Continuously structured models of population dynamics can be described by nonlinear abstract differential equations for the density of individuals with respect to some structuring variables, e.g., age, size or space.

With the aim at setting the theoretical background, let X be a Banach lattice and let u ∈ X represent the density of individuals. We introduce next a continuous-time evolution equation for the population density. For the theory of population dynamics we refer, among many others, to the monographs [10, 18, 19, 20, 21, 27] .

On the one hand, if we focus on ecological models (the class of models A in the sequel), then in general a nonlinear problem can be decomposed into birth and nonbirth terms as

where, for any u ∈ X, B(u) is the birth linear operator and M(u) is the linear operator accounting for other processes non related to birth (e.g., mortality or transition). Of course, this decomposition depends on what is exactly considered as a birth event, see an example in Section 2.1 and also [3, 9] . If we are concerned with the extinction of the population in (2.1), then we can focus on the steady-state u * ≡ 0. Its stability is analyzed by means of the (formally) linearized model

On the other hand, if we focus on epidemiological models (the class of models B in the sequel), then we firstly need to distinguish between infective individuals u ∈ X and (multi-type set of) non-infective individuals v ∈ X m , with m some positive integer. Here, the original nonlinear problem consists of m + 1 equations for the population densities of the two types of individuals. The one relevant to the infective individuals can be decomposed according to infection and non-infection terms as

The remaining ones have the general form

with F depending on the type of epidemics. For any v ∈ X m and u ∈ X, B(v, u) is the linear infection operator and M(v, u) is the linear operator describing processes other than infection (e.g., permanent or temporal recovery or transition processes). Again, the decomposition may differ depending on what is considered as an infection event.

Here, typically, the main concern is the early stage of the infection, i.e., the initial epidemic spread in (2.3)-(2.4). We can thus focus on the disease-free steady state (v * , u * ) ∈ X m + × {0} (where X + stands for the positive cone in X). Since the (formally) linearized system has triangular form, under common assumptions, the stability analysis of the disease-free equilibrium is reduced to the linear model

(2.5) Summarizing and simplifying the notation to make it uniform in (2.2) and (2.5), for either model A or B we are faced with the analysis of an abstract linear equation of the form

We recall that above X is a Banach lattice where the population density lives, B : X → X is a linear operator accounting for the birth process (meant as proper birth, or infection) and M : D(M) ⊆ X → X is a linear operator accounting for the remaining processes, which we call mortality for brevity (thus meant as proper mortality, or recovery or, in general, any stage transition as in [27, Part 2] ). Typically, the domain D(M) of M is a subset of a subspace Y ⊆ X where some degree of smoothness is present, characterized by additional constraints in general described through a linear map C : X → R p for some positive integer p:

Other common (biologically meaningful) assumptions are required. In particular, as in [4] , B is positive and bounded, while −M is the infinitesimal generator of a strongly-continuous semigroup {T(t)} t≥0 of positive linear operators. Its spectral bound s(−M) is strictly negative, which accounts for extinction in absence of birth. Consequently, 0 belongs to the resolvent of M and thus the latter is invertible with

For these and other aspects on (positive) linear operators see, e.g., [26, 30] .

It is worth to mention that although (2.6) represents a rather large family of population models (including in particular those with discrete structure when dim X < +∞), see, e.g., [10, Section 7.2]), not all the structured models can be cast into this form, see the discussion in [4] . However, in this standard case we can define the next generation operator as the operator BM −1 : X → X acting on the same (whole) space of the population density and, eventually, characterize the basic reproduction number as its spectral radius:

Let us remark that the next generation operator is a generalization to the infinitedimensional setting of the well-known next generation matrix for finite-dimensional models (i.e., dim X < +∞). Moreover, although the theoretical framework for the basic reproduction number has been well-established by many authors, see [10, 11, 17, 21] and the references therein to name a few, there is a lot of room for its efficient computation as already anticipated in the introduction. It follows from definition (2.8) that R 0 is a non-negative spectral value being BM −1 positive and bounded, see, e.g., [26] . It is actually the largest λ ≥ 0 for which Bφ − λMφ = ξ cannot be solved uniquely for φ ∈ D(M) once ξ ∈ X is given, see, e.g., [4] . Moreover, if it happens that this λ is a generalized eigenvalue, then

for some nontrivial generalized eigenfunction φ ∈ D(M). In particular, if BM −1 is also compact, then the Krein-Rutman theorem [24] ensures that R 0 is a positive eigenvalue (for slightly weaker, yet involved assumptions see [21, Section 10.2] ). At this point let us anticipate that the numerical method proposed in Section 3 relies on this assumption, so that (2.9) can be reduced to a standard (read finite-dimensional) generalized eigenvalue problem for matrices. Compactness is proved in [7] under reasonable assumptions on the specific ingredients of models A and B as described in the forthcoming Sections 2.1 and 2.2. A discussion on further aspects related to R 0 being or not a generalized eigenvalue and their relevance on the numerical outcome is left to Section 4.

Following definition (2.8), with any further assumption, we obtain the useful upper bound

As a first representative of prototype models to illustrate our numerical approach, let us consider a population of bacteria proliferating and moving along the intestine of an animal host, see, e.g., [3, 4] . Given that the intestine is portrayed as the line segment [0, l] of finite length l > 0, let us set X := L 1 (0, l) and let u(·, t) ∈ X be the spatial density of bacteria along the intestine at time t ≥ 0. The nonlinear problem is described as

is the population size weighted through a nonnegative weight σ. Above, c is the velocity of the flow, D is the diffusion coefficient, β is the fertility rate and µ is the mortality rate. All are non-negative functions of the position x ∈ [0, l] and of the size S. Note that transport, diffusion and vital processes are density-dependent accounting for limited resources (crowding effects) and spacespecific accounting for the heterogeneity of the intestine, see [3, 4, 20] .

For a population of bacteria, one can consider either symmetric division, i.e., when a mother cell divides then two daughter cells are born and the former disappears, or asymmetric division, i.e., one mother and one daughter. Both types can be considered simultaneously in (2.12) by setting 1] representing the probability of asymmetric division. For instance, for space-independent vital rates one would readily get that

since integrating over the whole space in (2.12), the system reduces to an ODE for the population size, see [4] . The latter is a clear example of the fact that the basic reproduction number depends on what is exactly considered as a birth event. In any case, here we assume symmetric division (ν = 0) for simplicity, and without loss of generality with respect to the numerical discretization.

In what follows we are concerned with the extinction of the bacterial population. Linearizing (2.12) around the extinction equilibrium u * ≡ 0 gives

0) to lighten the notation though with a little abuse. Moreover, we further assume c(x) > 0 for all x ∈ (0, l), as well as β and µ bounded with their sum not identically vanishing. Concerning diffusion, we consider either D(x) ≥D, x ∈ [0, l], for some positiveD (everywhere positive diffusion) or D ≡ 0 (complete lack of diffusion).

Eventually, with respect to the reference linear equation (2.6), it is natural to define from (2.13) the birth operator B : X → X as 14) and the mortality operator M :

Then, following (2.10), the upper bound

holds, as expected due to symmetric division.

As a second representative of prototype models to illustrate our numerical approach, let us consider the spread of an epidemics in an age-structured population. We split up the individuals of the population into three classes according to the stage of the infectious disease. If l > 0 denotes the maximum (chronological) age, we set X := L 1 (0, l) and let s(·, t), i(·, t) and r(·, t) in X denote the density with respect to age of, J o u r n a l P r e -p r o o f Journal Pre-proof respectively, susceptible, infected and removed individuals at time t ≥ 0, see, e.g., [18, 21] . The nonlinear problem is described as

Above, χ is the probability of infection transmission through an infectious contact and c is the density-dependent contact rate of a susceptible individual of age x with an infected individual of age y. Moreover, β 0 , µ, γ, δ and σ are, respectively, the fertility, (natural) mortality, removal, recovery and loss of immunity rates, all non-negative and with β 0 not identically vanishing. Finally, θ andθ are the probabilities of vertical transmission of infectiveness and immunity, respectively. Note that recovery, removal, immunity loss and vital processes are all age-specific during the lifespan [0, l].

System (2.18) includes different types of epidemic models. For instance, the SIS model is given by γ ≡ σ ≡ 0; the SIR model is given by δ ≡ σ ≡ 0; the SIRS model is given by δ ≡ 0. The exact interpretation of the removed class depends on the type of model considered. Moreover, we could also consider a vaccination rate and the resulting system would have an extra transition, from susceptible to recovered individuals in this case, or the consideration of exposed individuals, or even we could consider a multistrain epidemic model 1 . See [23] and [18, Chapter VI] , [21, , [27, Chapter 22] .

From now on, we assume that the population has reached a demographic steady state. Indeed, we assume zero demographic growth, i.e.,

, together with Π 0 (l) = 0, defines the survival probability, as well as initial conditions such that s 0 (x) 1 To consider the exposed individuals or the multistrain epidemics, we need u ∈ X n , with n some positive integer, see [7] . 8 

In what follows we are concerned with the early stage of the epidemics of an infectious disease of any type (SIS, SIR or SIRS). Linearizing (2.18) 

for the infected individuals. We let l 0 µ(x) dx = +∞ to avoid immortal individuals. In order to properly deal with this condition, we make the change of variables i(x, t) = u(x, t)Π 0 (x), thus removing the mortality term and reducing the problem to

and effective fertility rate β(

We recall that l 0 β(x) dx = 1 due to (2.19) . Note that here u(l, t) can be strictly positive allowing for the possibility of chronic infected individuals. 

with domain 

with

Above, we implicitly assume that θΠ 1 ≡ 1 in order for the last factor to be welldefined. Finally, following (2.10), we obtain the upper bound

representing the total number of infections taking into account vertical transmission. Eventually, the rhs of (2.26) being less than 1, can be used as a sufficient condition for the disease eradication.

Let N be a positive integer. Under the assumption that R 0 is a generalized eigenvalue, we propose a general discretization approach consisting in the collocation of (2.9) on a mesh of N + 1 distinct nodes 0 =: x N,0 < x N,1 < · · · < x N,N := l distributed on [0, l], l > 0. Let us soon observe that collocation is meaningless in the whole space X = L 1 (0, l), but the method is restricted to generalized eigenfunctions φ, which are, in general, smooth enough to guarantee pointwise definition, recall indeed (2.7).

In the sequel, let X N := R N+1 be the finite-dimensional counterpart of X and let Φ :

The choice of including x N,0 = 0 and x N,N = l in the collocation mesh is justified by the fact that, in general, φ belongs to some restricted domain (2.7) where the mortality operator M is suitably defined. For the models of interest in this work, in fact, this domain is characterized by boundary conditions at one or both the extrema 0, l of the domain of the function space X, see indeed (2.16) and (2.24). Therefore, the choice above allows for a simplified treatment as one can directly relate the boundary conditions to the collocation equations.

As the concerned mortality operators usually involve differentiation and/or integration, see and, for the same v above, w T N v is an approximation to l 0 f (x) dx. Both follow straightforwardly from approximating f with the N-degree interpolating polynomial

Even though at this moment it is not necessary to specify any particular choice of nodes, let us remark that in the case of Chebyshev extremal points, as we assume in Section 4, H N and w N can be obtained rather easily 5 . For these and other related aspects see, e.g., [29] . See also [12] for a recent review of results on H N .

The collocation method leads to a discrete version

of the generalized eigenvalue problem (2.9), where the structure of the matrix representation of the finite-dimensional operators B N and M N depends on the specific model as detailed in the following Section 3.1 for model A and Section 3.2 for model B. In any case, the discrete problem has finite dimension, as it is posed on X N . Of course, as a main outcome one expects that the eigenvalues of (3.2) approximate (part of) the eigenvalues of (2.9), the accuracy improving as N increases. This is shown to be the case in Section 4 by way of several numerical experiments, properly tuned to test the error behavior, to measure the convergence rate and to detect other peculiarities of the proposed approach. Note that the eigenvalues of (3.2) can be computed with standard techniques for finite-dimensional generalized eigenvalue problems (e.g., Matlab's eig, based on the well-known QZ algorithm, see, e.g., [15] ). Let us remark that we are anyway interested in the spectral radius of the next generation operator, so that we are mostly concerned with the dominant part of its spectrum.

Let us finally stress again that the use of the proposed methodology to approximate R 0 is based on the assumption that this number actually corresponds to a generalized eigenvalue. If this is not the case, already (2.9) looses sense. Nevertheless, in Section 4, we report on some tests under the latter hypothesis, where the scheme is still able to give reasonable approximations (even though with slower convergence than what experimented in the case of generalized eigenvalues).

Let us recall the main ingredients (2.14), (2.15) and (2.16) of the class of cell population models, model A in Section 2.1. The generalized eigenvalue problem (2.9) reads

The proposed collocation consists in looking for an N-degree polynomial p N satisfying

By using the Lagrange representation 

Note how the first and last rows of B N and M N realize the discrete boundary conditions (3.5), simulating those characterizing D(M) in (2.16), i.e., (3.4). 

The proposed collocation consists in looking for an N-degree polynomial p N satisfying 

A first series of numerical tests is presented in Section 4.1 with the aim of illustrating the convergence properties of the proposed techniques, as well as related aspects such as, e.g., the effect of the smoothness of the models' coefficients. Then in Section 4.2 we use the tested algorithms to perform some quantitative studies in the context of varying parameters, checking whether the outcomes confirm the theoretical expectations.

All the experiments are run on a MacBook Pro 2.3 GHz Intel Core i7 with 16 GB memory, through codes written in Matlab version R2019a (freely available at http://cdlab.uniud.it/software as anticipated in the Introduction). In these codes we always use Chebyshev extremal points as collocation nodes, i.e., x N,i = l[1 − cos(iπ/N)]/2, i = 0, 1, . . . , N. We refer to [29] for their numerous advantageous properties in the context of numerical interpolation, differentiation and quadrature.

The numerical properties of the proposed approach are mainly investigated by analyzing the behavior of the error |R 0,N − R 0 | between the approximated and the exact basic reproduction numbers. Here for exact we mean either the theoretical value of R 0 when this is explicitly available, or a reference value R 0,N computed with a suitably largeN when an exact value is not attainable.

A second error is also presented when possible, viz. p N − φ ∞,M . This is the error between the exact generalized eigenfunction (if available) and its collocation approximation, measured as the maximum absolute value of their differences on a mesh of M equidistant points in [0, l]. In particular, we always use M = 1 000. The collocation polynomial p N is reconstructed from the computed generalized eigenvector Φ associated to the dominant generalized eigenvalue of (3.2) through barycentric interpolation by applying the algorithm proposed in [5] . Of course, the same normalization is prescribed from time to time for both p N and φ.

Finally, the discrete generalized eigenvalue problem (3.2) is solved either by Matlab's eig or eigs. We also test the discretization For all the instances of model A we set l = 1, c(x) = l − x, x ∈ [0, l], and

The other ingredients differ as described next.

(A1) [the immortal case] LetD = 1, µ ≡ 0 and

With these choices the next generation operator is compact according to [7] and thus R 0 is a generalized eigenvalue. It is not difficult to recover R 0 = 2 from (2.11). Moreover, some calculations that we omit allow to show, starting from (2.9), that the corresponding generalized eigenfunction is

normalized as φ(0) = 1, and the eigenfunction in (4.2) is ψ = βφ. Let us remark that this case is not realistic, since we are neglecting the cell mortality, but it is useful as a numerical test. Independently of these choices, R 0 is unknown, as well as the corresponding generalized eigenfunction.

J o u r n a l P r e -p r o o f so that also Π 1 (x) = l−x l , x ∈ [0, l], follows. We assume absence of vertical transmission, i.e., θ = 0, which makes useless to specify β. With these choices the next generation operator is again compact according to [7] . Since K above is the product of functions in each of the variables x and y, the next generation operator in (2.25) becomes a rank one operator and we can explicitly compute

.

As far as α is concerned, we consider (B2.1) α = 1/4: life expectancy is 0.8l and R 0 = 2;

(B2.2) α = 1: life expectancy is 0.5l and R 0 = 1.59375.

Finally, and independently of α, the eigenfunction corresponding to R 0 is ψ(x) = x 2 (l − x) 2 , normalized as ψ(l/2) = l 4 /16, and the generalized eigenfunction is

which, for the choices above, turns out to be a polynomial of degree 5, viz. 

As a first experiment we test the convergence properties of the proposed approach for case (A1). The spectrally accurate behavior (namely the error decays faster than O(N −k ) for any natural k, [29] ) is evident in Figure 4.1 (left) for the approximation of both R 0 and the relevant generalized eigenfunction φ. To note that for the former the approximation is already very good for low even values of N, due to the (anti-) symmetry properties of the Chebyshev differentiation matrix H N [29] . Similar trends emerge also in the right panel, where (4. above, but for case (A2): the results are qualitatively the same.

In the previous tests both R 0 and the relevant generalized eigenfunction are analytically known. As the same is not possible for case (A3), we show in Figure 4 .3 only the error on R 0 with respect to a reference value R 0,N as explained at the beginning of Section 4.1. As expected, the convergence is spectrally accurate for case (A3.1), in whichD = 1 ensures compactness of the next generation operator (according to [7] ). Theoretically, alsoD = 10 −6 = 0 in case (A3.2) guarantees compactness, but it is clearly visible from the plot that much larger values of N are necessary to start appreciating the spectral accuracy. It is indeed reasonable to expect that the value of D affects the error constants, causing their increase asD → 0, still being the problem compact. When we deal instead with case (A3.3), in which the absence of diffusion causes the loss of compactness, convergence still occurs, even though at a fixed rate (seemingly linear). The fact is somehow surprising (and certainly merits future investigation), given that without compactness (2.9) may even become meaningless and we are thus using a finite-dimensional eigenvalue problem to approximate components of the spectrum possibly other than the point one. As a final remark for model A, let us mention that all the eigenfunctions computed with the proposed approach in case (A3) turned out to approximate seemingly smooth curves (as they are expected to be). To avoid redundancy, we give examples only for case (B3) in Section 4.2.

As far as case (B1) is concerned, the results of pending on α. First of all, let us recall that the generalized eigenfunction (4.7) is a polynomial of degree 5, justifying the sudden drop of the error to machine precision occurring with N = 6 in both the left and right panels. The theoretical proof of convergence fully elaborated in [7] shows that this indeed happens at N = 6 and not already at N = 5 as one may expect. As far as the approximation of R 0 is concerned, instead, the same behavior is ensured if enough regularity is provided, as it is the case for α = 1 in the right panel. In fact, for instance, it is evident from (4.5) that Π 0 is smooth for α ≥ 1, but for α < 1 it becomes rational and blows up at x = l. Indeed, the choice α = 1/4 in the left panel prevents the method to perform the standard spectral accuracy, and convergence of fixed order (seemingly 4) occurs. That the regularity of the model ingredients plays a role in the convergence analysis is a general fact, but for the specific problem currently investigated we leave the elaboration of the necessary details to [7] . 

To conclude, we give some examples of practical application of the proposed techniques to analyze the behavior of the chosen models in the presence of varying parameters. Below, the choices of the discretization index N are instructed by the results discussed in the previous section. Moreover, the intervals of variation of the concerned parameters are discretized by using uniform meshes of P points. Corresponding values of R 0 are repeatedly approximated with the chosen N for each of these points, or of their Cartesian product in the presence of two varying parameters. In the latter case we adopt standard contouring algorithms, viz. Matlab's contourf. Besides the graphical output, we give also some indication of the overall computational time with respect to the parameter(s) grid size P. Let us anticipate that all the tests concern the general instances of either model A or B, namely cases (A3) or (B3) as described in Section 4.1.

In Figure 4 .7 we investigate for case (A3) how R 0 varies as a function of either global fertility parameterβ (left panel) or global diffusion parameterD (right panel). For the formerμ = 1 and both cases (A3.1) and (A3.2) are considered. For the latter µ = 50 andβ = 10. As expected, R 0 grows with respect to the fertility from zero up to the asymptote R 0 = 2, since we assumed symmetric division (2.17), with larger diffusion values favoring the speed of increase. Instead interestingly, at fixed fertility, there is a single peak with respect to the diffusion, i.e. there exist an intermediate diffusion maximizing the basic reproduction number. At this point, a naive approach is that a population adopting such strategy would persist and become an unbeatable population. 11 Computationally speaking, both figures are obtained by repeatedly computing 11 The role of the diffusion from the evolutionary point of view is out of the scope of the present work.

J o u r n a l P r e -p r o o f R 0 with respect to the varying parameter upon discretizing the latter on a mesh of P = 1 000 points, procedure that took on average just few seconds (say 1 to 5). We can also combine the above results in a single plot, Figure 4 .8, showing the level curves of the surface R 0 = R 0 (D,β). The computation was performed on a grid of size 105 × 100, and took just few minutes (say less than 10). As for case (B3), in Figure 4 .9 we show the variation of R 0 as a function of the vertical transmission probability θ, confirming the biological expectation that this function is monotonically increasing as it is explicitly for case (B1). Notice that an increase of the vertical transmission can lead to the spread of the infection. The grid size is P = 1000 again, and the computation took less than 2 minutes.

Finally, always for case (B3), we show in 

We have considered structured population models whose linearization at the equilibrium admits the abstract formulation (2.6), where the birth and the mortality operators are decomposed. In this context we have proposed an efficient and flexible numerical method to compute R 0 , which is based on the pseudospectral discretization of both operators and on the solution of the resulting generalized matrix eigenvalue problem.

We have experimentally analyzed the convergence for two paradigmatic classes of models, arising in ecology and epidemiology. The numerical results give evidence that the convergence rate with respect to the discretization parameter N depends on the regularity of the model functions and so, when they are analytic, the spectral accuracy occurs. This typical convergence behaviour of the pseudospectral approximation of R 0 has been confirmed by the complete theoretical analysis carried out in [7] and it gives a great advantage in comparison to the fixed order methods proposed in [16, 25] , since we can obtain accurate results where N is about few tens. The limited computation cost allows also to analyze the behavior of the models in the presence of varying parameters.

Looking forward, further extensions of the numerical method will concern population models where the environment is time-periodically varying [22] .

Finally, let us also recall that R 0 is an alternative to the Malthusian parameter, i.e., the exponential population growth rate. Concerning computations, in general it can be more difficult to deal with the latter, i.e., the spectral bound r := s(B − M) of the complete generator, rather than with the basic reproduction number R 0 . For instance, we cannot ensure that r is a spectral value (not even an eigenvalue) and there can be no upper bound for r as it is always the case for R 0 (corresponding to a sufficient condition for population extinction or infection eradication). Moreover, the rank of B is typically lower than that of B − M, which is an advantage for R 0 . Finally, extra evolutionary aspects added to the models are more often related to R 0 than to r. Summarizing, besides the well-known sign relation between these quantities, see, e.g., [28] , the basic reproduction number offers diverse advantages over the Malthusian parameter, from both the theoretical and the numerical points of view.

Collocation techniques for structured populations modeled by delay equations

15 years or so of pseudospectral collocation methods for stability and bifurcation of delay equations

On the reproduction number of a gut microbiota model

A practical approach to R 0 in continuous-time ecological models

Barycentric Lagrange interpolation

Pseudospectral discretization of nonlinear delay equations: new prospects for numerical bifurcation analysis

Collocation of next-generation operators for computing the basic reproduction number of structured populations

Stability of linear delay differential equations -A numerical approach with MATLAB. SpringerBriefs in Control, Automation and Robotics

The many guises of R 0 (a didactic note)

Mathematical Tools for Understanding Infectious Disease Dynamics. Theoretical and Computational Biology

On the definition and the computation of the basic reproduction number R 0 in models for infectious diseases in heterogeneous populations

Pseudospectral discretization of delay differential equations in sun-star formulation: results and conjectures

Numerical computation of the basic reproduction number in population dynamics

PC-based sensitivity analysis of the basic reproduction number of population and epidemic models

Matrix computations. Johns Hopkins Studies in Mathematical Sciences

A theta-scheme approximation of basic reproduction number for an age-structured epidemic system in a finite horizon

A brief history of R 0 and a recipe for its calculation

Mathematical Theory of Age-Structured Population Dynamics

The Basic Approach to Age-Structured Population Dynamics -Models

An Introduction to Mathematical Population Dynamics -Along the trail of Volterra and Lotka. Number 79 in La matematica per il 3+2

Age-Structured Population Dynamics in Demography and Epidemiology

The basic reproduction number R 0 in time-heterogeneous environments

Contributions to the mathematical theory of epidemics III -Further studies of the problem of endemicity

Linear operators leaving invariant a cone in a banach space

Numerical approximation of the basic reproduction number for a class of age-structured epidemic models

Grundlehren der mathematischen Wissenschaften

Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity

Spectral methods in MATLAB. Software -Environment -Tools series

The Asymptotic Behaviour of Semigroups of Linear Operators

For all the instances of model B we set again l = 1 and the rest is defined below.(B1) [the age-independent case] 8 Let µ ≡μ = 28, β 0 ≡ µ (so that (2.19) holds), γ ≡γ = 1 and δ ≡δ = 1 and K(x, y) :=kμe −μy , x, y ∈ [0, l], fork = 52. It also follows thatFinally, let θ = 1/7. With these choices the next generation operator is compact according to [7] . Due to the finite maximum age there is a very small fraction of immortal individuals, but still we can explicitly obtain R 0 k /[(1 − θ)μ + γ +δ] = 2. It is evident from (2.25) that the corresponding eigenfunctions ψ in (4.2) are the constant functions given that K is independent of x. Then one can recover through (2.23) the generalized eigenfunctions Note that l 0 Π 0 (z) dz = l α+1 . Let moreover γ and δ satisfy 9 (4.6)