Efficient numerical computation of the basic reproduction number for structured populations
Dimitri Breda, Francesco Florian, Jordi Ripoll, Rossana Vermiglio
EEfficient numerical computation of the basicreproduction number for structured populations
Dimitri Breda , Francesco Florian , Jordi Ripoll and Rossana Vermiglio CDLab – Computational Dynamics LaboratoryDepartment of Mathematics, Computer Science and Physics – University of Udinevia delle scienze 206, 33100 Udine, [email protected] Institute of Mathematics–University of ZürichWinterthurerstrasse 190, CH-8057 Zürich, Switzerlandand CDLab – Computational Dynamics Laboratoryfrancesco.fl[email protected] Department of Computer Science, Applied Mathematics and Statistics – University of GironaCampus Montilivi, 17003 Girona, [email protected] CDLab – Computational Dynamics LaboratoryDepartment of Mathematics, Computer Science and Physics – University of Udinevia delle scienze 206, 33100 Udine, [email protected]
Abstract
As widely known, the basic reproduction number plays a key role in weighingbirth/infection and death/recovery processes in several models of population dy-namics. In this general setting, its characterization as the spectral radius of nextgeneration operators is rather elegant, but simultaneously poses serious obstaclesto its practical determination. In this work we address the problem numericallyby reducing the relevant operators to matrices through a pseudospectral collo-cation, eventually computing the sought quantity by solving finite-dimensionaleigenvalue problems. The approach is illustrated for two classes of models, re-spectively from ecology and epidemiology. Several numerical tests demonstrateexperimentally important features of the method, like fast convergence and influ-ence of the smoothness of the models’ coefficients. Examples of robust analysis ofinstances of specific models are also presented to show potentialities and ease ofapplication.
Keywords: structured population dynamics, stability analysis of equilibria, nextgeneration operator, spectral radius, spectral approximation, pseudospectral col-location
MSC: a r X i v : . [ m a t h . NA ] A p r Introduction
Aim of this paper is to give a concise, efficient and effective method for the numericalcomputation of the basic reproduction number , a fundamental tool in the analysis ofecological and epidemiological models of population dynamics, as well as in tacklingevolutionary aspects in those populations, see e.g. [4, Section 2.1] and the referencestherein. We refer to [17] as a starting reference on this quantity, commonly denoted R . The need for a numerical approach is soon justified since we are interested in popu-lation models posed on infinite-dimensional spaces, where R is usually characterizedas 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 readily practicablefor 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 are obtained by computing the (dominant) eigenvalues ofthe resulting matrices.As far as this general strategy is concerned, to the best of our knowledge [25] is thefirst scheme appearing in the literature. Therein an elementary discretization basedon the celebrated Euler’s method is proposed to approximate the basic reproductionnumber of a class of age-structured epidemics. Very recently a slight improvementhas been presented in [16], obtained by employing θ -methods, always in the contextof age-structured epidemics. Made exception for [14], which is anyway related to theapproach proposed here (see below), [16, 25] are the only two works concerning theproper numerical computation of R in continuously structured populations.Here, on the one hand, we intend to largely improve the seminal idea of [25] by ap-plying 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 convergenceideally of infinite order (known as spectral accuracy , see, e.g., [29]). This represents agreat advantage compared to the finite-order convergence of [16, 25], as it translatesinto much more accurate approximations obtained with much smaller matrices. Asa consequence the required computational load reduces, in terms of both time andmemory usage. The latter is a favorable feature when stability and bifurcation analy-ses are the final target in presence of varying or uncertain model parameters, and ofcourse 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, weillustrate two paradigmatic classes of models, namely spatially-structured cell popu-lations (whose prototype representative is referred to as model A in the sequel) andinfectious diseases in age-structured populations (whose prototype representative isreferred to as model B in the sequel). This choice potentially serves as a basis for ex-tending the proposed method to more involved systems of population dynamics (e.g.,with several structuring variables). If a structured epidemic model describing thenovel pandemia of coronavirus was available, such as Model B or any of its extensions,2his method could be applied to compute R for different values of the parameters.The contents of this research are organized as follows. In Section 2 we present thetwo classes of models mentioned above, introducing all the necessary ingredients andthe operators leading to the definition of R . The numerical treatment is proposedin Section 3, where explicit expressions for the discretizing matrices are constructivelygiven, also for the sake of implementation for those interested (Matlab codes are freelyavailable at http://cdlab.uniud.it/software ). In Section 4 we first illustrate thenumerical properties of the proposed techniques, analyzing experimentally error andconvergence, and then we finally perform some robust analysis on specific instancesof the models of interest.Let us remark that here we deliberately quit to tackle a theoretical investigationof the convergence features of the proposed methods. This would indeed requiresuitable tools, falling out of the scope of the current experimental study, as well asenough room for a proper description. We thus reserve to give formal proofs as wellas possible improvements in a related work, currently in preparation [8]. As far asother developments are concerned, we also plan to extend this approach in order totreat 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], wherean embryonic study of collocation techniques for the numerical computation of R isproposed (also the work [14] mentioned above is inspired from [13]). Moreover, let usalso note that pseudospectral techniques have become a reference tool in the numericaltreatment of infinite-dimensional dynamical systems in view of either stability andbifurcation analyses in the context of population dynamics, see, e.g., [6] or [7] and thereferences therein. Continuously structured models of population dynamics can be described by nonlin-ear abstract differential equations for the density of individuals with respect to somestructuring variables, e.g., age, size or space.With the aim at setting the theoretical background, let X be a Banach lattice andlet u ∈ X represent the density of individuals. We introduce next a continuous-timeevolution equation for the population density. For the theory of population dynamicswe 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 thesequel), then in general a nonlinear problem can be decomposed into birth and non-birth terms as u (cid:48) ( t ) = B ( u ( t )) u ( t ) − M ( u ( t )) u ( t ) , u ( t ) ∈ X , t ≥
0, (2.1)3here, for any u ∈ X , B ( u ) is the birth linear operator and M ( u ) is the linear operatoraccounting for other processes non related to birth (e.g., mortality or transition). Ofcourse, this decomposition depends on what is exactly considered as a birth event, seean example in Section 2.1 and also [3, 9].If we are concerned with the extinction of the population in (2.1), then we canfocus on the steady-state u ∗ ≡
0. Its stability is analyzed by means of the (formally)linearized model u (cid:48) ( t ) = B ( ) u ( t ) − M ( ) u ( t ) , u ( t ) ∈ X , t ≥
0. (2.2)On the other hand, if we focus on epidemiological models (the class of models B inthe 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 + u (cid:48) ( t ) = B ( v ( t ) , u ( t )) u ( t ) − M ( v ( t ) , u ( t )) u ( t ) , v ( t ) ∈ X m , u ( t ) ∈ X , t ≥
0. (2.3)The remaining ones have the general form v (cid:48) ( t ) = F ( v ( t ) , u ( t )) , v ( t ) ∈ X m , u ( t ) ∈ X , t ≥
0, (2.4)with F depending on the type of epidemics. For any v ∈ X m and u ∈ X , B ( v , u ) is thelinear infection operator and M ( v , u ) is the linear operator describing processes otherthan 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 ini-tial epidemic spread in (2.3)-(2.4). We can thus focus on the disease-free steady state ( v ∗ , u ∗ ) ∈ X m + × { } (where X + stands for the positive cone in X ). Since the (formally)linearized system has triangular form, under common assumptions, the stability anal-ysis of the disease-free equilibrium is reduced to the linear model u (cid:48) ( t ) = B ( v ∗ , 0 ) u ( t ) − M ( v ∗ , 0 ) u ( t ) , u ( t ) ∈ X , t ≥
0. (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 equationof the form u (cid:48) ( t ) = Bu ( t ) − Mu ( t ) , u ( t ) ∈ X , t ≥
0. (2.6)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 in-fection) and M : D ( M ) ⊆ X → X is a linear operator accounting for the remainingprocesses, which we call mortality for brevity (thus meant as proper mortality, or re-covery or, in general, any stage transition as in [27, Part 2]). Typically, the domain4 ( M ) of M is a subset of a subspace Y ⊆ X where some degree of smoothness ispresent, characterized by additional constraints in general described through a linearmap C : X → R p for some positive integer p : D ( M ) = { φ ∈ Y : C φ = } . (2.7)Other common (biologically meaningful) assumptions are required. In particular,as in [4], B is positive and bounded, while − M is the infinitesimal generator ofa strongly-continuous semigroup { T ( t ) } t ≥ of positive linear operators. Its spectralbound 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 M − = (cid:82) ∞ T ( t ) d t . 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 popula-tion 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, seethe discussion in [4]. However, in this standard case we can define the next genera-tion operator as the operator BM − : X → X acting on the same (whole) space of thepopulation density and, eventually, characterize the basic reproduction number as itsspectral radius: R : = ρ ( BM − ) . (2.8)Let us remark that the next generation operator is a generalization to the infinite-dimensional setting of the well-known next generation matrix for finite-dimensionalmodels (i.e., dim X < + ∞ ). Moreover, although the theoretical framework for thebasic 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 efficientcomputation as already anticipated in the introduction.It follows from definition (2.8) that R is a non-negative spectral value being BM − positive and bounded, see, e.g., [26]. It is actually the largest λ ≥ 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 B φ = λ M φ (2.9)for some nontrivial generalized eigenfunction φ ∈ D ( M ) . In particular, if BM − is alsocompact, then the Krein-Rutman theorem [24] ensures that R is a positive eigenvalue(for slightly weaker, yet involved assumptions see [21, Section 10.2]). At this point letus anticipate that the numerical method proposed in Section 3 relies on this assump-tion, so that (2.9) can be reduced to a standard (read finite-dimensional) generalizedeigenvalue problem for matrices. Compactness is proved in [8] under reasonable as-sumptions on the specific ingredients of models A and B as described in the forthcom-ing Sections 2.1 and 2.2. A discussion on further aspects related to R being or not ageneralized eigenvalue and their relevance on the numerical outcome is left to Section4. 5ollowing definition (2.8), with any further assumption, we obtain the useful upperbound R ≤ (cid:107) BM − (cid:107) X ← X : = sup ψ ∈ X \{ } (cid:107) BM − ψ (cid:107) X (cid:107) ψ (cid:107) X = sup φ ∈D ( M ) \{ } (cid:107) B φ (cid:107) X (cid:107) M φ (cid:107) X . (2.10)If, in addition, R = λ ≥ φ ∈ D ( M ) \ { } , then R = (cid:107) B φ (cid:107) X (cid:107) M φ (cid:107) X . (2.11)Finally, let us also recall that R is an alternative to the Malthusian parameter, i.e.,the exponential population growth rate. Concerning computations, in general it canbe more difficult to deal with the latter, i.e., the spectral bound r : = s ( B − M ) of thecomplete generator, rather than with the basic reproduction number R . For instance,we cannot ensure that r is a spectral value (not even an eigenvalue) and there can be noupper bound for r as it is always the case for R (corresponding to a sufficient conditionfor population extinction or infection eradication). Moreover, the rank of B is typicallylower than that of B − M , which is an advantage for R . Finally, extra evolutionaryaspects added to the models are more often related to R than to r . Summarizing,besides the well-known sign relation between these quantities, see, e.g., [28], the basicreproduction number offers diverse advantages over the Malthusian parameter, fromboth 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 intestineof an animal host, see, e.g., [3, 4]. Given that the intestine is portrayed as the linesegment [ l ] of finite length l >
0, let us set X : = L ( l ) and let u ( · , t ) ∈ X be thespatial density of bacteria along the intestine at time t ≥
0. The nonlinear problem isdescribed as ∂ t u ( x , t ) + ∂ x (cid:2) c ( x , S ( t )) u ( x , t ) − D ( x , S ( t )) ∂ x u ( x , t ) (cid:3) + µ ( x , S ( t )) u ( x , t ) = β ( x , S ( t )) u ( x , t ) c ( S ( t )) u ( t ) − D ( S ( t )) ∂ x u ( t ) = c ( l , S ( t )) u ( l , t ) − D ( l , S ( t )) ∂ x u ( l , t ) =
0, (2.12)where S ( t ) : = (cid:82) l σ ( x ) u ( x , t ) d x , t ≥
0, is the population size weighted through a non-negative 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 theposition x ∈ [ l ] and of the size S . Note that transport, diffusion and vital processesare density-dependent accounting for limited resources (crowding effects) and space-specific accounting for the heterogeneity of the intestine, see [3, 4, 20].6or a population of bacteria, one can consider either symmetric division, i.e., whena mother cell divides then two daughter cells are born and the former disappears, orasymmetric division, i.e., one mother and one daughter. Both types can be consideredsimultaneously in (2.12) by setting β − µ = (cid:2) νβ + ( − ν ) β (cid:3) − (cid:2) ( − ν ) β + µ (cid:3) forsome ν ∈ [
0, 1 ] representing the probability of asymmetric division. For instance, forspace-independent vital rates one would readily get that R = ( − ν ) β / [( − ν ) β + µ ] since integrating over the whole space in (2.12), the system reduces to an ODE forthe population size, see [4]. The latter is a clear example of the fact that the basicreproduction number depends on what is exactly considered as a birth event. In anycase, here we assume symmetric division ( ν =
0) for simplicity, and without loss ofgenerality 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 ∗ ≡ ∂ t u ( x , t ) + ∂ x (cid:2) c ( x ) u ( x , t ) − D ( x ) ∂ x u ( x , t ) (cid:3) +[ β ( x ) + µ ( x )] u ( x , t ) = β ( x ) u ( x , t ) c ( ) u ( t ) − D ( ) ∂ x u ( t ) = c ( l ) u ( l , t ) − D ( l ) ∂ x u ( l , t ) =
0, (2.13)where we implicitly set c ( x ) : = c ( x , 0 ) , D ( x ) : = D ( x , 0 ) , β ( x ) : = β ( x , 0 ) and µ ( x ) : = µ ( x , 0 ) to lighten the notation though with a little abuse. Moreover, we further assume c ( x ) > x ∈ ( l ) , as well as β and µ bounded with their sum not identicallyvanishing. Concerning diffusion, we consider either D ( x ) ≥ ˜ D , x ∈ [ l ] , for somepositive ˜ D (everywhere positive diffusion) or D ≡ B : X → X as ( B φ )( x ) : = β ( x ) φ ( x ) , x ∈ [ l ] , (2.14)and the mortality operator M : D ( M ) ⊆ X → X as ( M φ )( x ) : = (cid:2) c ( x ) φ ( x ) − D ( x ) φ (cid:48) ( x ) (cid:3) (cid:48) + [ β ( x ) + µ ( x )] φ ( x ) , x ∈ [ l ] , (2.15)with domain D ( M ) : = (cid:8) φ ∈ X : φ (cid:48) , ( c φ − D φ (cid:48) ) (cid:48) ∈ X and c ( x ) φ ( x ) − D ( x ) φ (cid:48) ( x ) = x = l (cid:9) . (2.16)Then, following (2.10), the upper bound R ≤ sup φ ∈D ( M ) \{ } (cid:82) l β ( x ) | φ ( x ) | d x (cid:82) l [ β ( x ) + µ ( x )] | φ ( x ) | d x ≤ .2 Model B: age-structured epidemics 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 splitup the individuals of the population into three classes according to the stage of theinfectious disease. If l > X : = L ( 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 ∂ t s ( x , t ) + ∂ x s ( x , t ) + µ ( x ) s ( x , t ) = − f ( x , t ) s ( x , t ) + δ ( x ) i ( x , t ) + σ ( x ) r ( x , t ) ∂ t i ( x , t ) + ∂ x i ( x , t ) + µ ( x ) i ( x , t ) = f ( x , t ) s ( x , t ) − [ γ ( x ) + δ ( x )] i ( x , t ) ∂ t r ( x , t ) + ∂ x r ( x , t ) + µ ( x ) r ( x , t ) = γ ( x ) i ( x , t ) − σ ( x ) r ( x , t ) s ( t ) = (cid:82) l β ( x )[ s ( x , t ) + ( − θ ) i ( x , t ) + ( − ˆ θ ) r ( x , t )] d xi ( t ) = θ (cid:82) l β ( x ) i ( x , t ) d xr ( t ) = ˆ θ (cid:82) l β ( x ) r ( x , t )) d x , (2.18)where f is the force of infection defined as f ( x , t ) : = (cid:82) l K ( x , y , N ( t )) i ( y , t ) d y , x ∈ [ l ] , t ≥
0, with N ( t ) the total population at time t and infection kernel K typicallydecomposed as K ( x , y , N ) = χ ( y ) · c ( x , y , N ) / N , x , y ∈ [ l ] , N >
0. Above, χ isthe probability of infection transmission through an infectious contact and c is thedensity-dependent contact rate of a susceptible individual of age x with an infectedindividual of age y . Moreover, β , µ , γ , δ and σ are, respectively, the fertility, (natural)mortality, removal, recovery and loss of immunity rates, all non-negative and with β not identically vanishing. Finally, θ and ˆ θ are the probabilities of vertical transmissionof infectiveness and immunity, respectively. Note that recovery, removal, immunityloss and vital processes are all age-specific during the lifespan [ l ] .System (2.18) includes different types of epidemic models. For instance, the SISmodel is given by γ ≡ σ ≡
0; the SIR model is given by δ ≡ σ ≡
0; the SIRS modelis given by δ ≡
0. The exact interpretation of the removed class depends on thetype of model considered. Moreover, we could also consider a vaccination rate andthe resulting system would have an extra transition, from susceptible to recoveredindividuals in this case, or the consideration of exposed individuals, or even we couldconsider a multistrain epidemic model. See [23] and [18, Chapter VI], [21, Chapters6–8], [27, Chapter 22].From now on, we assume that the population has reached a demographic steadystate. Indeed, we assume zero demographic growth, i.e., (cid:90) l β ( x ) Π ( x ) d x =
1, (2.19)where Π ( x ) : = exp (cid:0) − (cid:82) x µ ( y ) d y (cid:1) , x ∈ [ l ] , together with Π ( l ) =
0, defines thesurvival probability, as well as initial conditions such that s ( x ) + i ( x ) + r ( x ) = Π ( x ) / (cid:82) l Π ( y ) d y , x ∈ [ l ] , for some positive N . Therefore, the total population N ( t ) = (cid:82) l [ s ( x , t ) + i ( x , t ) + r ( x , t )] d x = N , t ≥
0, remains constant over time.In what follows we are concerned with the early stage of the epidemics of aninfectious disease of any type (SIS, SIR or SIRS). Linearizing (2.18) around the disease-free equilibrium s ∗ ( x ) = N Π ( x ) / (cid:82) l Π ( y ) d y , i ∗ ≡ r ∗ ≡ ∂ t i ( x , t ) + ∂ x i ( x , t ) = (cid:82) l χ ( y ) · c ( x , y , N ) N · i ( y , t ) d y · N Π ( x ) (cid:82) l Π ( y ) d y − [ µ ( x ) + γ ( x ) + δ ( x )] i ( x , t ) i ( t ) = θ (cid:82) l β ( x ) i ( x , t ) d x for the infected individuals. We let (cid:82) l µ ( x ) d x = + ∞ to avoid immortal individuals. Inorder to properly deal with this condition, we make the change of variables i ( x , t ) = u ( x , t ) Π ( x ) , thus removing the mortality term and reducing the problem to (cid:40) ∂ t u ( x , t ) + ∂ x u ( x , t ) = (cid:82) l K ( x , y ) u ( y , t ) d y − [ γ ( x ) + δ ( x )] u ( x , t ) u ( t ) = θ (cid:82) l β ( x ) u ( x , t ) d x (2.20)with effective infection kernel K ( x , y ) : = χ ( y ) c ( x , y , N ) · Π ( y ) (cid:82) l Π ( z ) d z ≥ x , y ∈ [ l ] , (2.21)and effective fertility rate β ( x ) : = β ( x ) Π ( x ) ≥ x ∈ [ l ] . We recall that (cid:82) l β ( x ) d x = u ( l , t ) can be strictly positive allowing for the possibilityof chronic infected individuals.Eventually, with respect to the reference linear equation (2.6), from the normalizedmodel (2.20) it is natural to define the birth operator B : X → X as ( B φ )( x ) : = (cid:90) l K ( x , y ) φ ( y ) d y , x ∈ [ l ] , (2.22)and the mortality operator M : D ( M ) ⊆ X → X as ( M φ )( x ) : = φ (cid:48) ( x ) + [ γ ( x ) + δ ( x )] φ ( x ) , x ∈ [ l ] , (2.23)with domain D ( M ) : = (cid:26) φ ∈ X : φ (cid:48) ∈ X and φ ( ) = θ (cid:90) l β ( x ) φ ( x ) d x (cid:27) . (2.24) We recall that it is actually the infection operator. We recall that it is actually the recovery and transition operator. BM − (and even obtain explicit expressions for some special cases, see Section 4).Indeed, by defining Π ( x ) : = exp (cid:0) − (cid:82) x [ γ ( y ) + δ ( y )] d y (cid:1) ≤
1, the next generationoperator becomes ( BM − ψ )( x ) = (cid:90) l K ( x , y ) Π ( y ) (cid:18) C + (cid:90) y ψ ( z ) Π ( z ) d z (cid:19) d y , x ∈ [ l ] , (2.25)with C : = θ (cid:90) l β ( x ) (cid:18) (cid:90) x Π ( x ) Π ( z ) ψ ( z ) d z (cid:19) d x · (cid:18) − θ (cid:90) l β ( x ) Π ( x ) d x (cid:19) − ≥ θ Π (cid:54)≡ R ≤ (cid:90) l (cid:90) l K ( x , y ) (cid:18) θ Π ( y ) − θ (cid:82) l β ( x ) Π ( x ) d x + (cid:19) d y d x , (2.26)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 forthe disease eradication. Let N be a positive integer. Under the assumption that R is a generalized eigenvalue,we propose a general discretization approach consisting in the collocation of (2.9) ona mesh of N + = : x N ,0 < x N ,1 < · · · < x N , N : = l distributed on [ l ] , l >
0. Let us soon observe that collocation is meaningless in the whole space X = L ( 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 + be the finite-dimensional counterpart of X and let Φ : = ( Φ , Φ , . . . , Φ 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 = x N , N = l in the collocation mesh is justifiedby the fact that, in general, φ belongs to some restricted domain (2.7) where the mor-tality 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, thechoice above allows for a simplified treatment as one can directly relate the boundaryconditions to the collocation equations.As the concerned mortality operators usually involve differentiation and/or inte-gration, see (2.15) and (2.23), let us also introduce the differentiation matrix and the10 uadrature weights associated to the collocation nodes. The first, denoted H N , has en-tries h N ; i , j : = (cid:96) (cid:48) N , j ( x N , i ) , i , j =
0, 1, . . . , N , where { (cid:96) N ,0 , (cid:96) N ,1 , . . . , (cid:96) N , N } is the Lagrange basis relevant to the chosen nodes: if f isa smooth function on [ l ] and v : = ( f ( x N ,0 ) , f ( x N ,1 ) , . . . , f ( x N , N )) T , then H N v is anapproximation to ( f (cid:48) ( x N ,0 ) , f (cid:48) ( x N ,1 ) , . . . , f (cid:48) ( x N , N )) T . The second, components of thevector w N : = ( w N ,0 , w N ,1 , . . . , w N , N ) T ∈ R N + , are given by w N , j : = (cid:90) l (cid:96) N , j ( x ) d x , j =
0, 1, . . . , N ,and, for the same v above, w TN v is an approximation to (cid:82) l f ( x ) d x . Both follow straight-forwardly from approximating f with the N -degree interpolating polynomial p N ( x ) : = N ∑ j = (cid:96) N , j ( x ) f ( x N , j ) , x ∈ [ l ] , (3.1)which indeed satisfies p ( x N , i ) = f ( x N , i ) , i =
0, 1, . . . , N . Even though at this momentit is not necessary to specify any particular choice of nodes, let us remark that in thecase of Chebyshev extremal points, as we assume in Section 4, H N and w N can beobtained rather easily . 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 B N Φ = λ M N Φ (3.2)of the generalized eigenvalue problem (2.9), where the structure of the matrix rep-resentation of the finite-dimensional operators B N and M N depends on the specificmodel as detailed in the following Section 3.1 for model A and Section 3.2 for modelB. 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 shownto be the case in Section 4 by way of several numerical experiments, properly tunedto test the error behavior, to measure the convergence rate and to detect other pecu-liarities of the proposed approach. Note that the eigenvalues of (3.2) can be computedwith 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 remarkthat 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. Here a N ; i , j denotes the ( i , j ) -th entry of a matrix A ∈ R ( N + ) × ( N + ) , i , j =
0, 1, . . . , N . This notationis also extended later on according to Matlab’s colon notation: for instance a N ; i ,0: N and a N ;0: N , j indicate,respectively, the i -th row and the j -th column of A . Usually these quantities are provided with reference to the normalized interval [ −
1, 1 ] . R is based on the assumption that this number actually corresponds to a generalizedeigenvalue. If this is not the case, already (2.9) looses sense. Nevertheless, in Section4, we report on some tests under the latter hypothesis, where the scheme is still ableto give reasonable approximations (even though with slower convergence than whatexperimented in the case of generalized eigenvalues). Let us recall the main ingredients (2.14), (2.15) and (2.16) of the class of cell populationmodels, model A in Section 2.1. The generalized eigenvalue problem (2.9) reads ( B φ )( x ) = λ ( M φ )( x ) , x ∈ [ l ] , (3.3)with generalized eigenfunction φ ∈ D ( M ) \ { } , i.e., c ( x ) φ ( x ) − D ( x ) φ (cid:48) ( x ) = x = l . (3.4)The proposed collocation consists in looking for an N -degree polynomial p N satisfying ( Bp N )( x N , i ) = λ ( Mp N )( x N , i ) , i =
1, . . . , N − c ( x N , i ) p N ( x N , i ) − D ( x N , i ) p (cid:48) N ( x N , i ) = i = N . (3.5)By using the Lagrange representation p N ( x ) = N ∑ j = (cid:96) N , j ( x ) Φ j , x ∈ [ l ] , (3.6)and by recalling the differentiation matrix H N introduced in Section 3, it is not difficultto recover (3.2) with B N ∈ R ( N + ) × ( N + ) given as B N ; i , j : = (cid:40) β ( x N , i ) , i = j =
1, . . . , N − i = N or i (cid:54) = j ,and M N ∈ R ( N + ) × ( N + ) given as M N ; i , j : = C N ;0,0 δ j − D N ;0,0 h N ;0, j , i = j =
0, 1, . . . , N , [ H N ( C N − D N H N ) + Σ N ] i , j , i =
1, . . . , N , j =
0, 1, . . . , N , C N ; N , N δ N , j − D N ; N , N h N ; N , j , i = N , j =
0, 1, . . . , N ,where δ i , j is the Kronecker’s delta and C N , D N , Σ N ∈ R ( N + ) × ( N + ) are defined as C N ; i , j : = (cid:40) c ( x N , i ) , i = j =
0, 1, . . . , N ,0, i (cid:54) = j ,12 N ; i , j : = (cid:40) D ( x N , i ) , i = j =
0, 1, . . . , N ,0, i (cid:54) = j , Σ N ; i , j : = (cid:40) β ( x N , i ) + µ ( x N , i ) , i = j =
0, 1, . . . , N ,0, i (cid:54) = j .Note how the first and last rows of B N and M N realize the discrete boundary condi-tions (3.5), simulating those characterizing D ( M ) in (2.16), i.e., (3.4). Let us recall the main ingredients (2.22), (2.23) and (2.24) of the class of age-structuredepidemics, model B in Section 2.2. The generalized eigenvalue problem (2.9) readsagain as (3.3), but with generalized eigenfunction φ ∈ D ( M ) \ { } , i.e., satisfying φ ( ) = θ (cid:90) l β ( x ) φ ( x ) d x . (3.7)The proposed collocation consists in looking for an N -degree polynomial p N satisfying ( Bp N )( x N , i ) = λ ( Mp N )( x N , i ) , i =
1, . . . , N ,together with p N ( ) = θ (cid:90) l β ( x ) p N ( x ) d x , (3.8)approximating the integral above as well as the one defining B in (2.22) by quadra-ture. By using again the Lagrange representation (3.6) and by recalling the (vector of)quadrature weights w N introduced in Section 3, it is not difficult to recover (3.2) with B N ∈ R ( N + ) × ( N + ) given as B N ; i , j : = (cid:40) i = j =
0, 1, . . . , N , w N , j K ( x N , i , x N , j ) , i =
1, . . . , N , j =
0, 1, . . . , N ,and M N ∈ R ( N + ) × ( N + ) given as M N ; i , j : = (cid:40) δ j − θ w N , j β ( x N , j ) , i = j =
0, 1, . . . , N , h N ; i , j + [ γ ( x N , i ) + δ ( x N , i )] δ i , j , i =
1, . . . , N , j =
0, 1, . . . , N ,where δ i , j is again the Kronecker’s delta. Note how the first row of both B N and M N realizes the (quadrature of the) discrete boundary condition (3.8) simulating the onecharacterizing D ( M ) in (2.24), i.e., (3.7). 13 Results
A first series of numerical tests is presented in Section 4.1 with the aim of illustratingthe convergence properties of the proposed techniques, as well as related aspects suchas, e.g., the effect of the smoothness of the models’ coefficients. Then in Section 4.2 weuse the tested algorithms to perform some quantitative studies in the context of vary-ing 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 16GB memory, through codes written in Matlab version R2019a (freely available at http://cdlab.uniud.it/software as anticipated in the Introduction). In these codeswe always use Chebyshev extremal points as collocation nodes, i.e., x N , i = l [ − cos ( i π / N )] /2, i =
0, 1, . . . , N . We refer to [29] for their numerous advantageous prop-erties in the context of numerical interpolation, differentiation and quadrature. The numerical properties of the proposed approach are mainly investigated by ana-lyzing the behavior of the error | R N − R | between the approximated and the exactbasic reproduction numbers. Here for exact we mean either the theoretical value of R when this is explicitly available, or a reference value R
0, ¯ N computed with a suitablylarge ¯ N when an exact value is not attainable.A second error is also presented when possible, viz. (cid:107) p N − φ (cid:107) ∞ , M . This is theerror between the exact generalized eigenfunction (if available) and its collocation ap-proximation, measured as the maximum absolute value of their differences on a meshof M equidistant points in [ l ] . In particular, we always use M = p N is reconstructed from the computed generalized eigenvector Φ associated to the dominant generalized eigenvalue of (3.2) through barycentric inter-polation by applying the algorithm proposed in [5]. Of course, the same normalizationis prescribed from time to time for both p N and φ .Finally, the discrete generalized eigenvalue problem (3.2) is solved either by Mat-lab’s eig or eigs . We also test the discretization B N M − N Ψ = λ Ψ , (4.1)of the eigenvalue problem BM − ψ = λψ (4.2)equivalent to (2.9), where ψ = M φ for φ in (2.9) and Ψ = M N Φ for Φ in (3.2). Theseeigenproblems are implemented as eig(B,M) for (3.2) and eig(B/M) for (4.1) (or thealternative versions through eigs ). As far as the models to test are concerned, below we list all the specific choices bygiving the defining rates and coefficients. In parallel we also give some of their keyfeatures. 14or all the instances of model A we set l = c ( x ) = l − x , x ∈ [ l ] , and D ( x ) : = ˜ D · (cid:20) l x ( l − x ) + (cid:21) ≥ ˜ D , x ∈ [ l ] .The other ingredients differ as described next.(A1) [the immortal case] Let ˜ D = µ ≡ β ( x ) : = ˜ β · (cid:20) l x ( l − x ) + (cid:21) , x ∈ [ l ] ,for ˜ β =
1. With these choices the next generation operator is compact accordingto [8] and thus R is a generalized eigenvalue. It is not difficult to recover R = φ ( x ) = e (cid:82) x c ( y ) D ( y ) d y , x ∈ [ l ] , (4.3)normalized as φ ( ) =
1, and the eigenfunction in (4.2) is ψ = βφ .(A2) [the proportional case] Let ˜ D = µ ( x ) : = ˜ µ xl , β ( x ) : = ˜ βµ ( x ) , x ∈ [ l ] ,for ˜ µ = β = R = β / ( ˜ β + ) = ψ = ( ˜ β + ) µφ .(A3) [the general case] Let µ be the same as in case (A2) and β be expressed as in case(A1), but with ˜ β =
10. We consider either(A3.1) the compact case: ˜ D = (A3.2) a “almost non-compact” case: ˜ D = − ;(A3.3) the non-compact case: ˜ D = Independently of these choices, R is unknown, as well as the correspondinggeneralized eigenfunction.For all the instances of model B we set again l = Here we refer to the compact case since when there is diffusion everywhere in space and undertechnical assumptions on the model parameters, the next generation operator is compact according to[8]. For the case of lack of diffusion R can be computed analytically and it turns out to be not aneigenvalue. age-independent case] Let µ ≡ ˜ µ =
28 (so that e − ˜ µ l (cid:39) β ≡ µ (so that (2.19)holds), γ ≡ ˜ γ = δ ≡ ˜ δ = K ( x , y ) : = ˜ k ˜ µ e − ˜ µ y , x , y ∈ [ l ] ,for ˜ k =
52. It also follows that β ( x ) = ˜ µ e − ˜ µ x , Π ( x ) = e − ( ˜ γ + ˜ δ ) x , x ∈ [ l ] .Finally, let θ = R (cid:39) ˜ k / [( − θ ) ˜ µ + ˜ γ + ˜ δ ] =
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 canrecover through (2.23) the generalized eigenfunctions φ ( x ) = e − x (cid:18) φ ( ) − ψ (cid:19) + ψ x ∈ [ l ] , (4.4)with φ ( ) (cid:39) ψ /182 following from (3.7), with ψ playing the role of normalizingfactor.(B2) [the proportionate case] Let K ( x , y ) : = ˜ kx ( l − x ) · Π ( y ) (cid:82) l Π ( z ) d z , x , y ∈ [ l ] ,for ˜ k =
16 065/64 and Π ( x ) : = (cid:18) l − xl (cid:19) α , x ∈ [ l ] , (4.5)for some α >
0. Note that (cid:82) l Π ( z ) d z = l α + . Let moreover γ and δ satisfy γ ( x ) + δ ( x ) = l − x , x ∈ [ l ] , (4.6)so that also Π ( x ) = l − xl , x ∈ [ l ] , follows. We assume absence of verticaltransmission, i.e., θ =
0, which makes useless to specify β . With these choicesthe next generation operator is again compact according to [8]. Since K aboveis the product of functions in each of the variables x and y , the next generationoperator in (2.25) becomes a rank one operator and we can explicitly compute R = k ( α + ) l ( α + )( α + )( α + )( α + ) .As far as α is concerned, we consider For this specific choice, and for the sake of pragmatism, here we exclude the last node x N , N = l from the computations to avoid overflow. In [8] the issue is tackled more rigorously from the numericalstandpoint. α = l and R = α =
1: life expectancy is 0.5 l and R = α , the eigenfunction corresponding to R is ψ ( x ) = x ( l − x ) , normalized as ψ ( l /2 ) = l /16, and the generalized eigenfunction is φ ( x ) = Π ( x ) (cid:90) x y ( l − y ) Π ( y ) d y , x ∈ [ l ] ,which, for the choices above, turns out to be a polynomial of degree 5, viz. φ ( x ) = x ( l − x )( l − x )
12 , x ∈ [ l ] (4.7)(normalized as φ ( l /2 ) = l /384).(B3) [the general case] Let K ( x , y ) : = ˜ ke − (cid:16) x − yl (cid:17) · Π ( y ) (cid:82) l Π ( z ) d z , x , y ∈ [ l ] ,for ˜ k = l = l and Π ( x ) : = e − α xl − x , x ∈ [ l ] , for α = l ). Some calculations give (cid:90) l Π ( x ) d x = l (cid:18) − α e α (cid:90) + ∞ α e − x x d x (cid:19) Let moreover γ and δ be as in case (B1). Here, infection events occur mostlybetween individuals of the same age and there are chronic infected individualssince e − ( ˜ γ + ˜ δ ) l >
0. Finally, let θ = β ( x ) : = b ( x ) Π ( x ) / (cid:82) l b ( y ) Π ( y ) d y , x ∈ [ l ] , for b ( x ) : = ( x / l ) e − x / l , which ensures (2.19) and maximum fertilityat x = l /3. For these choices the next generation operator is compact , but both R and the corresponding generalized eigenfunction are unknown. As a first experiment we test the convergence properties of the proposed approachfor 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 approximationof both R and the relevant generalized eigenfunction φ . To note that for the formerthe 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 trendsemerge also in the right panel, where (4.2) is solved through eig(B/M) instead of solv-ing (2.9) through eig(B,M) as in the left panel. Let us just remark that the choice of In general, under technical assumptions on the infection kernel, the next generation operator turnsout to be compact according to [8]. ig(B/M) seems slightly better from the point of view of algorithmic stability. Finally,let us mention that we have compared also with the use of eigs , obtaining indistin-guishable results. As these outcomes are unchanged in the remaining experiments, weshow only the results relevant to the use of eig(B/M) , although we measure alwaysthe error on the generalized eigenfunction φ . Figure 4.2 concerns the same analysis -16 -12 -8 -4 -16 -12 -8 -4 Figure 4.1: case (A1) – errors for increasing N on R = φ in (4.3), computedas (2.9) through eig(B,M) (left) and as (4.2) through eig(B/M) (right).above, but for case (A2): the results are qualitatively the same. -16 -12 -8 -4 -16 -12 -8 -4 Figure 4.2: case (A2) – errors for increasing N on R = φ in (4.3), computedas (2.9) through eig(B,M) (left) and as (4.2) through eig(B/M) (right).In the previous tests both R and the relevant generalized eigenfunction are ana-lytically known. As the same is not possible for case (A3), we show in Figure 4.3 onlythe error on R with respect to a reference value R N as explained at the beginningof Section 4.1. As expected, the convergence is spectrally accurate for case (A3.1),in which ˜ D = D = − (cid:54) = N are necessary to startappreciating the spectral accuracy. It is indeed reasonable to expect that the value of˜ D affects the error constants, causing their increase as ˜ D →
0, still being the problem18ompact. When we deal instead with case (A3.3), in which the absence of diffusioncauses the loss of compactness, convergence still occurs, even though at a fixed rate(seemingly linear). The fact is somehow surprising (and certainly merits future inves-tigation), given that without compactness (2.9) may even become meaningless and weare thus using a finite-dimensional eigenvalue problem to approximate componentsof the spectrum possibly other than the point one. -10 -5 Figure 4.3: case (A3) – error for increasing N on R (cid:39) ◦ , case(A3.1)), R (cid:39) • , case (A3.2)), R (cid:39) (cid:3) , case(A3.3)) computed as (4.2) through eig(B/M) (reference values R
0, ¯ N computed with¯ N = As far as case (B1) is concerned, the results of Figure 4.4 show spectral accuracy in theapproximation of both R 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 de-pending on α . First of all, let us recall that the generalized eigenfunction (4.7) is apolynomial of degree 5, justifying the sudden drop of the error to machine precisionoccurring with N = N = N = R is concerned,instead, the same behavior is ensured if enough regularity is provided, as it is the casefor α = Π issmooth for α ≥
1, but for α < x = l . Indeed, thechoice α =
10 15 20 2530 4010 -12 -8 -4 Figure 4.4: case (B1) – errors for increasing N on R (cid:39) φ in (4.4), computedas (4.2) through eig(B/M) (reference value R
0, ¯ N = N = -12 -8 -4 -12 -8 -4 Figure 4.5: case (B2) – errors for increasing N on R and on φ in (4.7), computed as(4.2) through eig(B/M) for α = R =
2) and α = R = R nor the relevantgeneralized eigenfunction are known. To conclude, we give some examples of practical application of the proposed tech-niques to analyze the behavior of the chosen models in the presence of varying pa-rameters. Below, the choices of the discretization index N are instructed by the resultsdiscussed in the previous section. Moreover, the intervals of variation of the concerned20 -12 -8 -4 Figure 4.6: case (B3) – error for increasing N on R , computed as (4.2) through eig(B/M) (reference value R
0, ¯ N = N = P points. Corresponding val-ues of R 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 lattercase we adopt standard contouring algorithms, viz. Matlab’s contourf . Besides thegraphical output, we give also some indication of the overall computational time withrespect to the parameter(s) grid size P . Let us anticipate that all the tests concern thegeneral instances of either model A or B, namely cases (A3) or (B3) as described inSection 4.1.In Figure 4.7 we investigate for case (A3) how R varies as a function of eitherglobal fertility parameter ˜ β (left panel) or global diffusion parameter ˜ D (right panel).For the former ˜ µ = µ =
50 and ˜ β =
10. As expected, R grows with respect to the fertility from zeroup to the asymptote R =
2, since we assumed symmetric division (2.17), with largerdiffusion 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 intermediatediffusion maximizing the basic reproduction number. At this point, a naive approachis that a population adopting such strategy would persist and become an unbeatablepopulation. Computationally speaking, both figures are obtained by repeatedly com-puting R with respect to the varying parameter upon discretizing the latter on a meshof P = R = R ( ˜ D , ˜ β ) . The computation was performed on a gridof size 105 × R as a function of the verti-cal transmission probability θ , confirming the biological expectation that this functionis monotonically increasing as it is explicitly for case (B1). Notice that an increaseof the vertical transmission can lead to the spread of the infection. The grid size is The role of the diffusion from the evolutionary point of view is out of the scope of the present work. Figure 4.7: case (A3) – R as a function of ˜ β for varying ˜ D and ˜ µ = R N computed with N =
30 for case (A3.1) and N =
500 for case (A3.2)) and R as afunction of ˜ D for ˜ β =
10 and ˜ µ =
50 (right, values R N computed with N = Figure 4.8: case (A3) – level curves of R as a function of ( ˜ D , ˜ β ) for ˜ µ =
50 (values R N computed with N = P = Figure 4.9: case (B3) – R as a function of θ (values R N computed with N = {
0, 0.25, 0.5, 0.75, 1 } , but independently of θ we can appreciate how the proposed ap-proach approximate smooth curves as expected. Figure 4.10: case (B3) – normalized eigenfunctions (left) and generalized eigenfunc-tions (right) for several values of θ , computed with N = Acknowledgement:
DB and RV are members of INdAM Research group GNCS andare partially supported by the INdAM GNCS project “Approssimazione numerica diproblemi di evoluzione: aspetti deterministici e stocastici” (2018). JR is part of theCatalan research group 2017 SGR 1392 and has partially received support from theSpanish Ministry of Science and Innovation, reference MTM2017-84214-C2, and fromthe University of Girona, references MPC UdG 2016/047 and PONT2019/08.
References [1] A. Andò and D. Breda. Collocation techniques for structured populations mod-eled by delay equations. In M. Aguiar, C. Brauman, B. Kooi, A. Pugliese, N. Stol-lenwerk, and E. Venturino, editors,
Current Trends in Dynamical Systems in Biologyand Natural Sciences , volume 21 of
SEMA SIMAI series , pages 43–62. Springer, 2020.In print.[2] A. Andò, D. Breda, D. Liessi, S. Maset, F. Scarabel, and R. Vermiglio. 15 yearsor so of pseudospectral collocation methods for stability and bifurcation of delayequations. In G. Valmorbida, W. Michiels, and P. Pepe, editors,
Incorporatingconstraints on the Analysis of Delay and Distributed Parameter Systems , Adv. DelaysDyn. series. Springer. To appear.[3] C. Barril, A. Calsina, and J. Ripoll. On the reproduction number of a gut micro-biota model.
Bull. Math. Biol. , 79:2727–2746, 2017.234] C. Barril, A. Calsina, and J. Ripoll. A practical approach to R in continuous-timeecological models. Math. Meth. Appl. Sci. , 41(18):8432–8445, 2017.[5] J.-P. Berrut and L. N. Trefethen. Barycentric Lagrange interpolation.
SIAM Rev. ,46(3):501–517, 2004.[6] D. Breda, O. Diekmann, M. Gyllenberg, F. Scarabel, and R. Vermiglio. Pseu-dospectral discretization of nonlinear delay equations: new prospects for numer-ical bifurcation analysis.
SIAM J. Appl. Dyn. Sys. , 15(1):1–23, 2016.[7] D. Breda, S. Maset, and R. Vermiglio.
Stability of linear delay differential equations– A numerical approach with MATLAB . SpringerBriefs in Control, Automation andRobotics. Springer, New York, 2015.[8] D. Breda, J. Ripoll, and R. Vermiglio. Collocation of next-generation operators forcomputing the basic reproduction number of structured populations. In prepara-tion.[9] J. M. Cushing and O. Diekmann. The many guises of R (a didactic note). J. Theor.Biol. , 404:295–302, 2016.[10] O. Diekmann, J. A. P. Heesterbeek, and T. Britton.
Mathematical Tools for Un-derstanding Infectious Disease Dynamics . Theoretical and Computational Biology.Princeton University Press, Princeton and Oxford, 2013.[11] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition andthe computation of the basic reproduction number R in models for infectiousdiseases in heterogeneous populations. J. Math. Biol. , 28:365–382, 1990.[12] O. Diekmann, F. Scarabel, and R. Vermiglio. Pseudospectral discretization of de-lay differential equations in sun-star formulation: results and conjectures.
DiscreteContin. Dyn. S. - S , 8:95–105, 2007. To appear.[13] F. Florian. Numerical computation of the basic reproduction number in popula-tion dynamics. Master’s thesis, University of Udine, 2018.[14] F. Florian and R. Vermiglio. PC-based sensitivity analysis of the basic reproduc-tion number of population and epidemic models. In M. Aguiar, C. Brauman,B. Kooi, A. Pugliese, N. Stollenwerk, and E. Venturino, editors,
Current Trends inDynamical Systems in Biology and Natural Sciences , volume 21 of
SEPA SIMAI series ,pages 205–222. Springer, 2020. In print.[15] G. Golub and C. Van Loan.
Matrix computations . Johns Hopkins Studies in Mathe-matical Sciences. Johns Hopkins University Press, Baltimore, fourth edition, 2013.[16] W. Guo, M. Ye, X. Li, A. Meyer-Baese, and Q. Zhang. A theta-scheme approxi-mation of basic reproduction number for an age-structured epidemic system in afinite horizon.
Math. Biosci. Eng. , 16(5):4107–4121, 2019.2417] J. A. P. Heesterbeek. A brief history of R and a recipe for its calculation. ActaBiother. , 50:189–204, 2002.[18] M. Iannelli.
Mathematical Theory of Age-Structured Population Dynamics . AppliedMathematics Monographs (C.N.R.). Giardini Editori e Stampatori, Pisa, Italy,1995.[19] M. Iannelli and F. Milner.
The Basic Approach to Age-Structured Population Dynamics– Models, Methods and Numerics . Lecture Notes on Mathematical Modelling in theLife Sciences. Springer, The Netherlands, 2017.[20] M. Iannelli and A. Pugliese.
An Introduction to Mathematical Population Dynamics– Along the trail of Volterra and Lotka . Number 79 in La matematica per il 3+2.Springer, New York, 2014.[21] H. Inaba.
Age-Structured Population Dynamics in Demography and Epidemiology .Springer, New York, 2017.[22] H. Inaba. The basic reproduction number R in time-heterogeneous environ-ments. J. Math. Biol. , 79(2):731–764, 2019.[23] W. O. Kermack and A. G. McKendrick. Contributions to the mathematical theoryof epidemics III – Further studies of the problem of endemicity.
Proc. Roy. Soc.Lond. A , 141A:94–122, 1933. Reprinted in
Bull. Math. Biol. , 53(1/2), (1991), 89–18.[24] M. G. Krein and M. A. Rutman. Linear operators leaving invariant a cone in abanach space.
Uspehi Matem. Nauk (N. S.) 3 , 1(23):4–95, 1948. (in Russian). Amer.Math. Soc. Transl., 26:128pp, 1950 (in English).[25] T. Kuniya. Numerical approximation of the basic reproduction number for a classof age-structured epidemic models.
Appl. Math. Lett. , 73:106–112, 2017.[26] H. H. Shaefer.
Banach lattices and positive operators . Grundlehren der mathematis-chen Wissenschaften. Springer-Verlag, Berlin Heidelberg, 1974.[27] H. R. Thieme.
Mathematics in Population Biology . Theoretical and ComputationalBiology. Princeton University Press, Princeton and Oxford, 2003.[28] H. R. Thieme. Spectral bound and reproduction number for infinite-dimensionalpopulation structure and time heterogeneity.
SIAM J. Appl. Math. , 70:188–211,2009.[29] L. N. Trefethen.
Spectral methods in MATLAB . Software - Environment - Toolsseries. SIAM, Philadelphia, 2000.[30] J. van Neerven.
The Asymptotic Behaviour of Semigroups of Linear Operators , vol-ume 88 of