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.
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. We thus reserve to give formal proofs as well as possible improvements in a related work, currently in preparation [8] . 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 timeperiodically 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 [7] 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 [8] 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
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.
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, 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. 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.,
where Π 0 (x) := exp −
x 0 µ(y) dy , x ∈ [0, l], together with Π 0 (l) = 0, defines the survival probability, as well as initial conditions such that s 0 (x) + i 0 (x) + r 0 (x) = N 0 Π 0 (x)/ l 0 Π 0 (y) dy, x ∈ [0, l], for some positive N 0 . Therefore, the total population
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 
with domain
Moreover, we can explicitly compute the (normalized) next generation operator BM −1 (and even obtain explicit expressions for some special cases, see Section 4). Indeed, by defining Π 1 (x) := exp − x 0 [γ(y) + δ(y)] dy ≤ 1, the next generation operator becomes
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 Φ := (Φ 0 , Φ 1 , . . . , Φ N ) T ∈ X N with Φ i representing the numerical approximation of φ(x N,i ), i = 0, 1, . . . , N.
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
which indeed satisfies p(x N,i ) = f (x N,i ), i = 0, 1, . . . , N. 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 4 . 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 
As far as the models to test are concerned, below we list all the specific choices by giving the defining rates and coefficients. In parallel we also give some of their key features.
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 [8] 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 Independently of these choices, R 0 is unknown, as well as the corresponding generalized eigenfunction.
For all the instances of model B we set again l = 1 and the rest is defined below.
(B1) [the age-independent case] Let µ ≡μ = 28 (so that e −μl 0), β 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 that
Finally, let θ = 1/7. With these choices the next generation operator is compact according to [8] . 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
with φ(0) ψ/182 following from (3.7), with ψ playing the role of normalizing factor.
(B2) [the proportionate case] Let
fork = 16 065/64 and
for some α > 0. Note that l 0 Π 0 (z) dz = l α+1 . Let moreover γ and δ satisfy
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 [8] . 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.
(normalized as φ(l/2) = 5l 5 /384).
Let moreover γ and δ be as in case (B1). Here, infection events occur mostly between individuals of the same age and there are chronic infected individuals since e −(γ+δ)l > 0. Finally, let θ = 0.5 and β(x) := b(x)Π 0 (x)/ l 0 b(y)Π 0 (y) dy, x ∈ [0, l], for b(x) := (x/l) 2 e −6x/l , which ensures (2.19) and maximum fertility at x = l/3. For these choices the next generation operator is compact 8 , but both R 0 and the corresponding generalized eigenfunction are unknown.
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. 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 [8] ). 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 Figure 4 .4 show spectral accuracy in the approximation of both R 0 and the relevant generalized eigenfunction as expected.
Slightly dissimilar is the situation for case (B2), where we still get convergence, but with trends different from what experimented so far and, moreover, possibly depending 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 [8] 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 [8] . Finally, the analysis of case (B3) reveals spectral accuracy again, Figure 4 .6. Recall that in this case, as for case (A3) of the preceding section, neither R 0 nor the relevant generalized eigenfunction are known.
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. 9 Computationally speaking, both figures are obtained by repeatedly computing 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 