This is the accepted version of the journal article: Barril, Carles; Calsina i Ballesta, Àngel; Cuadrado Gavilán, Sílvia; [et al.]. «On the basic reproduction number in continuously structured populations». Mathematical methods in the applied sciences, Vol. 44, Issue 1 (January 2021), p. 799-812. DOI 10.1002/mma.6787 This version is available at https://ddd.uab.cat/record/268839 under the terms of the license On the basic reproduction number in continuously structured populations Carles Barril* , Àngel Calsina* , Sı́lvia Cuadrado* , and Jordi Ripoll** * Departament de Matemàtiques, Universitat Autònoma de Barcelona ** Departament d’Informàtica, Matemàtica Aplicada i Estadı́stica, Universitat de Girona carlesbarril@mat.uab.cat, acalsina@mat.uab.cat, silvia@mat.uab.cat, jripoll@imae.udg.edu July 23, 2020 Abstract In the framework of population dynamics, the basic reproduction number R0 is, by definition, the expected number of offspring that an individual has during its lifetime. In constant and time periodic environments it is calculated as the spectral radius of the so-called next-generation operator ([12, 19]). In continuously structured popula- tions defined in a Banach lattice X with concentrated states at birth one cannot define the next-generation operator in X. In the present paper we present an approach to compute the basic reproduction number of such models as the limit of the basic repro- duction number of a sequence of models for which R0 can be computed as the spectral radius of the next-generation operator. We apply these results to some examples: the (classical) size-dependent model, a size structured cell population model, a size struc- tured model with diffusion in structure space (under some particular assumptions) and a (physiological) age-structured model with diffusion in structure space. Key words: next generation operator; basic reproduction number ; physiologically struc- tured population. 1 Introduction In a deterministic model of population growth, the basic reproduction number R0 is defined as the expected number of offspring that an individual has throughout his life (in an epi- demiological model, as the expected number of new infections a newly infected individual will produce). The main interest of the basic reproduction number is that, as it is intuitively clear, a small population will (begin to) establish in a given environment if R0 > 1 whereas it will become extinct whenever R0 < 1. Similarly, provided some infected individual is present, an epidemic outbreak will arise if R0 > 1 whereas the opposite strict inequality guarantees that the epidemic will not occur. In structured populations, when the birth event can happen in different individual states (of size, phenotype, spatial position, etc.) one talks of “typical” individual but it is not always clear what “typical” individual means and so, which is the expected number of offspring of a ”typical” individual. 1 On the other hand, in deterministic modeling one could be concerned with a definition that involves a probability concept as expectation and would like to have a deterministic definition as well. Both drawbacks are solved by defining R0 as the spectral radius of the so-called next- generation operator, which maps a (distribution of) population to the (distribution of) pop- ulation of their offspring along the whole life span of the former ([11], [27], [12], [19]). Nevertheless, in some models the definition of the birth operator and that of the next- generator operator depend on the choice of the sometimes arbitrary concept of birth event ([10], [4]). Hence more than one basic reproduction number can be meaningfully defined in some models ([10], [4]). When the population dynamics can be described by a linear system of ODEs, the next- generation operator is a matrix for which the spectral radius can be computed and therefore there is no problem with the definition above leaving apart that it may be non-uniquely determined, as commented in the previous paragraph. Fortunately (we should better say inevitably) the sign of R0 − 1 coincides with that of the real part of the so-called Malthus constant (the first eigenvalue of the matrix defining the system of ODEs), independently of the choice of the decomposition of the latter in a mortality/transition operator and a birth operator provided some natural assumptions hold ([11], [10]). The sign relation also holds in infinite dimension [12]. However, in continuously structured populations with concentrated state at birth, the birth rate shows as a boundary condition instead of a birth operator. This makes it difficult (or impossible) to obtain the next-generation operator and thus adapt to these type of models in infinite dimensional spaces (which we will call type II models) the above definition of R0 . In order to overcome this difficulty, in the present paper we consider some examples where R0 is defined for type II models as a limit of basic reproduction numbers of approximate models with distributed states at births for which the next generation operator can be defined as a bounded linear operator, which we call type I models. More explicitly, we will consider diffusion-convection partial differential equations of the form of conservation laws for the density of individuals with respect to some one-dimensional continuous variable which stands for a physiological state as size or age or for an external state as spatial position, with non-flux boundary conditions and with an incoming flow of new individuals via a birth operator whose range is the span of a function in the space of states giving the distribution of the offspring (for which the next generation operator is well defined). Next we will assume a sequence of type I models such that this distribution of offspring concentrates at some point in the closure of the domain of the structuring variable. This sequence tends to a conservation law where the birth term appears as a boundary condition and consequently the next generation operator is not definable as a bounded linear operator in the state space. Many of the classic age and size structured population models are of this second class which we call type II models. For this we define its basic reproduction number as the limit (it exists uniquely under suitable hypotheses) of the basic reproduction numbers of the sequence of type I models tending to it. A different approach to overcome the difficulty to define the next generation operator in the case of concentrated states at birth is to consider the space of measures as state space. This is done for discrete-time population dynamics in the recent reference [28]. The paper is organized as follows: in Section 2 we give, in a general framework, the definition of R0 for type II models as the limit of the sequence of basic reproduction numbers of type I models given by conservation laws for the density of individuals with respect to some one-dimensional continuous variable. Section 3 is devoted to the application of the 2 results in Section 2 to some examples: the (classical) size-dependent model, a size structured cell population model, a size structured model with diffusion in structure space (under some particular assumptions) and a (physiological) age-structured model with diffusion in structure space. For this last example, the computation of the basic reproduction number for the sequence of type I models is done in two different ways, one of them included in the Appendix. At the end of Section 3, a brief discussion on the dependence of R0 on the diffusion coefficient is undertaken with the conclusion that, for some age specific birth functions, there is an optimal value of the diffusion coefficient which maximizes R0 . Therefore, the presence of a moderate diffusion in age tends to increase R0 and hence the survival chances of a population. 2 General framework Some models in population dynamics (which we call type I models, see [4]) can be described by non-linear abstract differential equations for the density of individuals with respect to some structuring variables. Specifically, let X be a Banach lattice and let u(t) ∈ X be the distribution of individuals at time t ≥ 0, then the original non-linear ecological problem can be decomposed as u0 (t) = B(u(t))u(t) − M (u(t))u(t) , u(0) = u0 ∈ X (1) where for any u, B(u) is the linear birth operator and for any u, M (u) is the linear operator corresponding to non-birth terms (mortality and transitions in general). Of course, as men- tioned before, this decomposition depends on what is considered as a birth event ([10, 5]). For the study of possible extinction or settlement of the population we can focus on the so-called extinction steady-state u∗ = 0 and its stability which is analyzed by means of the following formal linearized model u0 (t) = Bu(t) − M u(t) , u(0) = u0 ∈ X . where we set B := B(0) and M := M (0) with some abuse of notation. In order to proceed, we need the following assumptions: B : DB ⊂ X → X is a positive linear operator, and M : DM ⊂ X → X is such that −M is the generator of a strongly continuous semigroup of positive linear operators T (t) := e−M t , whose spectral bound is strictly negative s(−M ) < 0 (so that the population goes to extinctionRwhen there are no ∞ births). Therefore, 0 belongs to the resolvent set of M and there exists 0 T (t) dt = M −1 . Moreover we assume that DM ⊂ DB and that the linear operator BM −1 is bounded. That is guaranteed for instance when B is bounded but this is not a necessary condition (see [5], Sect. 5). Finally we assume that B − M is the infinitesimal generator of a positive semigroup. The operator BM −1 can be interpreted as the next generation operator. For the finite dimensional case, see for instance [10]. In the infinite dimensional case, if B is bounded see [12] and Section 3 in [5]. If B is not bounded, BM −1 is still interpreted as the next- generation operator, see [27] and Section 5 in [5]. For type I models, we have that the next-generation operator BM −1 is defined on the same Banach space X as the population density. The basic reproduction number R0 (which is an alternative to the Malthusian parameter) is then given by the spectral radius of the positive bounded linear operator BM −1 ([12, 5]). Moreover, R0 is a non-negative spectral value which is actually the largest λ ≥ 0 for which the operator B − λM does not have a bounded inverse, see e.g. [4]. 3 Often R0 cannot be computed explicitly and then the following upper bound is useful to find a sufficient condition for extinction kBM −1 φkX   −1 R0 ≤ kBM kX = sup : φ ∈ X \ {0} = kφkX   kBψkX = sup : ψ ∈ DM \ {0} . kM ψkX It is important to notice that, although type I models contain a large family of popu- lation models (in particular, those of discrete structure, for which dim X < ∞), often some structured models cannot be cast into this form. Indeed, beyond the “standard” case considered above where R0 = ρ(BM −1 ), we can consider another family of population models (which we call type II models as in [4]) where the population distribution at time t ≥ 0 is such that u(t) ∈ X and the distribution of newborns per time-unit does not belong to X, and so the model cannot be written in the form (1). For instance, in the case of age-structured populations, all newborns are concentrated at age 0 which translates as a boundary condition for an equation of the form u0 (t) = −M (u(t))u(t). More in general this is very often the case when the individual state space, i.e., the set where the structuring variable lives (its domain) is a continuum, for instance a real interval, and the set of states at birth is discrete. In the forthcoming we will consider models given by conservation laws for densities of individuals living on an interval I, where often the newborns show as a flow crossing the boundary of this domain. Notice that in the case of type II models one cannot, in principle, define a next gener- ation operator. In order to define R0 for type II models, following the approach developed in [4], we will consider a specially-built sequence of type I models in X = L1 (I) u0 (t) = Bk u(t) − M u(t) , u(0) = u0 ∈ X, k ∈ N . (2) The linear operator M : DM ⊂ X → X, as before, is such that −M is the generator of a strongly continuous semigroup of positive linear operators whose spectral bound is strictly negative s(−M ) < 0 and Bk : DBk ⊂ X → X are positive linear operators, all of finite rank m < ∞ such that Bk M −1 is bounded and Bk − M is the generator of a strongly continuous semigroup. To fix ideas let us assume that m = 1. We further assume that, for any k, Bk u = (Lu) ϕk for some positive linear functional L and some positive ϕk ∈ X. Notice that ϕk represents the distribution of offspring. Moreover we assume that the sequence M −1 ϕk converges in X to a certain function ψ∞ ∈ DL = DBk , and that LM −1 ϕk converges to Lψ∞ . Notice that the last hypothesis (and the fact that ψ∞ ∈ DL ) automatically hold if B is bounded. Then we will have that the sequence of basic reproduction numbers corresponding to the approximate models has a limit, which we claim can be named as the basic reproduction number of the limit type II model: R0 := lim R0,k = lim ρ(Bk M −1 ) = lim ρ((LM −1 ·)ϕk ) = lim LM −1 ϕk = Lψ∞ , (3) k→∞ k→∞ k→∞ k→∞ where ρ(·) stands for the spectral radius of the positive bounded linear operators involved and where we have used that range(Bk ) = span{ϕk }. Of course this definition needs a proof of correctness in the sense that it does not depend on the sequence ϕk chosen to approximate the limit model. More precisely, we have to clarify 4 what we do understand by limit model of a sequence of type I models of the form (2) with Bk u = (Lu) ϕk and check that two sequences of functions ϕk and ϕ̂k giving rise to the same limit model are such that M −1 ϕk and M −1 ϕ̂k tend to the same limit. Thus let us consider a sequence of conservation laws in the form Z b Z b Z b d u(s, t)ds = Φ(a, t) − Φ(b, t) − µ(s)u(s, t)dt + Lu(·, t) ϕk (s)ds (4) dt a a a for the rate of change of the population number of individuals whose individual state (size, spatial position, age, etc.) belongs to the (arbitrary) interval (a, b) ⊂ I. Here Φ(x, t) is the flow crossing the point x at time t and we assume that it is given by a linear combination of the density u and its partial derivative us , µ is the state-dependent per capita death rate (so the integral above gives the rate of loss of individuals with individual state within the interval [a, b]). We also assume that the new individuals R are being born according to a fixed distribution of individual states ϕk ∈ X, with I ϕk (s)ds = 1. Hence the last term stands for the population birth rate. This leads, for a smooth density u, to the following sequence of partial differential equations  ∂t u(x, t) = −∂x γ(x)u(x, t) − D(x)∂x u(x, t) − µ(x)u(x, t) + Lu(·, t) ϕk (x) (5) (where D(x) > 0 or D(x) ≡ 0), with non-flux boundary conditions at the (possibly empty) boundary of I. In (5) we can firstly identify an unbounded linear operator, formally written  (M u)(x) = ∂x γ(x)u(x) − D(x)∂x u(x) + µ(x)u(x), with non-flux boundary conditions, for transition between states and mortality and secondly rank one linear operators (Bk u)(x) = Lu ϕk (x) corresponding to reproduction and verify that the model is of the type defined in (2). Let us further assume that the sequence of states distributions at birth ϕk tends to con- centrate, for instance around a point in the closure of I. More precisely, that there exists Rb an individual state x0 ∈ I such that for any [a, b] ⊂ I, limk→∞ a ϕk (s)ds = 1 whenever Rb x0 ∈ [a, b] and limk→∞ a ϕk (s)ds = 0 otherwise. Assuming the second case (i.e. x0 ∈ / [a, b]) and making k go to infinity in (4) we simply obtain the partial differential equation  ∂t u(x, t) = −∂x γ(x)u(x, t) − D(x)∂x u(x, t) − µ(x)u(x, t) (6) for x 6= x0 . On the other hand, if [a, b] contains x0 , making k go to infinity in (4) we first obtain the integral equation Z b Z b d u(s, t)ds = Φ(a, t) − Φ(b, t) − µ(s)u(s, t)dt + Lu(·, t) dt a a and afterwards, making that the interval [a, b] collapses to x0 , we obtain the condition 0 = lim− Φ(a, t) − lim+ Φ(b, t) + Lu(·, t). (7) a→x0 b→x0 This can be seen as a boundary condition when x0 belongs to the boundary of I. For instance if x0 is the left endpoint of I we necessarily have a = x0 and lima→x− Φ(a, t) = Φ(x0 , t) = 0 0 5 for the non-flux boundary condition and (7) reduces to a (in general) Robin type boundary condition lim Φ(b, t) = γ(x0 )u(x0 , t) − D(x0 )∂x u(x0 , t) = Lu(·, t). (8) b→x+ 0 Notice, as particular cases, that the latter supplementing the PDE (6) with D(x) ≡ 0 corresponds to the well known equation of the linear size dependent population dynamics when individuals are all born with the same size (and in particular to the equation of the linear age dependent population dynamics if, furthermore, γ(x) ≡ 1) whereas when D(x) is strictly positive the PDE (6) could be interpreted as a linear convection-diffusion equation (purely diffusive if γ(x) ≡ 0) for a population living in a one-dimensional habitat such that, in order to reproduce, adult individuals travel instantly to the point x0 (a common “nest”) and instantly go back after, to where they were living. On the other hand, if x0 belongs to the interior of I, (7) causes a discontinuity in the flux at the point x0 . Biologically this could correspond to a situation like the latter, but where the “nest” is not at the boundary of the domain. From the point of view of the pde (6) one can consider it on two subdomains of I separated by x0 and with “boundary conditions” at x0 given by   lim+ γ(x)u(x, t) − D(x)∂x u(x, t) − lim− γ(x)u(x, t) − D(x)∂x u(x, t) = Lu(·, t). (9) x→x0 x→x0 As one could expect, we define the limit model of a sequence of type I models of the form (5) (a particular case of the abstract form (2)) as the one given by the partial differential equation (6) plus boundary conditions of the form (8) or (9). Next we show that different sequences ϕk leading to the same limit model (i.e. concentrating at the same point x0 ) share the same limit of the sequence M −1 ϕk and thus that the basic reproduction number of the limit model (given by (3)) is well defined. Indeed, let us denote by G(x, s) the Green’s function of the operator M , i.e., such that one can write Z (M −1 ϕ)(x) = G(x, s)ϕ(s)ds, I (for a comprehensive treatise on Green’s functions one can see [26]). Then, due to the hypotheses on ϕk , Z ψ∞ (x) = lim (M −1 ϕk )(x) = lim G(x, s)ϕk (s)ds = G(x, x0 ), k→∞ k→∞ I independently of the choice of the sequence. Moreover we can now write R0 := lim R0,k = lim LM −1 ϕk = Lψ∞ = LG(·, x0 ). (10) k→∞ k→∞ 3 Examples 3.1 A size dependent model As a first example let us consider the type II model associated to the classical (linear) size dependent population dynamics with potentially unbounded size. Then, (6) reduces to  ∂t u(x, t) = −∂x γ(x)u(x, t) − µ(x)u(x, t) =: (M u)(x, t) x ∈ (x0 , +∞) 6 for u(·, t) ∈ L1 (x0 , ∞), where x0 denotes the individual’s size at birth and with boundary condition Z ∞ γ(x0 )u(x0 , t) = β(s)u(s, t)ds =: (Lu)(t) x0 (a particular form of (8)). We assume that the individual growth rate γ(x) and the mortality rate µ(x) are smooth and bounded below by a positive number and that β(x) is bounded and such that β(x) ≥R0. Moreover, to prevent that the individuals reach infinite size in ∞ 1 finite time we assume x0 γ(x) dx = ∞. We can choose a sequence ϕk (x) concentrating at x0 and consider this model as a limit of type I models of the form (5) with D(x) ≡ 0. Notice that −M and Bk − M (where Bk u = (Lu)ϕk ) are both generators of strongly continuous semigroups. Using the variation of constants formula it is straightforward to compute Rs µ(σ) Rx µ(s) Rx dσ − ds M −1 v (x) 1  = x0 v(s)e x0 γ(σ) ds γ(x) e x0 γ(s) R∞ Rx µ(σ) 1 − dσ = x0 γ(x) e s γ(σ) H(x − s)v(s)ds, Rx µ(σ) 1 − dσ where H is the Heaviside function. Hence, G(x, s) = γ(x) e s γ(σ) H(x − s), and we have 1 − Rxx µ(σ) dσ ψ∞ (x) = G(x, x0 ) = e 0 γ(σ) . γ(x) Using (10) we obtain Z ∞ β(x) − Rxx µ(s) ds R0 = Lψ∞ = e 0 γ(s) dx (11) x0 γ(x) which coincides with the basic reproduction number in size dependent population dynamics obtained in the literature ([18], [19]) as the number of newborns that an individual is ex- pected to produce during his reproductive life. 3.2 A size-structured cell population model Let us start by considering a cell population structured by size x and such that cells divide when they reach a given size (normalized to 1) giving rise to two daughter cells whose size is not necessarily the same. Indeed, let us assume x ∈ (0, 1) and a size distribution of the newborn cells given by a probability density function ϕ(x). Conservation of the mass during cell division forces that ϕ(x) = ϕ(1 − x), i.e. the distribution of newborn cells is symmetric with respect to half the mother cell size. In this way we obtain a type I model given by the PDE  ∂t u(x, t) = −∂x (γ(x)u(x, t)) − µ(x)u(x, t) + 2γ(1)u(1, t)ϕ(x), x ∈ 0, 1 (12) with zero flux at x = 0 for the density of cells u(x, t) and where γ is the individual growth rate, µ the mortality rate and the last term is the (distributed) birth rate. We assume the same hypotheses on γ and µ of the previous section. The PDE is of the form (5) with D(x) ≡ 0 as before. Notice that in this case the birth operator (Bu)(x) = (Lu)ϕ(x) = 2γ(1)u(1)ϕ(x), despite being of rank 1, is only relatively bounded with respect to the linear operator (−M u)(x) = −∂x (γ(x)u(x)) − µ(x)u(x) with non-flux boundary condition at the left hand end x0 = 0 (i.e. B is not bounded as a linear operator in L1 but BM −1 is bounded, which 7 amounts to that LM −1 is a bounded linear form). In order to be in the hypotheses of Section 2 we need B − M to be the infinitesimal gener- ator of a strongly continous positive semigroup. The fact that B is relatively bounded is in general not enough to ensure that B − M is a generator of a strongly continuous semigroup [13]. Nevertheless it suffices that B is a bounded linear operator in the domain of (−M ) endowed with the graph norm (see [13], Cor. III.1.5). For (12) this is satisfied provided that ϕ ∈ D(−M ) = {u ∈ W 1,1 (0, 1) : u(0) = 0}. Indeed, we have for u ∈ DM , and using that 0 belongs to the resolvent set of the operator M , kBukDM = kM Buk = kM (Lu)ϕk =| Lu | kM ϕk =| LM −1 M u | kM ϕk ≤ kM ϕk kL−1 M k kM uk = CkukDM , where the norms without subscript are L1 −norms. Since the mortality/transition operator M is exactly the same as in the previous section, the Green’s function is R x µ(σ) e− s γ(σ) dσ G(x, s) = H(x − s). γ(x) So, the basic reproduction number for this model is given by Z 1 Z 1 R 1 µ(σ) R0 = LM −1 ϕ = 2γ(1)(M −1 ϕ)(1) = 2γ(1) G(1, s)ϕ(s)ds = 2 e− s γ(σ) dσ ϕ(s)ds. 0 0 A bit less realistic but apparently simpler, one can alternatively assume that the division is perfectly symmetric. To deal with this assumption, similarly to what we have done in the previous example, we can choose a sequence of functions ϕk ∈ D(−M ) concentrating 1 R1 at x0 = 1/2 (for instance ϕk (x) = ak x(1 − x)e−k|x− 2 | with ak such that 0 ϕk = 1) and consider the model with symmetric division as a limit of a sequence of type I models as the ones described above (of the form (5)). As limiting equation we get ∂t u(x, t) = −∂x (γ(x)u(x, t)) − µ(x)u(x, t), (13) with “boundary condition” (according to (9))   1 γ( ) lim u(x, t) − lim u(x, t) = Lu(·, t) = 2γ(1)u(1, t), (14) 2 x→( 12 )+ x→( 12 )− which can be considered a birth term giving the influx of new cells (double of the number of cells that disappear because they divide). We have that Z 1 M −1 ϕk = G(·, s)ϕk (s)ds 0 tends to G(·, 21 ) in L1 (0, 1) and uniformly on (0, 1) r ( 12 − ε, 12 + ε) (for any ε < 12 ). Therefore, the basic reproduction number for the limit model can be computed as Z 1 R µ(σ) −1 1 − 11 dσ R0 = lim R0,k = lim LM ϕk = lim 2γ(1) G(1, s)ϕk (s)ds = 2γ(1)G 1, = 2e 2 γ(σ) . 0 2 Finally notice that, for any initial condition, the solution of the initial value problem for (13)-(14) will vanish identically on the interval (0, 12 ) after a finite time which means that we can consider equivalently the PDE (13) on the interval ( 12 , 1) with boundary condition γ( 12 )u( 12 , t) = 2γ(1)u(1, t). 8 3.3 Size-structured particular model with diffusion Let us consider a closed population living in a specific habitat where individuals are classified according to some biometric measure (physiological size) which we assume to evolve in a diffusive way, i.e. random fluctuation in the individual growth. For some examples on structured population models in structure space see for instance [17], [24], [15]. Let X = L1 (x0 , ∞), x0 ≥ 0, and let u(·, t) ∈ X be the density with respect to size at time t ≥ 0. The non-linear problem is described as:      ∂t u(x, t) + ∂x γ(x, S(t))u(x, t) − D∂x u(x, t) + µ(x, S(t))u(x, t) = 0  R∞ γ(x0 , S(t))u(x0 , t) − D∂x u(x0 , t) = x0 β(x, S(t))u(x, t) dx (15)   R∞ S(t) = x0 σ(x)u(x, t) dx  where S(t) is a weighted population size with σ(x) ≥ 0, γ(x, S(t)) > 0 is the individual growth rate, i.e. dx dt = γ(x, S(t)), x(0) = x0 , D > 0 is the diffusion coefficient, β(x, S(t)) ≥ 0 is the fertility rate and µ(x, S(t)) ≥ 0 is the mortality rate. Notice that all newborns are assumed to have size x0 ≥ 0, and growth and vital processes here are density-dependent accounting for limited resources (crowding effects) and size-specific (heterogeneity of the population). Finally, notice that we would get the classical size-structured model if we had D ≡ 0. The linearization around the trivial steady-state is given by: (   ∂t u(x, t) + ∂x γ(x)u(x, t) − D∂x u(x, t) + µ(x)u(x, t) = 0 R∞ (16) γ(x0 )u(x0 , t) − D∂x u(x0 , t) = x0 β(x)u(x, t) dx where we set γ(x) := γ(x, 0), β(x) := β(x, 0) and µ(x) := µ(x, 0) with an abuse of notation. We assume the same hypotheses on γ(x), µ(x) and β(x) as in Section 3.1. Remark. Assuming the nonlinearity constant S in (15), and having then γ(x, S), β(x, S) and µ(x, S) as functions of x uniquely, one could define the so-called net reproduction func- tion for the nonlinear model R(S). This is useful if one is interested in finding nontrivial equilibria of (15) since the total weighted population at equilibrium is given by the solution of R(S) = 1 (see [14]). We can choose ϕk (x) = k 1[x0 ,x0 +1/k] (x) and consider model (16) as a limit of the following sequence of type I models in X = L1 (x0 , ∞): ∂t u(x, t) + ∂x γ(x)u(x, t) − D∂x u(x, t) + µ(x)u(x, t) = x0 β(x)u(x, t)dx · k 1[x0 ,x0 +1/k] (x) (   R∞ γ(x0 )u(x0 , t) − D∂x u(x0 , t) = 0 , k∈N. Here, the birth operators are rank one bounded linear operators of the form Bk : X → X, Z ∞ β(x)φ(x)dx · k 1[x0 ,x0 +1/k] (x)  Bk φ (x) = x0 and the transition/mortality operator is given by the unbounded linear operator M : DM ⊂ X → X, M φ (x) = (γ(x)φ(x) − Dφ0 (x))0 + µ(x)φ(x)  with domain DM = {φ ∈ X : φ0 , (γφ − Dφ0 )0 ∈ X, and γ(x0 )φ(x0 ) − Dφ0 (x0 ) = 0}. 9 The fact that −M is the generator of a positive semigroup goes back to the work by Feller, Hille and Yosida (see [16] and the references therein). Since Bk is bounded, Bk − M is the infinitesimal generator of a positive semigroup. For each k ≥ 1, R0,k is computed as the largest λ ≥ 0 for which the operator Bk − λM does not have a bounded inverse, that is, we have to study the problem Bk ψ − λM ψ = ξ, for ψ ∈ D, ξ ∈ X, which is, β(x)ψ(x)dx · k 1[x0 ,x0 +1/k] (x) − λ (γ(x)ψ(x) − Dψ 0 (x))0 + µ(x)ψ(x) = ξ(x) R∞   x0 γ(x0 )ψ(x0 ) − Dψ 0 (x0 ) = 0 and integrating the first equation along the whole size-span we get Z ∞ Z ∞   β(x) − λµ(x) ψ(x) dx = ξ(x) dx. x0 x0 There is a special case in which the basic reproduction number has an explicit expression. Namely, the “proportional vital rates case”: β(x) = β̃ · µ(x) , β̃ ≥ 0 (17)  R∞ R∞ where the relation above becomes β̃ − λ x0 µ(x)ψ(x) dx = x0 ξ(x) dx and therefore the relation cannot be inverted if λ = β̃ = R0,k . Finally, we get R0 = limk→∞ R0,k = β̃ which in this case is independent of the size-dependent growth and diffusion processes. Notice that (11) gives the same R0 under hypothesis (17): Z ∞ β(x) − Rxx µ(y) dy R0 = e 0 γ(y) dx = β̃ . x0 γ(x) 3.4 Age structured model with diffusion in age In this section we consider a linear age-structured population model for the dynamics of a closed population of individuals (animals, cells, ...) structured by physiological age a ≥ 0 with diffusion D ≥ 0: ( ∂t u(a, t) + ∂a u(a, t) + µ u(a, t) = D ∂aa u(a, t), R∞ (18) u(0, t) − D ∂a u(0, t) = 0 β(a)u(a, t) da. We interpret physiological age as an internal variable associated to the physiological de- velopment of an individual, normally distributed and with expected value equal to the chronological age. The physiological (or biological) age of an individual, as opposed to chronological age or time since birth, varies in a diffusive way. For papers on diffusion age-structured problems see [7],[21]. Notice that (18) is a particular case of (16). The main goal is to compute the basic reproduction number RD 0 using (3) for model (18) and to compare it with the classical age-structured population model corresponding to D = 0. Looking for exponential solutions of system (18) one obtains a characteristic equation from which the following threshold can be derived [9]: Z ∞ √ D 2 a R̃0 = √ β(a) e(1− 1+4Dµ) 2D da . (19) 1 + 1 + 4Dµ 0 10 which tends to the classical expression when the diffusion coefficient tends to zero (by an easy application of the Lebesgue dominated convergence theorem) Z ∞ D lim R̃0 = β(a) e−µa da . D→0 0 This quantity R̃D D 0 might not be equal to the basic reproduction number R0 that we D will compute using the framework set in Section 2 but the sign of R̃0 − 1 must coincide with the one of RD D 0 − 1. In order to obtain R0 we start by computing the inverse of the transition/mortality operator for system (18) which is defined as M φ = φ0 + µφ − Dφ00 . Its inverse can be computed as the solution v ∈ L1 (0, ∞) of the boundary value problem ( v 0 (a) + µ v(a) − Dv 00 (a) = φ(a), (20) v(0) = Dv 0 (0), for a given datum φ ∈ L1 (0, ∞). Let us call 1 p 1 p λ1 = (1 + 1 + 4Dµ) > 0 , λ2 = (1 − 1 + 4Dµ) < 0. 2D 2D The general solution of the second order inhomogeneous linear ODE above can be written as:  Z ∞ Z a  1 v(a) = √ eλ1 a (c1 + e−λ1 s φ(s)ds) + eλ2 a (c2 + e−λ2 s φ(s)ds) . 1 + 4Dµ a 0 By imposing that v is integrable we getR ∞c1 = 0, and using the boundary condition at a = 0 straightforwardly we obtain c2 = − λλ12 0 e−λ1 s φ(s)ds, and thus the solution of (20) can be written as: Z ∞ Z a Z ∞  1 λ1 (a−s) λ2 (a−s) λ2 λ2 a−λ1 s v(a) = √ e φ(s)ds + e φ(s)ds − e φ(s)ds . 1 + 4Dµ a 0 λ1 0 Therefore we can write Z ∞ (M −1 φ)(a) = G(a, s)φ(s)ds, (21) 0 where   1 λ1 (a−s) λ2 (a−s) λ2 λ2 a−λ1 s G(a, s) = √ e H(s − a) + e H(a − s) − e 1 + 4Dµ λ1 Along the lines of [4] and Section 2 we introduce a sequence of “type I” models which approximate (18): β(a)u(a, t) da k 1[0,1/k] (a) ( R∞ ∂t u(a, t) + ∂a u(a, t) + µ u(a, t) = D ∂aa u(a, t) + 0 (22) u(0, t) − D ∂a u(0, t) = 0. Notice that (22) has a biological meaning on its own. That is, it corresponds to assuming that the newborns size (or physiological age) is distributed uniformly on a small interval 11 around the minimal size 0. Similarly to Subsection 3.2, this is an assumption at least as reasonable as assuming that the size of the newborns is concentrated at a given value, for instance 0. Here the birth operator Bk u := 0 β(a)u(a, t) da k 1[0,1/k] has a one-dimensional range R∞ R∞ spanned by the characteristic function of the interval [0, 1/k]. So, Lv := 0 β(a)v(a)da is a bounded linear form given by the L∞ function β. Consequently, the basic reproduction number for (22) can be computed as Z ∞ R0,k = ρ Bk M −1 = LM −1 (k 1[0,1/k] ) = β(a)M −1 (k 1[0,1/k] )(a) da.  0 In this case, a straightforward computation using (21) gives M −1 (k 1[0,1/k]  h ) (a) =   i 1 (1 − ak) F (λ1 ( k1 − a)) − F (λ2 ( k1 − a)) 1[0,1/k] (a) + F ( λk2 ) − λ2 λ1  = √1+4µD λ1 F ( k ) eλ2 a (23) 1−e−x where F (x) := x for x 6= 0, F (0) = 1. We now proceed to compute the basic reproduction number for the limit model (18). So, let us compute the L1 −limit of the sequence M −1 (k 1[0,1/k] ). First notice that F (x) is an analytic function and that its first derivative is negative and its second derivative is positive. Since F is decreasing and λ2 < λ1 , the coefficient of 1[0,1/k] in (23) is non-positive, implying   1 λ2 λ2 λ1 M (k 1[0,1/k] ) (a) ≤ √ −1 F ( ) − F ( ) eλ2 a .  1 + 4µD k λ1 k On the other hand, let us consider the function G(x) = F ( λλ21 x) − λλ12 F (x), whose derivative G0 (x) = λλ12 F 0 ( λλ21 x) − F 0 (x) is positive for x > 0 since λλ21 < 0 and F 00 (x) > 0. Therefore,  F ( λk2 ) − λλ21 F ( λk1 ) = G( λk1 ) decreases with k. So we obtain the following bound, independent of k, 1 λ2 M −1 (k 1[0,1/k] ) (a) ≤ √ F (λ2 ) − F (λ1 ) eλ2 a .   1 + 4µD λ1 λ (1− λ2 )eλ2 a M −1 (k 1[0,1/k] ) (a) converges pointwise to ψ∞ (a) := √1+4µD  1 and hence the convergence 1 is also in L by the Lebesgue dominated convergence theorem. Notice that the limit function does not belong to the domain of the operator −M because it does not satisfy the boundary condition. Also notice that ψ∞ (a) = G(a, 0). Along the lines of Section 2, we compute the basic reproduction number of system (18) as RD 0 = lim R0,k = Lψ∞ = LG(., 0) λ √ R∞ (1− λ2 )eλ2 a 2 R∞ a = 0 β(a) √1+4µD 1 da = 1+√1+4Dµ 0 β(a) e(1− 1+4Dµ) 2D da , which is the basic reproduction number for the diffusion age-structured population model (18) and coincides with the conjecture given at the beginning of the subsection. For an alternative computation of RD 0 using the solution semigroup we refer the reader to the appendix. For a constant fertility rate, i.e. β(a) ≡ β̄, the basic reproduction number is given by β̄ RD0 = , for all D ≥ 0, µ 12 so it is independent of the diffusion process. For age specific fertilities, an interesting issue is the dependence of RD 0 on the diffusion coefficient D. As expectable, RD 0 does not depend on D if the fertility function is a constant. On the other hand, we already know that it tends to the basic reproduction number without diffusion when D tends to 0. It is also very easy to show that it tends to 0 for D tending to ∞ whenever the fertility function β is integrable, albeit it does it slowly (as √1D ). More precisely, √ √ 2 D Z ∞ √ −2µ a 1 Z ∞ lim D RD 0 = lim √ β(a)e 1+ 1+4Dµ da = √ β(a)da, D→∞ D→∞ 1 + 1 + 4Dµ 0 µ 0 as can be seen using the Lebesgue dominated convergence theorem. As for the biological interpretation, an integrable fertility function combined with a large value of the parameter D entails too many unfertile individuals. However this does not mean that RD 0 is always a monotone (decreasing) function of D. For instance the particular case β(a) = β0 a2 e−a (see [29] for a discussion on the choice of 32β0 D3 this birth function) leads to RD0 = √ √ , which tends to (2D + 1 + 4µD − 1)3 (1 + 1 + 4µD) 2β0 (µ+1)3 when D goes to 0 (the value of R0 for D = 0) but, for µ > 1/2, it has a maximum 8β0 2β0 value 27µ (greater than (µ+1) 3 ) at D = 4µ − 2 (see Figure 1). Biologically, this could be explained because for small values of D individuals die before reaching an age with high fertility whereas for large values of D individuals grow old too fast being then fertile for a short period of their life. On the other hand, a very small fertility function in early ages like the one above leads to think in the case of a maturation age below which the individuals do not reproduce. √ For β0 1− 2D 1+4µD instance, if β(a) vanishes for a < 1 and it is a constant β0 for a > 1, then RD 0 = µ e is an increasing function of D ranging from βµ0 e−µ at D = 0 to βµ0 (the basic reproduction number without maturation delay) at infinity. Then, the presence of a moderate diffusion in age tends to increase R0 and hence the survival chances of a population. Figure 1: A graph of RD 0 showing values larger than 1 for moderate diffusion coefficient. Conclusions In this paper we have considered continuously structured population models given by con- servation laws and we have shown useful formulas for the basic reproduction number R0 based on the linear functional L (birth event), the Green’s function G(x, s) (the integral 13 kernel of M −1 , the inverse of the linear operator M giving transitions, i. e., events other than births) and the offspring distributions ϕ and ϕk . Namely, for type I models: Z R0 = ρ(BM −1 ) = (LM −1 ϕ) = L G(·, s)ϕ(s)ds . I For type I models approximating type II models: Z R0,k = ρ(Bk M −1 ) = (LM −1 ϕk ) = L G(·, s)ϕk (s)ds , k∈N. I And finally, for type II models as limit of type I models: Z R0 := lim R0,k = lim L G(·, s)ϕk (s)ds = LG(·, x0 ) , k→∞ k→∞ I where states at birth are concentrated at x0 ∈ I. Notice that we do not need to choose specific ϕk necessarily since we can apply the last formula above for type II models. This limit approach, in contrast to heuristic approaches, allows us to formally show that what we have computed is actually the basic reproduction number of the model. Nevertheless, if one is interested in type I models (because they are biologically relevant by themselves) or, in a more general setting where we have just a definition: R0 := lim R0,k = lim ρ(Bk M −1 ) , k→∞ k→∞ the choice of the ingredients (i.e. how to define operators Bk approximating B) would be important for several reasons: - Some ϕk would describe more accurately the birth distribution. - The rate of convergence of the sequence would be an indication of the goodness of the chosen sequence. - Simplicity in the computation of R0,k . - The constrains of the model as in the example of symmetric division in 3.2 where we need ϕk to be in the domain of the transition/mortality operator. Acknowledgements This work has been partially supported by the projects MTM2017-84214-C2-2-P and RED2018- 102650-T of the Spanish government. We thank two anonymous referees for their useful comments. This work does not have any conflicts of interest. 4 Appendix: An alternative computation of the basic reproduction number for the age structured model with diffusion using the solution semigroup An explicit expression for the solution semigroup of the system   ∂t v(a, t) + ∂a v(a, t) + µ v(a, t) = D ∂aa v(a, t),  v(0, t) − D ∂a v(0, t) = 0, (24)  v(a, 0) = v0 (a)  14 which is the initial value problem for model (18) with zero fertility rate, is possible by means of the Green’s function. First notice that this initial boundary value problem can be rewritten as an initial boundary problem for the heat equation satisfying a slightly different Robin type boundary a 1 condition by means of the change of dependent variable v(a, t) = e 2D −(µ+ 4D )t u(a, t). Indeed, a straightforward computation shows that v is a solution of (24) if and only if u solves   ∂t u(a, t) = D ∂aa u(a, t),  u(0, t) − 2D ∂a u(0, t) = 0, (25) a − 2D  u(a, 0) = e v0 (a)  A Green’s function for problem (25) can be found in the book [8] (page 602), cataloged as X30 : (a−s)2 (a+s)2 t a+s e− + e− e 4D + 2D    4Dt 4Dt a+s+t GX30 (a, s, t) = √ − 1 − erf √ . 2 πDt 2D 2 Dt So, the solution of (25) is given by Z ∞ s u(a, t) = GX30 (a, s, t)e− 2D v0 (s)ds, 0 and consequently, that of (24) by Z ∞ (T (t)v0 )(a) := v(a, t) = G0 (a, s, t)v0 (s)ds 0 where a 1 s G0 (a, s,t) = e 2D −(µ+ 4D )t e− 2D GX30 (a, s, t) 2(a−s)−t (a−s)2 (a+s)2 a    e− 4Dt +e− 4Dt = e−µt e √πDt 4D 2 − eD 2D 1 − erf a+s+t √ 2 Dt With the aim of writing the next generation operator, after interchanging the integration order we can write Z ∞ Z ∞Z ∞ (T (t)v)(a)dt = G0 (a, s, t) dt v(s) ds 0 0 0 Now, proceeding in a similar manner as we did in Section 3.4, we consider system (22) and, setting vk = k 1[0,1/k] , we can write Z ∞ Z 1/k Z ∞ Z ∞ (T (t)vk )(a)dt = k G0 (a, s, t) dt ds → G0 (a, 0, t) dt, 0 0 0 0 where (a−t)2 a ! e−   −µt 4Dt eD a+t G0 (a, 0, t) = e · √ − 1 − erf √ . πDt 2D 2 Dt Arguing as above, Z ∞ Z ∞ Z ∞ RD 0 = lim R0,k = lim β(a) (T (t)vk )(a) dt = β(a) G0 (a, 0, t) dt da, 0 0 0 15 and the only thing left is the computation of the integral of G0 (a, 0, t) with respect to t on (0, ∞). Indeed, a standard argument based on “completing squares” gives (a−τ )2 a √ ∞ e− e 2D (1− 1+4Dµ) Z 4Dτ −µτ e · √ dτ = 2 √ , (26) 0 πDτ 1 + 4Dµ whereas for the second term we have, Z ∞ a    a Z ∞ Z ∞ eD a+τ eD 2 2 e−µτ 1 − erf √ dτ = e−µτ √ e−x dx dτ 0 2D 2 Dτ 2D 0 a+τ √ π 2 Dτ which, after changing integration order equals a Z ∞ Z 2Dx2 −a+2Dx√x2 − Da eD −x2 √ e e−µτ dτ dx D π √a/D √ a 2Dx2 −a−2Dx x2 − D √ √ a ! ∞ ∞ e(1+µD) D Z Z −(1+2µD)x2 +2µDx a x2 − D −(1+2µD)x2 −2µDx a x2 − D = √ √ e dx − √ e dx µD π a/D a/D a ∞ Z √ Da ! e(1+µD) D (s2 + a )(s2 +(4µD+1) a ) Z D D 1 a/D  = √ √a + e− 4s2 − ds, µD π D 0 2 2s2 a/D a/D where we replaced x = 12 (s + > ad in the first integral and x = 1 p s ), s 2 (s + s ), s ∈ pa (0, d ) in the second one. This finally gives a Z ∞ Z ∞ e(1+µD) D (s2 + a )(s2 +(4µD+1) a ) a − (s2 + Da )(s2 +(4µD+1) a )  − D D D √ e 4s 2 ds − e 4s 2 ds 2µD π 0 0 Ds2 √ √  √ √ a a 1+4µD a 1+4µD ∞ = e4µD e 2D erf 2s + a 1+4µD + e− 2D erf 2s − a 1+4µD 2D  2Ds 2Ds 0 √ √  a√1+4µD √  1 s a 1+4µD  − a 1+4µD s a 1+4µD  ∞ + √1+4µD e 2D erf 2 + 2Ds −e 2D erf 2 − 2Ds 0 a √ e 2D (1− 1+4µD) 1  = 1− √ , 2µD 1 + 4µD which can be checked by direct differentiation. Subtracting the last expression from (26) one easily obtains (19). References [1] B. Abdellaoui, T.M. Touaoula. Decay solution for the renewal equation with dif- fusion. Nonlinear Differential Equations and Applications, 17: 271–288, 2010. [2] A.S. Ackleh, J. Z. Farkas. On the net reproduction rate of continuous structured populations with distributed states at birth. Computers and Mathematics with Applica- tions, 66: 1685–1694, 2013. 16 [3] N. Bacaër, S. Guernaoui., The epidemic threshold of vector-borne diseases with seasonality. The case of cutaneous leishmaniasis in Chichaoua, Morocco. Journal of Mathematical Biology, 53 (3): 421–436, 2006. [4] C. Barril, À. Calsina, J. Ripoll. A practical approach to R0 in continuous-time ecological models. Mathematical Methods in the Applied Sciences, 41 (18): 8432–8445, 2017. [5] C. Barril, À. Calsina, J. Ripoll. On the reproduction number of a gut microbiota model. Bulletin of Mathematical Biology, 79: 2727–2746, 2017. [6] B. Basse, B.C. Baguley, E.S. Marshall, W.R. Joseph, B. Van Brunt, G. Wake, D.J.N. Wall. A mathematical model for analysis of the cell cycle in cell lines derived from human tumors. Journal of Mathematical Biology, 47 (4): 295–312, 2003. [7] J.W. Brewer. The age-dependent eigenfunctions of certain Kolmogorov equations of engineering, economics, and biology. Applied Mathematical Modeling, 13:47–57, 1989. [8] K.D. Cole, J.V. Beck, A. Haji-Sheikh, B. Litkouhi. Heat conduction using Green’s functions. Second edition. Series in Computational and Physical Processes in Mechanics and Thermal Sciences. CRC Press, Boca Raton, FL, 2011. [9] J.Chu, A. Ducrot, P. Magal, S Ruan. Hopf bifurcation in a size-structured popu- lation dynamic model with random growth. J. Differential Equations 247 (3): 956–1000, 2009. [10] J.M. Cushing, O. Diekmann. The many guises of R0 (a didactic note). Journal of Theoretical Biology, 404: 295–302, 2016. [11] O. Diekmann, J. A. P. Heesterbeek, M. G. Roberts. The construction of next generation matrices for compartmental epidemic models. Journal of the Royal Society Interface, 7 (47): 873–85, 2010. [12] O. Diekmann, J.A.P. Heesterbeek, J.A.J.Metz. On the definition and the com- putation of the basic reproduction ratio R0 in models for infectious diseases in hetero- geneous populations. Journal of Mathematical Biology, 28 (4): 365–382, 1990. [13] K-J. Engel, R. Nagel, A short course on operator semigroups, Springer, 2006. [14] J. Z. Farkas, T. Hagen, Stability and regularity results for a size-structured popula- tion model. Journal of Mathematical Analysis and Applications, 328 : 119–136, (2007). [15] J. Z. Farkas, P. Hinow, Physiologically structured populations with diffusion and dynamic boundary conditions. Mathematical Biosciences and Engineering, 8 (2): 503– 513, (2011). [16] W. Feller, The parabolic differential equations and the associated semi-groups of transformations. Annals of Mathematics (2), 55: 468–519, 1952. [17] K. P. Hadeler, Structured populations with diffusion in state space. Mathematical Biosciences and Engineering, 7 (1): 37–49, (2010). [18] M. Iannelli, A. Pugliese. An introduction to mathematical population dynamics. Along the trail of Volterra and Lotka. Unitext, 79. La Matematica per il 3+2. Springer, Cham, 2014. 17 [19] H. Inaba Age-structured population dynamics in demography and epidemiology. Springer, Singapore, 2017. [20] H. Inaba, On a new perspective of the basic reproduction number in heterogeneous environments. Journal of Mathematical Biology, 65 (2): 309–348, 2012. [21] B.K. Kakumani, S.K. Tumuluri. On a nonlinear renewal equation with diffusion, Mathematical Methods in the Applied Sciences, 39 (4): 697–708, 2016. [22] Y.H. Lee, L. Sherbakov, J. Taber, J. Shi. Bifurcation diagrams of population models with nonlinear diffusion. Journal of Computational and Applied Mathematics, 194: 357–367, 2006. [23] P. Michel, T.M. Touaoula. Asymptotic behavior for a class of the renewal nonlinear equation with diffusion. Mathematical Methods in the Applied Sciences, 36(3):323–335, 2012. [24] A. Pugliese, F. Milner. A structured population model with diffusion in structure space. Journal of Mathematical Biology, 77 (6-7): 2079–2102, (2018). [25] H.H. Schaefer, Banach lattices and positive operators, Springer, 1974. [26] I. Stakgold, M. Holst. Green’s functions and boundary value problems. Third edi- tion. Pure and Applied Mathematics (Hoboken). John Wiley & Sons, Inc., Hoboken, NJ, 2011. [27] H. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, SIAM Journal on Applied Mathematics, 70 (1): 188– 211, 2009. [28] H. Thieme, Discrete-time population dynamics on the state space of measures, Math- ematical Biosciences and Engineering, 17 (2): 1168–1217, 2020. [29] G.F. Webb. Theory of nonlinear age-dependent population dynamics. Monographs and Textbooks in Pure and Applied Mathematics, 89. Marcel Dekker, Inc., New York, 1985. 18