On the basic reproduction number in continuously structured populations

In the framework of population dynamics, the basic reproduction number R_0 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. In continuously structured populations 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 reproduction number of a sequence of models for which R_0 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 structured model with diffusion in structure space (under some particular assumptions) and a (physiological) age-structured model with diffusion in structure space.


Introduction
In a deterministic model of population growth, the basic reproduction number R 0 is defined as the expected number of offspring that an individual has throughout his life (in an epidemiological 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 R 0 > 1 whereas it will become extinct whenever R 0 < 1. Similarly, provided some infected individual is present, an epidemic outbreak will arise if R 0 > 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.
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 R 0 as the spectral radius of the so-called nextgeneration operator, which maps a (distribution of) population to the (distribution of) population of their offspring along the whole life span of the former ( [11], [26], [12], [18]).
Nevertheless, in some models the definition of the birth operator and that of the nextgenerator 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 nextgeneration 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 R 0 − 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 R 0 . In order to overcome this difficulty, in the present paper we consider some examples where R 0 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 [27].
The paper is organized as follows: in Section 2 we give, in a general framework, the definition of R 0 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 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 R 0 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 R 0 . Therefore, the presence of a moderate diffusion in age tends to increase R 0 and hence the survival chances of a population.

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 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 mentioned 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 where we set B(0) := B and M (0) := M with some abuse of notation. In order to proceed, we need the following assumptions: B : D B ⊂ X → X is a positive linear operator, and M : D M ⊂ 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 extinction when 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 D M ⊂ D B 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 nextgeneration operator, see [26] 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 R 0 (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, R 0 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].
Often R 0 cannot be computed explicitly and then the following upper bound is useful to find a sufficient condition for extinction It is important to notice that, although type I models contain a large family of population 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 R 0 = ρ(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 u (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 generation operator.
In order to define R 0 for type II models, following the approach developed in [4], we will consider a specially-built sequence of type I models in X = L 1 (I) The linear operator M : D M ⊂ 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 B k : D B k ⊂ X → X are positive linear operators, all of finite rank m < ∞ such that B k M −1 is bounded and B k − 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, B k 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 ψ ∞ ∈ D L = D B k , and that LM −1 ϕ k converges to Lψ ∞ . Notice that the last hypothesis (and the fact that ψ ∞ ∈ D L ) authomatically 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: where ρ(·) stands for the spectral radius of the positive bounded linear operators involved and where we have used that range(B k ) = 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 what we do understand by limit model of a sequence of type I models of the form (2) with B k 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 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 u s , µ 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 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 , with non-flux boundary conditions at the (possibly empty) boundary of I. In (5) we can firstly identify an unbounded linear operator, formally written with non-flux boundary conditions, for transition between states and mortality and secondly rank one linear operators (B k 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 concentrate, for instance around a point in the closure of I. More precisely, that there exists an individual state x 0 ∈ I such that for any [a, b] ⊂ I, lim k→∞ b a ϕ k (s)ds = 1 whenever x 0 ∈ [a, b] and lim k→∞ b a ϕ k (s)ds = 0 otherwise. Assuming the second case (i.e. x 0 / ∈ [a, b]) and making k go to infinity in (4) we simply obtain the partial differential equation On the other hand, if [a, b] contains x 0 , making k go to infinity in (4) we first obtain the integral equation and afterwards, making that the interval [a, b] collapses to x 0 , we obtain the condition This can be seen as a boundary condition when x 0 belongs to the boundary of I. For instance if x 0 is the left endpoint of I we necessarily have a = x 0 and lim a→x − 0 Φ(a, t) = Φ(x 0 , t) = 0 for the non-flux boundary condition and (7) reduces to a (in general) Robin type boundary condition lim 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 x 0 (a common "nest") and instantly go back after, to where they were living.
On the other hand, if x 0 belongs to the interior of I, (7) causes a discontinuity in the flux at the point x 0 . 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 x 0 and with "boundary conditions" at x 0 given by 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 x 0 ) share the same limit of the sequence M −1 ϕ k and thus that the basic reproduction number of the limit model (given by (3) (for a comprehensive treatise on Green's functions one can see [25]).
Then, due to the hypotheses on ϕ k , independently of the choice of the sequence. Moreover we can now write 3 Examples

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 for u(·, t) ∈ L 1 (x 0 , ∞), where x 0 denotes the individual's size at birth and with boundary condition (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) ≥ 0. Moreover, to prevent that the individuals reach infinite size in finite time we assume where H is the Heaviside function. Hence, G(x, s) Using (10) we obtain which coincides with the basic reproduction number in size dependent population dynamics obtained in the literature ( [17], [18]) as the number of newborns that an individual is expected to produce during his reproductive life.

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 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, as before D(x) ≡ 0. 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 x 0 = 0 (i.e. B is not bounded as a linear operator in L 1 but BM −1 is bounded, which 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 generator 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 ∈ D M , and using that 0 belongs to the resolvent set of the operator M , where the norms without subscript are L 1 −norms.
Since the mortality/transition operator M is exactly the same as in the previous section, the Green's function is So, the basic reproduction number for this model is given by 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 at x 0 = 1/2 (for instance ϕ k (x) = a k x(1 − x)e −k|x− 1 2 | with a k such that 1 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
with "boundary condition" (according to (9)) which can be considered a birth term giving the influx of new cells (double of the number of cells that disappear because they divide).

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 [16], [23], [14].
Let X = L 1 (x 0 , ∞), x 0 ≥ 0, and let u(·, t) ∈ X be the density with respect to size at time t ≥ 0. The non-linear problem is described as: ≥ 0 is the fertility rate and µ(x, S(t)) ≥ 0 is the mortality rate. Notice that all newborns are assumed to have size x 0 ≥ 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: 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 . We can choose ϕ k (x) = k1 [x0,x0+1/k] (x) and consider model (15) as a limit of the following sequence of type I models in X = L 1 (x 0 , ∞): Here, the birth operators are rank one bounded linear operators of the form B k : X → X, and the transition/mortality operator is given by the unbounded linear operator M : The fact that −M is the generator of a positive semigroup goes back to the work by Feller, Hille and Yosida (see [15] and the references therein). Since B k is bounded, B k − M is the infinitesimal generator of a positive semigroup.
For each k ≥ 1, R 0,k is computed as the largest λ ≥ 0 for which the operator B k − λM does not have a bounded inverse, that is, we have to study the problem B k ψ − λM ψ = ξ, for ψ ∈ D, ξ ∈ X, which is, and integrating the first equation along the whole size-span we get There is a special case in which the basic reproduction number has an explicit expression. Namely, the "proportional vital rates case": where the relation above becomes β − λ ∞ x0 µ(x)ψ(x) dx = ∞ x0 ξ(x) dx and therefore the relation cannot be inverted if λ =β = R 0,k . Finally, we get R 0 = lim k→∞ R 0,k =β which in this case is independent of the size-dependent growth and diffusion processes.

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: We interpret physiological age as an internal variable associated to the physiological development 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], [20]. Notice that (17) is a particular case of (15). The main goal is to compute the basic reproduction number R D 0 using (3) for model (17) and to compare it with the classical age-structured population model corresponding to D = 0.
Looking for exponential solutions of system (17) one obtains a characteristic equation from which the following threshold can be derived [9]: which tends to the classical expression when the diffusion coefficient tends to zero (by an easy application of the Lebesgue dominated convergence theorem) This quantityR D 0 might not be equal to the basic reproduction number R D 0 that we will compute using the framework set in Section 2 but the sign ofR D 0 − 1 must coincide with the one of R D 0 − 1. In order to obtain R D 0 we start by computing the inverse of the transition/mortality operator for system (17) which is defined as M φ = φ + µφ − Dφ . Its inverse can be computed as the solution v ∈ L 1 (0, ∞) of the boundary value problem for a given datum φ ∈ L 1 (0, ∞). Let us call The general solution of the second order inhomogeneous linear ode above can be written as: By imposing that v is integrable we get c 1 = 0, and using the boundary condition at a = 0 straightforwardly we obtain c 2 = − λ2 λ1 ∞ 0 e −λ1s φ(s)ds, and thus the solution of (19) can be written Therefore we can write where G(a, s) Along the lines of [4] and Section 2 we introduce a sequence of "type I" models which approximate (17): Notice that (21) 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 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 B k u := In this case, a straightforward computation using (20) gives We now proceed to compute the basic reproduction number for the limit model (17). So, let us compute the L 1 −limit of the sequence M −1 (k1 [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 (22) is non-positive, implying On the other hand, let us consider the function So we obtain the following bound, independent of k, )e λ 2 a √ 1+4µD and hence the convergence is also in L 1 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 (17) as which is the basic reproduction number for the diffusion age-structured population model (17) and coincides with the conjecture given at the beginning of the subsection. For an alternative computation of R D 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 so it is independent of the diffusion process. For age specific fertilities, an interesting issue is the dependence of R D 0 on the diffusion coefficient D. As expectable, R D 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 (an expectable result too), albeit it does it slowly (as 1 √ D ). However this does not mean that R D 0 is always a monotonous (decreasing) function of D. For instance the particular case β(a) = β 0 a 2 e −a (see [28] for a discussion on the choice of this birth function) leads to , which tends to 2β0 (µ+1) 3 when D goes to 0 (the value of R 0 for D = 0) but, for µ > 1/2, it has a maximum value 8β0 27µ (greater than 2β0 (µ+1) 3 ) at D = 4µ − 2. 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 instance, if β(a) vanishes for a < 1 and it is a constant β 0 for a > 1, then 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 R 0 and hence the survival chances of a population.

Acknowledgements
This work has been partially supported by the project MTM2017-84214-C2-2-P of the Spanish government.
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 which is the initial value problem for model (17) 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 condition by means of the change of dependent variable v(a, t) = e a 2D −(µ+ 1 4D )t u(a, t). Indeed, a straightforward computation shows that v is a solution of (23) if and only if u solves      ∂ t u(a, t) = D ∂ aa u(a, t), u(0, t) − 2D ∂ a u(0, t) = 0, u(a, 0) = e − a 2D v 0 (a) A Green's function for problem (24) can be found in the book [8] (page 602), cataloged as X30 : So, the solution of (24) is given by u(a, t) = With the aim of writing the next generation operator, after interchanging the integration order we can write Now, proceeding in a similar manner as we did in Section 3.4, we consider system (21) and, setting v k = k1 [0,1/k] , we can write Arguing as above, G 0 (a, 0, t) dt da, and the only thing left is the computation of the integral of G 0 (a, 0, t) with respect to t on (0, ∞).
Indeed, a standard argument based on "completing squares" gives whereas for the second term we have, which can be checked by direct differentiation. Subtracting the last expression from (25) one easily obtains (18).