# On competition through growth reduction

Carles Barril, Àngel Calsina, Odo Diekmann, József Z. Farkas

###### Abstract

We consider a population organised hierarchically with respect to size in such
a way that the growth rate of each individual depends only on the presence of
larger individuals. As a concrete example one might think of a forest, in
which the incidence of light on a tree (and hence how fast it grows) is
affected by shading of taller trees. The model is formulated as a delay
equation, more specifically a scalar renewal equation, for the population
birth rate. After discussing the well-posedness of the model, we analyse how
many stationary birth rates the equation can have in terms of the functional
parameters of the model. In particular we show that, under reasonable and
rather general assumptions, only one stationary birth rate can exist besides
the trivial one (associated to the state in which there are no individuals and
the population birth rate is zero). We give conditions for this non-trivial
stationary birth rate to exist and we analyse its stability using the
principle of linearised stability for delay equations. Finally we relate the
results to an alternative formulation of the model taking the form of a
quasilinear partial differential equation for the population size-density.

We dedicate this paper to Professor Glenn F. Webb, a friend, mentor and
distinguished scientist. Over the past 50 years Glenn has made tremendous
contributions to a wide variety of research domains, ranging from semigroup
theory to cancer modelling. Structured population dynamics has been a central
fixture to Glenn’s research interest for decades. Indeed, Glenn is recognised
as one of the worldwide leading experts of age- and size-structured population
dynamics, an exciting field of research, which has enjoyed tremendous growth
in recent decades. We are happy to have the opportunity to make a small
contribution to this field and this special issue honouring Glenn and
celebrating his achievements.

##  1 Introduction

In terms of numbers, the dynamics of a population is generated by mortality
and reproduction. In structured population models [23, 19, 17], individuals
are characterized by variables such as age, size or other (physiological)
characteristics. In that case, development/maturation needs to be modelled too
(a trivial task in the case of age, but certainly not in general!).

As explained in detail in [7], density dependence is most easily incorporated
in a two step procedure: i) first introduce the environmental condition via
the requirement that individuals are independent from one another when this
condition is prescribed as a function of time; ii) next model feedback by
specifying how, at any particular time, the environmental condition is
influenced by the population size and composition. In the inspiring book [20]
detailed ecological motivation is presented for including in this feedback
loop the impact of density dependence on development and maturation.

Here our aim is to investigate the consequences of density dependence directly
that only affects development (fertility is affected indirectly, since it
depends on the developmental stage of the individual). We do so for a one-
dimensional i-state (i.e., the variable capturing the relevant differences
among individuals ‘lives’ on the real line), so for an i-state space that
comes equipped with an order relation. In fact we shall assume that the
presence of ‘larger’ individuals has a negative impact on the growth rate of
‘smaller’ individuals (as a motivating example one might think of trees and
shading, with the i-state interpreted as ‘height’; but please note that we
ignore spatial structure and that, consequently, the model is but a
caricature).

The organisation of the paper is as follows. In Section 2 we present the
biological assumptions of the model and we deduce a scalar nonlinear renewal
equation for the population birth rate (the so called delay formulation). In
Section 3 a dynamical systems framework for the renewal equation is outlined.
In Section 4 we give conditions guaranteeing the existence of a non-zero
stationary birth rate. In Section 5 we apply the principle of linearised
stability for delay equations [11] to prove that, for a certain two-parameter
family of fertility functions, such a stationary birth rate (whenever it
exists) is locally asymptotically stable. We also show that, under natural
hypotheses on the ingredients, the zero stationary birth rate is a global
attractor when it is the only stationary birth rate.

In Appendix A a technical result needed in Section 5 is shown. In Appendix B a
more classical formulation of the model, taking the form of a first order PDE
involving non-local functionals, is presented. There we show that the
conditions guaranteeing the existence of stationary population densities (with
respect to height) coincide with the conditions guaranteeing non-trivial
stationary birth rates in the delay formulation. This makes sense since both
formulations model the same phenomena (although they are independently derived
from biological assumptions). Such a phenomenological relation between the two
formulations suggests that the stability results for the delay formulation can
be translated to the PDE formulation (as indeed is done in [3]). Although this
issue is not addressed rigorously in the present paper, some comments are
included in the concluding remarks section.

##  2 The model and its delay formulation

Individuals are fully characterized by a variable $x$, taking values in
$\mathbb{R}_{+}$. In general, $x$ is called i-state but here, for clarity, we
call it ‘height’, the point being that we motivate our assumptions about
interaction in terms of competition for light (this phenomenon is also
addressed mathematically in [16, 24, 18], among others). Indeed, we assume
that the growth rate $g$ of an individual of height $x$ does not depend on $x$
directly, but only indirectly, as it depends on the amount of light the
individual receives per unit of time. We assume that the latter, in turn, is
fully determined by the number $E(x,t)$ of individuals that are taller than
$x$ (we call $E$ an interaction variable, since it mediates how the
environmental condition, here light intensity, is influenced by the extant
population). We assume that the per capita death rate $\mu$ and the per capita
reproduction rate $\beta$ only depend on the height $x$. In fact we assume
that $\mu$ is constant, i.e., independent of $x$, while $\beta$ is a non-
decreasing function of $x$.

We assume that all individuals are born with the minimal height $x_{m}$ and
that $g$ is positive (we do not impose an upper bound on height). We assume
that a density function $u=u(x,t)$ exists such that the integral of $u$ with
respect to the first variable over an interval gives the number of individuals
with size within this interval at time $t$. This allows us to write

$E(x,t)=\int_{x}^{\infty}u(s,t)ds,$ (2.1)

so that the height of an individual evolves according to

$X^{\prime}(t)=g(E(X(t),t)).$ (2.2)

Let $B(t)$ denote the population birth rate at time $t$. Then $B$ equals the
influx at $x_{m}$, which originates from reproduction by the extant
population:

$B(t)=\int_{0}^{\infty}\beta(y)u(y,t)dy$ (2.3)

Let $n(t,\cdot)$ denote the age density. We do not need to write a PDE and
solve it in order to conclude that

$n(t,a)=B(t-a)e^{-\mu a}$ (2.4)

This allows us to rewrite (2.3) as

$B(t)=\int_{0}^{\infty}\beta(S(a,t))B(t-a)e^{-\mu a}da$ (2.5)

with $x=S(a,t)$ specifying the height of an individual having age $a$ at time
$t$ (and hence being born at time $t-a$).

We refer to section III.4 of [19], entitled “Integration along
characteristics, transformation of variables, and the following of cohorts
through time”, for general considerations about switching between size- and
age-densities. Here the situation is relatively simple, since the individuals
taller than you are exactly those that are older than you, i.e., were born
earlier than you. Or, in a formula

$E(x,t)=\int_{\tau}^{\infty}B(t-\alpha)e^{-\mu\alpha}d\alpha$ (2.6)

when $x=S(\tau,t)$.

Next note that an individual that was born at time $t-a$ has age $\tau$ at
time $t-a+\tau$. The height $y=y(\tau):=S(\tau,t-a+\tau)$ of such an
individual evolves according to

$\displaystyle\frac{dy}{d\tau}(\tau)$ $\displaystyle=g(E(y(\tau),t-a+\tau))$
(2.7) $\displaystyle=g(E(S(\tau,t-a+\tau),t-a+\tau))$
$\displaystyle=g\left(\int_{\tau}^{\infty}B(t-a+\tau-\alpha)e^{-\mu\alpha}d\alpha\right)$

Noting that $y(0)=x_{m}$ we obtain by integration that

$\displaystyle S(a,t)$ $\displaystyle=y(a)$ (2.8)
$\displaystyle=x_{m}+\int_{0}^{a}g\left(\int_{\tau}^{\infty}B(t-a+\tau-\alpha)e^{-\mu\alpha}d\alpha\right)d\tau$

Inserting (2.8) into (2.5) we obtain

$B(t)=\int_{0}^{\infty}\beta\bigg{(}x_{m}+\int_{0}^{a}g\left(\int_{\tau}^{\infty}e^{-\mu\alpha}B_{t}(\tau-a-\alpha)\mbox{d}\alpha\right)\,\mbox{d}\tau\bigg{)}\,\,e^{-\mu
a}B_{t}(-a)\,\mbox{d}a,$ (2.9)

where

$B_{t}(\theta):=B(t+\theta).$ (2.10)

Notice that (2.9) can also be written as

$B(t)=\int_{0}^{\infty}\beta\left(x_{m}+\int_{0}^{a}g\left(e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}B_{t}(-s)\mbox{d}s\right)\,\mbox{d}\tau\right)\,\,e^{-\mu
a}B_{t}(-a)\,\mbox{d}a.$ (2.11)

######  Remark 2.1.

Notice that (2.9) can be understood directly from the description of the
physical situation: it states that the birth rate at time $t$ is given by the
addition for the age $a$ of the mothers, with density equal to the birth rate
at time $t-a$ times their survival probability, their size specific fertility
with size given by the birth size plus the integral with respect to $\tau$ of
the individual growth rate at age $\tau$, which depends on how many
individuals are larger (integral of the density of individuals of age
$\alpha>\tau$, i.e. $s=\alpha+a-\tau>a$).

##  3 The dynamical systems framework

Equation $(\ref{scalar2})$ provides the delay formulation of the model (see
Appendix B for the alternative PDE formulation). Here the state variable is
the population birth rate history $B_{t}:=B(t+\cdot)$, instead of the
population density $u(\cdot,t)$ with respect to height (characterized in
Appendix B). Specifically, one can consider the state space (of the weighted
birth rate histories)

$\mathcal{X}=L_{\rho}^{1}(-\infty,0):=\left\\{\phi\in
L_{loc}^{1}(-\infty,0):\left|\left|\phi\right|\right|_{\mathcal{X}}=\int_{-\infty}^{0}e^{\rho
s}|\phi(s)|ds<\infty\right\\},$

for some $\rho>0$ (so $\mathcal{X}$ contains, in particular, constant
functions, and therefore the possible steady states) and the delay equation
$B(t)=\mathcal{F}(B_{t})$ with $\mathcal{F}:\mathcal{X}\rightarrow\mathbb{R}$
defined by

$\mathcal{F}(\phi)=\int_{0}^{\infty}\beta\left(x_{m}+\int_{0}^{a}g\left(e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\phi(-s)\,\mbox{d}s\right)\,\mbox{d}\tau\right)\,\,e^{-\mu
a}\phi(-a)\,\mbox{d}a.$ (3.1)

We denote by $\mathcal{X}^{+}$ the standard positive cone in $\mathcal{X}$.

As discussed in [10, 11], the delay equation $B(t)=\mathcal{F}(B_{t})$,
together with an initial history $B_{0}=\phi\in\mathcal{X}$, can be
interpreted as an abstract Cauchy problem with a semilinear structure:

$\left\\{\begin{aligned}
&\frac{d}{dt}\varphi(t)=A\varphi(t)+\mathcal{F}(\varphi(t))\delta_{0}\\\
&\varphi(0)=\phi\in\mathcal{X}\end{aligned}\right.,$ (3.2)

where $A$ is the generator of the linear semigroup defined as
$T_{A}(t)\phi:=\phi(t+\cdot)\mathds{1}_{-}(t+\cdot)$. Notice that the mapping
$t\mapsto T_{A}(t)\phi$ tells us how the population birth rate history would
evolve without considering birth (and growth and mortality), since all these
processes are summarised in the $\mathcal{F}$ function. The previous setting
makes it possible to analyse the well posedness of the problem and some
dynamical properties by means of a generalised variation of constants
equation. The standard variation of constants equation cannot be applied in a
straightforward manner (as in [21]) since the semilinear term of the problem
(namely $\phi\mapsto\mathcal{F}(\phi)\delta_{0}$) does not take values in
$\mathcal{X}$, but in the space of measures.

Here the theory included in the references mentioned above ([10] and [11])
applies provided that $\mathcal{F}$ is continuously differentiable, which is
stated in the following theorem, and proved in Appendix A. We assume that $g$
is smoothly extended to the whole of $\mathbb{R}$, implying that the right
hand side of (3.1) is defined on the whole Banach space $\mathcal{X}$ (so even
for non positive $\phi$). Of course negative birth rates do not have
biological meaning, but they allow us to work on the whole space (recall that
the positive cone of $L^{1}$ has empty interior).

######  Theorem 3.1.

Assume that $g:\mathbb{R}\rightarrow\mathbb{R}$ and
$\beta:\mathbb{R}^{+}\rightarrow\mathbb{R}$ have a bounded and globally
Lipschitzian first derivative. Also assume that $g$ is bounded and positive
and bounded away from 0 and that $\beta$ is non-negative. Then the map
$\mathcal{F}:\mathcal{X}\rightarrow\mathbb{R}$ defined in (3.1) is
continuously differentiable with bounded derivative provided that
$\rho<\mu/5.$

######  Theorem 3.2.

Existence and uniqueness Under the hypotheses of the previous theorem, for any
$\phi\in\mathcal{X},$ there exists a unique $B\in L_{loc}^{1}(\mathbb{R})$
such that $B(t)=\phi(t)$ for $t<0$ and $B(t)$ satisfies (2.11) for $t\geq 0$.
Moreover, $B$ belongs to the positive cone of $L_{loc}^{1}(\mathbb{R})$
whenever $\varphi\in\mathcal{X}^{+}$.

###### Proof.

It is an immediate consequence of Theorem 3.1 (notice that a bounded
derivative implies global Lipschitz continuity), Theorem 3.15 in [11] (which
implies the equivalence between (2.11) and (3.2)) and Theorem 2.2 in [11]
(which implies the existence and uniqueness of mild solutions of (3.2) and the
generation of a nonlinear semigroup $\Sigma(t;\phi)$ satisfying
$\Sigma(t;\phi)=B_{t}$). The facts that the linear semigroups in Theorem 2.2
of [11] are positive and $\mathcal{F}$ maps the positive cone of $\mathcal{X}$
to $\mathbb{R}^{+}$ imply that $B$ belongs to the positive cone whenever
$\phi\in\mathcal{X}^{+}$.

Let $B\in\mathbb{R}$ be a stationary population birth rate, i.e. $B$ satisfies
$B=\mathcal{F}(\bar{B})$ where $\bar{B}\in\mathcal{X}$ is defined by
$\bar{B}(\theta)=B$ for (almost) all $\theta\in(-\infty,0)$. The following
theorem determines the local stability of $\bar{B}$ in terms of properties of
$D\mathcal{F}(\bar{B})$. Since $D\mathcal{F}(\bar{B})$ is a bounded linear
operator from $\mathcal{X}$ to $\mathbb{R}$, the Riesz Representation Theorem
implies that $D\mathcal{F}(\bar{B})$ can be written as

$D\mathcal{F}(\bar{B})\phi=\int_{0}^{\infty}k(s)\phi(-s)ds=:\left<k,\phi\right>$

with $k$ an element of the dual space of $\mathcal{X}$, represented by

$\mathcal{X}^{\prime}=L_{\rho}^{\infty}(0,\infty):=\left\\{f\in
L^{\infty}(0,\infty):\left|\left|f\right|\right|_{\mathcal{X^{\prime}}}=\sup_{s\in(0,\infty)}e^{\rho
s}|f(s)|<\infty\right\\}.$

######  Theorem 3.3.

(Theorem 3.15 in [11]) Under the hypotheses of Theorem 3.1, let
$\bar{B}\in\mathcal{X}$ be a stationary state of (2.11) and let
$k\in\mathcal{X}^{\prime}$ represent $DF(\bar{B})$. Consider the
characteristic equation

$0=1-\hat{k}(\lambda),$ (3.3)

where $\hat{k}$ is the Laplace transform of $k$ (i.e.
$\hat{k}(\lambda)=\int_{0}^{\infty}e^{-\lambda s}k(s)ds$).

  * (a)

If all the roots of the characteristic equation (3.3) have negative real part,
then the stationary state $\bar{B}$ is locally asymptotically stable.

  * (b)

If there exists at least one root with positive real part, then the steady
state $\bar{B}$ is unstable.

##  4 Existence and characterization of steady states

A stationary solution of the problem can be found by simply assuming that $B$
in (2.11) is independent of $t$. Of course there is a trivial stationary
solution $B=0$ that corresponds to the absence of individuals. When dealing
with non-trivial stationary solutions of (2.11), we make the following abuse
of notation to ease readability: we use $\bar{B}$ to denote a constant
function in $\mathcal{X}$ and $\bar{B}\in\mathbb{R}$ as the image it takes (so
that we let the context tell whether $\bar{B}$ refers to the constant function
or the value it takes). With this in mind, and taking into account (2.11), it
follows that a non-trivial equilibrium $\bar{B}\in\mathcal{X}$ is a constant
function whose image is a non-zero solution of

$1=\int_{0}^{\infty}\beta\left(x_{m}+\int_{0}^{a}g\left(B\frac{e^{-\mu\tau}}{\mu}\right)\,\mbox{d}\tau\right)\,\,e^{-\mu
a}\,da=:R(B).$ (4.1)

Under natural hypotheses concerning $\beta$ and $g$, which essentially amount
to assuming that larger sizes correspond to larger fertilities, that more
competition (more individuals higher in the hierarchy than the one we are
observing) means slower growth, and that the first generation progeny of an
individual is finite (more precisely, that $R(0)<\infty$), we readily obtain
the following theorem.

######  Theorem 4.1.

Under the hypotheses of Theorem 3.1 and the assumption that $\beta$ is a
strictly increasing function on $[x_{m},\infty)$, and that $R(0)<\infty$ and
$g$ is a strictly decreasing function on $[0,\infty)$ vanishing at infinity,
there exists a non-trivial equilibrium of (2.11) if and only if

$R_{0}:=R(0)=\int_{0}^{\infty}\beta\big{(}x_{m}+g(0)a\big{)}\,e^{-\mu
a}\,da>1\qquad\text{and}\qquad\frac{\beta(x_{m})}{\mu}<1,$

and there is at most one such non-trivial equilibrium.

###### Proof.

Under the hypotheses, a (double) application of the Lebesgue dominated
convergence theorem yields that $R$ is a well defined continuous and strictly
decreasing function on $[0,\infty)$, tending to $\beta(x_{m})/\mu$ at
infinity.

######  Remark 4.2.

As usual, $R_{0}$ can be interpreted as the so-called basic reproduction
number, i.e., the expected total number of offspring of an individual
experiencing the virgin environment, i.e., when there are no individuals
older/larger than it.

###  4.1 Age and size equilibrium profiles

The age density of a steady state is given by $\bar{n}(a)=\bar{B}e^{-\mu a}$
(see (2.4)).
Let us set

$\bar{S}(a)=x_{m}+\int_{0}^{a}g\big{(}\int_{\tau}^{\infty}\bar{B}e^{-\mu\alpha}\mbox{d}\alpha\big{)}\mbox{d}\tau=x_{m}+\int_{0}^{a}g\big{(}\bar{B}\frac{e^{-\mu\tau}}{\mu}\big{)}\mbox{d}\tau$
(4.2)

(the size of an individual of age $a$ at the nontrivial equilibrium, see
(2.8)).
The density $\bar{u}(x)$ with respect to size of the same population
distribution can then be computed by taking into account the equality

$\int_{\alpha_{1}}^{\alpha_{2}}\bar{n}(a)da=\int_{\bar{S}(\alpha_{1})}^{\bar{S}(\alpha_{2})}\frac{\bar{n}\big{(}\bar{S}^{-1}(x)\big{)}}{\bar{S}^{\prime}\big{(}\bar{S}^{-1}(x)\big{)}}dx=\int_{\bar{S}(\alpha_{1})}^{\bar{S}(\alpha_{2})}\bar{u}(x)dx,$

which follows from the change of variable $x=\bar{S}(a)$ and the
interpretation of $\bar{n}$ and $\bar{u}$. Thus, we find

$\bar{u}(x)=\frac{\bar{n}\big{(}\bar{S}^{-1}(x)\big{)}}{\bar{S}^{\prime}\big{(}\bar{S}^{-1}(x)\big{)}}=\frac{\bar{B}e^{-\mu\bar{S}^{-1}(x)}}{g\big{(}\bar{B}\frac{e^{-\mu\bar{S}^{-1}(x)}}{\mu}\big{)}},$
(4.3)

which is an alternative expression to (B.2).

##  5 Stability of steady states

The linearisation of (2.9) around the origin is simply (see A.2 in Appendix
A),

$y(t)=D\mathcal{F}(0)y_{t}=\int_{0}^{\infty}\beta\big{(}x_{m}+g(0)a\big{)}\,e^{-\mu
a}y(t-a)\mbox{d}a=:\int_{0}^{\infty}k(a)y_{t}(-a)\mbox{d}a$ (5.1)

(as indeed one can understand by using only the interpretation: it describes
the linear population model corresponding to the virgin environment $E=0$).

######  Theorem 5.1.

Under the hypotheses of Theorem 3.1, the trivial equilibrium is (locally)
exponentially stable if $R_{0}<1$ and unstable if $R_{0}>1.$

###### Proof.

Clearly the kernel $k\in L_{\rho}^{\infty}(0,\infty)$ corresponds to the Riesz
representation of $D\mathcal{F}(0)$. Then, according to Theorem 3.3, the
stability of the steady state is determined by the sign of the real part of
the zeroes of the characteristic equation $\hat{k}(\lambda)=1,$ where
$\hat{k}$ stands for the Laplace transform of $k.$
$\hat{k}(\lambda)$ is defined (at least) for $Re(\lambda)>-\rho.$ Moreover,
since the kernel $k$ is positive, $\hat{k}$ is for real $\lambda$ a decreasing
function with limit 0 at infinity. Hence there is at most one real solution
$\hat{\lambda}$ of the characteristic equation, which indeed exists and is
positive if $\hat{k}(0)=R_{0}>1.$ So then the trivial equilibrium is unstable.
When $\hat{k}(0)=R_{0}<1,$ if there is a real root, it is negative. Moreover,
if a non-real $\lambda$ is a root of the characteristic equation, then
$1=\hat{k}(\lambda)=Re(\hat{k}(\lambda))<\hat{k}(Re\lambda)$, which implies,
by the fact that $\hat{k}$ tends to $0$, that there is a real root
$\hat{\lambda}$ larger than $Re(\lambda)$. As such a real root is necessarily
negative, the trivial equilibrium is locally exponentially stable.

######  Theorem 5.2.

If $R_{0}<1$ and the hypotheses of Theorem 4.1 hold, all solutions of (5.1)
tend exponentially to $0$ as $t\rightarrow\infty$.

###### Proof.

For a given solution let us write (cf. (2.8))

$S(a,t)=x_{m}+\int_{0}^{a}g\left(\int_{\tau}^{\infty}e^{-\mu\alpha}B_{t}(\tau-a-\alpha)\mbox{d}\alpha\right)\,\mbox{d}\tau,$

the size at time $t$ of an individual of age $a.$ From (2.9) we can write

$\displaystyle B(t)$ $\displaystyle=\int_{0}^{\infty}\beta(S(a,t))e^{-\mu
a}B(t-a)\,\mbox{d}a$ (5.2)
$\displaystyle=\int_{-\infty}^{0}\beta(S(t-s,t))e^{-\mu(t-s)}B(s)\,\mbox{d}s+\int_{0}^{t}\beta(S(t-s,t))e^{-\mu(t-s)}B(s)\,\mbox{d}s$
$\displaystyle=:f(t)+\int_{0}^{t}\beta(S(t-s,t))e^{-\mu(t-s)}B(s)\,\mbox{d}s$
$\displaystyle\leq
f(t)+\int_{0}^{t}\beta(x_{m}+g(0)(t-s))e^{-\mu(t-s)}B(s)\,\mbox{d}s.$

The kernel $k(a)=\beta(x_{m}+g(0)a)e^{-\mu a}$ of the linear Volterra integral
equation

$y(t)=f(t)+\int_{0}^{t}k(t-s)y(s)ds$ (5.3)

has a nonnegative resolvent $r$ (meaning that
$r(t)=k(t)+\int_{0}^{t}k(t-s)r(s)ds$ and $y(t)=f(t)+\int_{0}^{t}r(t-s)f(s)ds$)
(see Theorem 2.3.4 in [14]). Then by a generalized Grönwall lemma, one obtains
$B(t)\leq y(t)$ where $y(t)$ is the solution of (5.3).
Indeed, using the usual notation for convolution, (5.2) can be written as
$B\leq f+k*B$ and so $B=f-g+k*B$ for $g=f+k*B-B\geq 0.$ Then we have

$B=f-g+r*(f-g)=f+r*f-(g+r*g)=y-(g+r*g)\;\Rightarrow\;B\leq y,$

since $r$ and $g$ are nonnegative (cf. Theorem 9.8.2 in [14]). The claim
follows since $y(t)$ tends exponentially to $0$ when $R_{0}<1$ by Theorem 3.12
in [11] and the final part of the proof of Theorem 5.1.

Let us recall the notation

$\bar{S}(a)=x_{m}+\int_{0}^{a}g\left(\bar{B}\frac{e^{-\mu\tau}}{\mu}\right)\,\mbox{d}\tau,$
(5.4)

for the size of an individual of age $a$ at the non-trivial equilibrium.

Let us now compute the linearisation of (2.9) around the nontrivial
equilibrium $\bar{B}$ using (5.4). For this we set $B(t)=\bar{B}+y(t)$ and
write (formally)

$\begin{array}[]{ll}\bar{B}+y(t)\\\
=\displaystyle\int_{0}^{\infty}\left(\beta\big{(}\bar{S}(a)\big{)}+\beta^{\prime}\big{(}\bar{S}(a)\big{)}\int_{0}^{a}g^{\prime}\left(\bar{B}\frac{e^{-\mu\tau}}{\mu}\right)\right.\\\
\left.\quad\quad\quad\times\displaystyle\int_{\tau}^{\infty}e^{-\mu\alpha}y_{t}(-a+\tau-\alpha)\mbox{d}\alpha\mbox{d}\tau+o(y_{t})\right)e^{-\mu
a}\big{(}\bar{B}+y_{t}(-a)\big{)}\mbox{d}a,\end{array}$ (5.5)

which, using the steady state condition (4.1) and neglecting higher order
terms, leads to

$\displaystyle y(t)=$ $\displaystyle
D\mathcal{F}(\bar{B})y_{t}=\int_{0}^{\infty}\beta\big{(}\bar{S}(a)\big{)}\,e^{-\mu
a}y(t-a)\mbox{d}a$ (5.6)
$\displaystyle+\int_{0}^{\infty}\beta^{\prime}\big{(}\bar{S}(a)\big{)}\,e^{-\mu
a}\bigg{(}\,\int_{0}^{a}\bar{B}g^{\prime}\left(\bar{B}\frac{e^{-\mu\tau}}{\mu}\right)\int_{\tau}^{\infty}e^{-\mu\alpha}y(t-a+\tau-\alpha)\,\mbox{d}\alpha\mbox{d}\tau\bigg{)}\,\mbox{d}a$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\beta\big{(}\bar{S}(a)\big{)}e^{-\mu
a}y(t-a)\mbox{d}a+\int_{0}^{\infty}\beta^{\prime}\big{(}\bar{S}(a)\big{)}e^{-\mu
a}$
$\displaystyle\quad\times\left(\int_{0}^{a}\bar{B}g^{\prime}\big{(}\bar{B}\frac{e^{-\mu\tau}}{\mu}\big{)}e^{-\mu(-a+\tau)}\int_{a}^{\infty}e^{-\mu\sigma}y(t-\sigma)\mbox{d}\sigma\,\mbox{d}\tau\right)\mbox{d}a.$

######  Remark 5.3.

See Appendix A for a rigorous derivation of (5.6). There, $\mathcal{F}$ is
written essentially as the composition of simpler functions and then the chain
rule is applied.

Changing the order of integration, the expression within parentheses inside
the last integral in (5.6) can be rewritten as:

$\displaystyle\int_{0}^{a}\bar{B}g^{\prime}\left(\bar{B}\frac{e^{-\mu\tau}}{\mu}\right)\int_{a}^{\infty}e^{-\mu(-a+\tau+\sigma)}y(t-\sigma)\mbox{d}\sigma\mbox{d}\tau$
(5.7) $\displaystyle=$
$\displaystyle\int_{a}^{\infty}\int_{0}^{a}\bar{B}g^{\prime}\left(\bar{B}\frac{e^{-\mu\tau}}{\mu}\right)e^{-\mu\tau}\mbox{d}\tau
e^{-\mu(\sigma-a)}y(t-\sigma)\mbox{d}\sigma$ $\displaystyle=$
$\displaystyle\int_{a}^{\infty}\left(g\left(\frac{\bar{B}}{\mu}\right)-g\left(\bar{B}\frac{e^{-\mu
a}}{\mu}\right)\right)e^{-\mu(\sigma-a)}y(t-\sigma)\mbox{d}\sigma.$

Thus, changing the integration order again, the second term on the right hand
side of (5.6) reads

$\int_{0}^{\infty}\int_{0}^{\sigma}\beta^{\prime}(\bar{S}(a))\left[g\left(\frac{\bar{B}}{\mu}\right)-g\left(\bar{B}\frac{e^{-\mu
a}}{\mu}\right)\right]\mbox{d}a\,e^{-\mu\sigma}\,y(t-\sigma)\mbox{d}\sigma.$

Hence, (5.6) is of the form $y(t)=\int_{0}^{\infty}k(a)y(t-a)\mbox{d}a$ with
the kernel

$k(a)=\beta\big{(}\bar{S}(a)\big{)}\,e^{-\mu a}+e^{-\mu
a}\int_{0}^{a}\beta^{\prime}(\bar{S}(\alpha))\left[g\left(\frac{\bar{B}}{\mu}\right)-g\left(\frac{\bar{B}e^{-\mu\alpha}}{\mu}\right)\right]\,\mbox{d}\alpha.$

Since

$\displaystyle\int_{0}^{a}\beta^{\prime}(\bar{S}(\alpha))g\left(\frac{\bar{B}e^{-\mu\alpha}}{\mu}\right)\mbox{d}\alpha=\int_{0}^{a}\beta^{\prime}(\bar{S}(\alpha))\bar{S}^{\prime}(\alpha)\mbox{d}\alpha$
$\displaystyle=$
$\displaystyle\beta(\bar{S}(a))-\beta(\bar{S}(0))=\beta(\bar{S}(a))-\beta(x_{m})$

the kernel $k$ simplifies to

$k(a)=\beta(x_{m})e^{-\mu a}+g\left(\frac{\bar{B}}{\mu}\right)e^{-\mu
a}\int_{0}^{a}\beta^{\prime}(\bar{S}(\alpha))\mbox{d}\alpha,$

which leads to the characteristic equation

$1=\hat{k}(\lambda)=\frac{\beta(x_{m})}{\mu+\lambda}+\frac{1}{\mu+\lambda}g\left(\frac{\bar{B}}{\mu}\right)\int_{0}^{\infty}\beta^{\prime}(\bar{S}(a))e^{-(\mu+\lambda)a}\,\mathrm{d}a.$
(5.8)

Without the loss of generality, we will assume in the rest of this section
that the minimum size is $x_{m}=0$.

###  5.1 Stability of non-trivial steady states for a non-trivial example

We will assume in the following that the per capita fertility is given by:
$\beta(s)=\beta_{0}\max(0,s-x_{A}),$ where $x_{A}\geq 0$ is the adult size at
which individuals start to reproduce. Let us define $\bar{a}$ by

$\int_{0}^{\bar{a}}g\left(\frac{\bar{B}\exp(-\mu\tau)}{\mu}\right)\mbox{d}\tau=x_{A},$
(5.9)

i.e., $\bar{a}$ is the age at which individuals begin to reproduce given the
environmental condition associated to the equilibrium.

We have from (4.1) and $x_{m}=0$ that

$\displaystyle R(B)=$
$\displaystyle\int_{0}^{\infty}\beta\left(\int_{0}^{a}g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\,\mbox{d}\tau\right)\,\,e^{-\mu
a}\,\mbox{d}a$ (5.10) $\displaystyle=$
$\displaystyle\beta_{0}\int_{0}^{\infty}\max\left(0,\int_{0}^{a}g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\mbox{d}\tau-
x_{A}\right)e^{-\mu a}\mbox{d}a$ $\displaystyle=$
$\displaystyle\beta_{0}\int_{\bar{a}}^{\infty}\left(\int_{0}^{a}g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\mbox{d}\tau-\int_{0}^{\bar{a}}g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\mbox{d}\tau\right)e^{-\mu
a}\mbox{d}a$ $\displaystyle=$
$\displaystyle\beta_{0}\int_{\bar{a}}^{\infty}\int_{\bar{a}}^{a}g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\mbox{d}\tau
e^{-\mu a}\mbox{d}a$ $\displaystyle=$
$\displaystyle\beta_{0}\int_{\bar{a}}^{\infty}\int_{\tau}^{\infty}e^{-\mu
a}\mbox{d}a\,g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\mbox{d}\tau$
$\displaystyle=$
$\displaystyle\beta_{0}\int_{\bar{a}}^{\infty}\frac{e^{-\mu\tau}}{\mu}g\left(\frac{Be^{-\mu\tau}}{\mu}\right)\mbox{d}\tau=\frac{\beta_{0}}{\mu^{2}}\int_{0}^{e^{-\mu\bar{a}}}g\left(\frac{B}{\mu}\zeta\right)d\zeta,$

and hence, $R_{0}=R(0)=\frac{\beta_{0}g(0)}{\mu^{2}}e^{-\mu\bar{a}}.$

On the other hand, the characteristic equation (5.8) reduces to

$1=\beta_{0}\frac{g\left(\frac{\bar{B}}{\mu}\right)}{\lambda+\mu}\int_{\bar{a}}^{\infty}e^{-(\lambda+\mu)a}\mbox{d}a=\beta_{0}\frac{g(\bar{B}/\mu)}{(\mu+\lambda)^{2}}e^{-(\lambda+\mu)\bar{a}}=\frac{\mu^{2}}{(\mu+\lambda)^{2}}R_{0}\frac{g(\bar{B}/\mu)}{g(0)}e^{-\lambda\bar{a}},$
(5.11)

which, in the special case $x_{A}=0$ (or, equivalently, $\bar{a}=0$) allows to
identify the (two) roots as

$\lambda=-\mu\pm\sqrt{\beta_{0}g\left(\bar{B}/\mu\right)}=\mu\left(-1\pm\sqrt{R_{0}\frac{g(\bar{B}/\mu)}{g(0)}}\,\right).$
(5.12)

So, under this assumption, we are able to explicitly formulate the
characteristic equation and even to explicitly compute its roots. From the
condition for the existence of a nontrivial equilibrium (4.1) and (5.10) with
$\bar{a}=0,$ we have

$1=R(\bar{B})=\frac{\beta_{0}}{\mu^{2}}\int_{0}^{1}g\left(\frac{\bar{B}}{\mu}\zeta\right)d\zeta>\frac{\beta_{0}}{\mu^{2}}\min_{\zeta\in[0,1]}g\left(\frac{\bar{B}}{\mu}\zeta\right)=\frac{\beta_{0}}{\mu^{2}}g\left(\frac{\bar{B}}{\mu}\right)$

which implies that both eigenvalues given by (5.12) are negative. Hence,
Theorem 3.3 ensures that under these hypotheses the nontrivial steady state is
always asymptotically stable. We next show that also for $x_{A}>0$, if the
individual growth rate is decreasing, the non-trivial equilibrium is locally
asymptotically stable whenever it exists, excluding the possibility of a Hopf
bifurcation.

######  Theorem 5.4.

Let $x_{m}=0,\,\,\beta(s)=\beta_{0}\max\\{0,s-x_{A}\\},$ $R_{0}>1$ and let $g$
be decreasing. Then the nontrivial steady state is locally asymptotically
stable.

###### Proof.

We show that for all $x_{A}>0$ the characteristic equation has no purely
imaginary roots, preventing the presence of Hopf bifurcations and, via a
continuity argument, extending the result of the case $x_{A}=0$ to the case
$x_{A}>0.$ Indeed, we put $\lambda=i\omega$ with $\omega\in\mathbb{R}$ in
(5.11) to obtain

$(\mu+i\omega)^{2}=\beta_{0}g(\bar{B}/\mu)e^{-(\mu+i\omega)\bar{a}},$

which, taking real and imaginary parts, leads to

$\begin{array}[]{rcl}\mu^{2}-\omega^{2}=c\cos(\omega\bar{a})\\\ \\\
2\omega\mu=-c\sin(\omega\bar{a})\end{array}$ (5.13)

with $c=\beta_{0}g(\bar{B}/\mu)e^{-\mu\bar{a}}.$ Notice that (4.1), (5.10) and
(5.13) rule out that $\omega$ can be $0$. Solving the second equation for $c$
and substituting the result in the first, we obtain a quadratic equation for
$\mu$, with $\omega>0$ as a free parameter:

$\mu^{2}+2\omega\frac{\cos(\omega\bar{a})}{\sin(\omega\bar{a})}\mu-\omega^{2}=\left(\mu+\omega\cot\left(\frac{\omega\bar{a}}{2}\right)\right)\left(\mu-\omega\tan\left(\frac{\omega\bar{a}}{2}\right)\right)=0.$

Choosing the first root,

$\mu=-\omega\cot\left(\frac{\omega\bar{a}}{2}\right),$ (5.14)

we find

$c=-\frac{2\omega\mu}{\sin(\omega\bar{a})}=\frac{2\omega^{2}}{2\sin\left(\frac{\omega\bar{a}}{2}\right)\cos\left(\frac{\omega\bar{a}}{2}\right)}\frac{\cos\left(\frac{\omega\bar{a}}{2}\right)}{\sin\left(\frac{\omega\bar{a}}{2}\right)}=\left(\frac{\omega}{\sin\left(\frac{\omega\bar{a}}{2}\right)}\right)^{2},$
(5.15)

whereas the second one gives

$c=-\frac{2\omega\mu}{\sin(\omega\bar{a})}=-\frac{2\omega^{2}}{2\sin\left(\frac{\omega\bar{a}}{2}\right)\cos\left(\frac{\omega\bar{a}}{2}\right)}\frac{\sin\left(\frac{\omega\bar{a}}{2}\right)}{\cos\left(\frac{\omega\bar{a}}{2}\right)}=-\left(\frac{\omega}{\cos\left(\frac{\omega\bar{a}}{2}\right)}\right)^{2},$

which is incompatible with $c$ being positive. Now, using (5.15), the
definition of $c$, the fact that $g$ is decreasing, (5.10), the condition of
steady state $R(B)=1$ and (5.14) we have

$\left(\frac{\omega}{\sin\left(\frac{\omega\bar{a}}{2}\right)}\right)^{2}=c=\beta_{0}g\left(\frac{\bar{B}}{\mu}\right)e^{-\mu\bar{a}}<\beta_{0}\int_{0}^{e^{-\mu\bar{a}}}g\left(\frac{\bar{B}}{\mu}\zeta\right)d\zeta=\mu^{2}=\left(\omega\cot\left(\frac{\omega\bar{a}}{2}\right)\right)^{2},$

which implies $1<\left(\cos(\frac{\omega\bar{a}}{2})\right)^{2},$ impossible
for any real $\omega$. So roots cannot enter the right half plane by crossing
the imaginary axis. As the right hand side of (5.11) tends to zero for
$|\lambda|\rightarrow\infty$ when $\text{Re}\,\lambda\geq 0$, they cannot
enter from infinity either. Essentially from Rouché’s Theorem it now follows
that there are no roots in the right half plane for arbitrary $x_{A}>0$, see
Lemma XI.2.8 in [8] or Lemma 9.17.4 in [12].

###  5.2 Semi-explicit expression for a particular case

In this section we assume that the minimum size is $x_{m}=0$ and that the per
capita fertility is proportional to the size: $\beta(s)=\beta_{0}s$ (i.e.,
$x_{A}=0$). In addition we consider that the individual growth rate is of the
form $g(z)=\frac{g_{0}}{1+z/z_{0}}$ where $g_{0}>0$ and $z_{0}>0$ (recall that
$z$ represents the environment that an individual experiences, which is given
by the number of individuals that are larger than it).

In this situation, (5.10) gives

$R(B)=\frac{\beta_{0}}{\mu^{2}}\int_{0}^{1}\frac{g_{0}}{1+\frac{B}{\mu}\frac{\zeta}{z_{0}}}\mbox{d}\zeta=\frac{\beta_{0}g_{0}}{\mu^{2}}\frac{\ln\left(1+B/(\mu
z_{0})\right)}{B/(\mu z_{0})}=R_{0}\frac{\ln\left(1+B/(\mu
z_{0})\right)}{B/(\mu z_{0})}.$

Therefore, the birth rate at the nontrivial equilibrium (which necessarily
exists if $R_{0}>1$ as discussed in Section 4) is the unique positive solution
$\bar{B}$ of the equation $\frac{\ln\left(1+B/(\mu z_{0})\right)}{B/(\mu
z_{0})}=\frac{1}{R_{0}}.$ This allows an explicit expression for $\bar{B}$ in
terms of the Lambert function $W_{-1}$ as

$\bar{B}=\mu
z_{0}\left(-R_{0}W_{-1}\bigg{(}-\frac{\exp(-1/R_{0})}{R_{0}}\bigg{)}-1\right).$
(5.16)

Indeed, take $z=-(1+B/(\mu z_{0}))/R_{0}<-1/R_{0}$ in the preceding equation,
which gives $ze^{z}=-(1/R_{0})e^{-1/R_{0}}$. Then the (only) solution to this
equation is the Lambert function $W_{-1}$ ( i.e., the inverse function of the
(monotonously decreasing) function $f(z)=ze^{z}$ restricted to the interval
$(-\infty,-1)$) evaluated at $-(1/R_{0})e^{-1/R_{0}}$.

More interestingly, an explicit expression can also be obtained for the
density with respect to size in the steady state. Indeed, (5.4) gives in this
case,

$\bar{S}(a)=\int_{0}^{a}\frac{g_{0}}{1+\frac{\bar{B}e^{-\mu\tau}}{\mu
z_{0}}}\mbox{d}\tau=\frac{g_{0}}{\mu}\ln\left(\frac{\mu z_{0}e^{\mu
a}+\bar{B}}{\mu z_{0}+\bar{B}}\right),$ (5.17)

which leads to

$\bar{S}^{-1}(x)=\frac{1}{\mu}\ln\left(\frac{(\mu
z_{0}+\bar{B})e^{\frac{\mu}{g_{0}}x}-\bar{B}}{\mu z_{0}}\right)$

and to

$\bar{S}^{\prime}\big{(}\bar{S}^{-1}(x)\big{)}=g_{0}\frac{(\mu
z_{0}+\bar{B})\exp{\left(\frac{\mu}{g_{0}}x\right)}-\bar{B}}{(\mu
z_{0}+\bar{B})\exp{\left(\frac{\mu}{g_{0}}x\right)}}.$

By (4.3) we finally obtain

$\bar{u}(x)=\frac{\bar{B}\exp(-\mu\bar{S}^{-1}(x))}{\bar{S}^{\prime}\big{(}\bar{S}^{-1}(x)\big{)}}=\frac{\mu
z_{0}\bar{B}(\mu
z_{0}+\bar{B})\exp\left(\frac{\mu}{g_{0}}x\right)}{g_{0}\Big{(}\left(\mu
z_{0}+\bar{B}\right)\exp\left(\frac{\mu}{g_{0}}x\right)-\bar{B}\Big{)}^{2}}.$
(5.18)

Moreover, an easy integration gives the following expression for the
population number above an individual of size $x$

$\int_{x}^{\infty}\bar{u}(s)\mbox{d}s=\frac{\bar{B}z_{0}}{(\mu
z_{0}+\bar{B})e^{\frac{\mu}{g_{0}}x}-\bar{B}}.$ (5.19)

##  6 Concluding remarks

The principle of linearised stability (PLS for short), widely used in the
theory of ODEs, says that the stability of a stationary state is determined by
the stability properties of the linearised semigroup. This principle has also
been proved to hold in dynamical systems of infinite dimension with a
“semilinear” structure (namely semilinear PDEs and DE, see [15, 21, 8]) via
the variation of constants formula. In this article we used the PLS to analyse
rigorously the local stability of stationary birth rates of (2.11). As a
consequence of such an analysis we found that for reasonable and rather
general biological functional responses (see the hypotheses of Theorem 5.4),
the non-trivial stationary birth rate of (2.11) is locally asymptotically
stable.

The PLS, as stated above, cannot be applied to the PDE formulation presented
in Appendix B. The reason is that, as explained in detail in [3], the
nonlinear semigroup associated to (B.1) is not differentiable, and hence it
cannot be linearised. This does not mean, however, that, if the PDE system
(B.1) is linearised “formally” around a stationary distribution $\bar{u}$, the
stability of $\bar{u}$ can’t be determined from the stability of the
linearised system. In fact we expect that such is possible, but a proof is, as
far as we know, is still missing.

As explained in [3], a way to prove this result would be to establish an
“equivalence” between orbits of the delay formulation (in the state space of
weighted birth rate histories, i.e. $\mathcal{X}$) and orbits of the PDE
formulation (in the state space of integrable functions of height, i.e.
$L^{1}(x_{m},\infty)$). By an “equivalence” we specifically mean to find a
continuous function
$\mathcal{L}_{\text{DE}}^{\text{PDE}}:\mathcal{X}\rightarrow
L^{1}(x_{m},\infty)$ mapping orbits in $\mathcal{X}$ to orbits in
$L^{1}(x_{m},\infty)$ and vice-versa (i.e. an analogous continuous function
$\mathcal{L}_{\text{PDE}}^{\text{DE}}:L^{1}(x_{m},\infty)\rightarrow\mathcal{X}$),
so that stability results can be translated from one formulation to the other.
In [3] we found that for these functions to exist, one needed to work in a
(exponentially) weighted space of integrable functions of height,
$L^{1}_{w}(x_{m},\infty)$, where the proper value of $w$ depended on the
weight $\rho$ chosen for $\mathcal{X}$ (working with the unweighted space
$L^{1}(x_{m},\infty)$ was possible if $\rho$ was chosen to be equal to the
mortality rate $\mu$, since that implied $w=0$). In fact, in that paper the
phase spaces for both the PDE and the DE included a component with information
on the environmental condition. These additional components allowed to
establish a surjective function (with the desired properties mentioned above)
mapping states from the delay formulation to states of the PDE formulation
(and vice-versa by taking a pseudoinverse of that function). As we are about
to see, the analogous function associated to the (simpler) phase spaces used
in this paper fails to be surjective (precluding any attempt of extending the
results of [3] to the present work).

Natural candidates for $\mathcal{L}_{\text{DE}}^{\text{PDE}}$ and
$\mathcal{L}_{\text{PDE}}^{\text{DE}}$ may be obtained using the biological
interpretation of the functions involved in (2.11) and (B.1). Indeed, take
$\phi\in\mathcal{X}$ a birth rate history and $u_{0}\in L^{1}(x_{m},\infty)$ a
‘corresponding’ population height-distribution and define

$\displaystyle X(\tau;\phi):=S(-\tau,0;\phi)$
$\displaystyle=x_{m}+\int_{0}^{-\tau}g\left(\int_{\sigma}^{\infty}\phi(\tau+\sigma-\alpha)e^{-\mu\alpha}d\alpha\right)d\sigma$
(6.1)
$\displaystyle=x_{m}+\int_{0}^{-\tau}g\left(\int_{-\infty}^{\tau}\phi(\theta)e^{-\mu(\tau+\sigma-\theta)}d\theta\right)d\sigma$

for $\tau\in(-\infty,0]$ (i.e. the size at time 0 of an individual born at
$\tau$ given the birth rate history $\phi$, see (2.8)) and $T(x;\phi)$, for
$x\in[x_{m},\infty)$, as the inverse of $X(\cdot;\phi)$ (which exists if $g$
is bounded and decreasing and gives the time at birth of an individual with
size $x$ at time 0 given the birth rate history $\phi$). Then we have

$\int_{x_{m}}^{x}u_{0}(x)dx=\int_{T(x;\phi)}^{0}\phi(\theta)e^{\mu\theta}d\theta$
(6.2)

because being younger means being smaller, and hence the individuals smaller
than $x$ must coincide with those born after $T(x;\phi)$ that have survived.
Then, differentiation with respect to $x$ gives

$u_{0}(x)=-\phi(T(x;\phi))e^{\mu T(x;\phi)}T^{\prime}(x;\phi),$ (6.3)

which gives a natural candidate for $\mathcal{L}_{\text{DE}}^{\text{PDE}}$.
Similarly, by rewriting (6.2) as

$\int_{x_{m}}^{X(\tau;\phi)}u_{0}(x)dx=\int_{\tau}^{0}\phi(\theta)e^{\mu\theta}d\theta,$

differentiation with respect to $\tau$ gives

$u_{0}(X(\tau;\phi))X^{\prime}(\tau;\phi)=-\phi(\tau)e^{\mu\tau}.$ (6.4)

Unlike (6.3), the above equation is problematic in that it does not give an
explicit formula for $\phi$ in terms of $u_{0}$. It turns out that the above
equation does not define implicitly $\phi\in\mathcal{X}$ for each $u_{0}\in
L^{1}(x_{m},\infty)$ (which is equivalent to say that
$\mathcal{L}_{\text{DE}}^{\text{PDE}}$ defined through (6.3) is not
surjective). To see this choose, as a counterexample, $\mu=0$, $g(E)=1-E$ for
$E<1/2$ (it doesn’t matter what $g$ does for $E\geq 1/2$, besides being
decreasing) and $u_{0}(x)=1$ for $x\in(x_{m},x_{m}+1)$ and 0 otherwise. Then
formula (6.1) simplifies to $X(\tau;\phi)=x_{m}-\tau
g\left(\int_{-\infty}^{\tau}\phi(\theta)d\theta\right)$ and equation (6.4)
implies

$\phi(\tau)=\frac{g\left(\int_{-\infty}^{\tau}\phi(\theta)d\theta\right)}{1-\tau
g^{\prime}\left(\int_{-\infty}^{\tau}\phi(\theta)d\theta\right)}$

if $X(\tau;\phi)<1$ and $\phi(\tau)=0$ otherwise. This forces the support of
$\phi$ to be $(-1,0)$, so that $X(\tau;\phi)<x_{m}+1$ for $\tau\in(-1,0)$, and
thus $\phi$ solves (6.4) only if it satisfies

$\phi(\tau)=\frac{1-\int_{-1}^{\tau}\phi(\theta)d\theta}{1+\tau}$

as long as $\int_{-1}^{\tau}\phi(\theta)d\theta<1/2$. Since the right hand
side of this equation has a non-integrable singularity for $\tau\downarrow-1$,
this relation contradicts that $\phi\in\mathcal{X}$. The fact that equation
(6.4) fails to define a birth rate history in $\mathcal{X}$ as a function of
$u_{0}$ means that there are reasonable population densities with respect to
size (such as the indicator function used in the example) that cannot be
obtained by prescribing an integrable birth rate history.

As already mentioned, this situation deviates from what we had in [3], where
an explicit formula for $\mathcal{L}_{\text{PDE}}^{\text{DE}}$ was derived
thanks to the additional environmental variable that was considered as part of
the phase space (and somehow provided more room to play with). Since the
scalar renewal equation presented in Section 2 was obtained precisely by
expressing the environmental variable in terms of the birth rate history (and
thus restricting the set of admissible environmental histories), a way to
overcome this difficulty would be to work with an extended version of the
delay formulation in which the environmental history is a proper element of
the phase space (and thus there is also a delay equation for it). In addition,
such an extended version would allow us to analyse more general environmental
feedbacks. For instance environmental feedbacks of the form
$E(x,t)=\int_{x}^{\infty}\alpha(y)u(y,t)dy$ (compare with (2.1)), where the
impact of larger individuals depends on their size. Such situations cannot be
formulated in terms of only a renewal equation for the birth rate. Indeed,
since the environmental history is needed to give the size individuals will
have in the future, the environmental condition felt by an individual is no
longer determined only by the individuals born before him but it depends also
on the environmental history itself. The drawback of an extended formulation
is that then the environmental history $t\mapsto E(\cdot,t)$ takes values in
an infinite dimensional space, which makes the analysis of the
differentiability (analogue of Theorem 3.1) much more involved (the theory to
deal with these cases is developed in [10]).

What could be the implications of such a non-equivalence between the two
formulations, and specifically of the fact that there are population densities
that cannot be obtained naturally from a birth rate history? It seems that the
non-equivalence does not imply differences in the number of stationary states
and attractors in general found in each formulation. In fact we expect a one-
to-one correspondence between orbits in the $\omega$-limit sets of the two
formulations (such a correspondence would be a consequence of the relation
between solutions of the RE and solutions of the PDE given in subsection B.1
of Appendix B). What might be affected by the non-equivalence is the stability
behaviour of the corresponding $\omega$-limit sets. A priori (with what we
have shown in this paper) we cannot rule out the possibility that a stationary
population density of the PDE formulation is unstable, while the corresponding
stationary birth rate of the delay formulation is stable. The reason is that
there are states arbitrarily close to such a stationary population density
that cannot be related to any birth history from a neighbourhood of the
stationary birth rate. Further work is needed to rule out this kind of
discrepancy between the two formulations (or, alternatively, to give a
specific example where the discrepancy takes place, although we doubt that
such an example exists).

##  Appendix A Differentiability

######  Theorem A.1.

(Theorem 3.1) Assume that $g$ and $\beta$ have a bounded and globally
Lipschitzian first derivative (with common constant $2C$). Also assume that
$g$ is bounded, positive and bounded away from $0$. Then the map
$\mathcal{F}:\mathcal{X}\rightarrow\mathbb{R}$ defined in (3.1) is
continuously differentiable with bounded derivative provided that the
parameter $\rho$ in the definition of $\mathcal{X}$ satisfies $\rho<\mu/5.$

###### Proof.

First notice that the hypotheses imply the following estimate for any $z\geq
0$ and $h>-z$:

$\begin{array}[]{lr}|g(z+h)-g(z)-g^{\prime}(z)h|\\\
=\left|\int_{z}^{z+h}g^{\prime}(s)ds-g^{\prime}(z)h\right|\leq\left|\int_{z}^{z+h}|g^{\prime}(s)-g^{\prime}(z)|\,ds\right|\leq
2C\left|\int_{z}^{z+h}|z-s|ds\right|\leq Ch^{2}\end{array}$ (A.1)

and analogously for $\beta$.

The statement of the theorem amounts to showing that

$\phi\rightarrow(\tilde{\mathcal{F}}(\phi))(a)=e^{-\mu
a}\,\beta\bigg{(}x_{m}+\int_{0}^{a}g\big{(}e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\phi(-s)ds\big{)}\,d\tau\bigg{)}$

is a continuously differentiable map from the positive cone of the Banach
space

$\mathcal{X}=\left\\{\phi\in
L_{loc}^{1}(-\infty,0):\left|\left|\phi\right|\right|_{\mathcal{X}}:=\int_{-\infty}^{0}e^{\rho
s}|\phi(s)|ds<\infty\right\\}$

to its dual identified with the Banach space

$\mathcal{X}^{\prime}=\left\\{f\in
L_{loc}^{\infty}(0,\infty):||f||_{\mathcal{X}^{\prime}}:=esssup_{a\in[0,\infty)}e^{\rho
a}|f(a)|<\infty\right\\}$

with the duality product $\langle
f,\phi\rangle=\int_{0}^{\infty}f(a)\phi(-a)\mbox{d}a.$

Indeed, we can write
$\mathcal{F}(\phi)=\langle\tilde{\mathcal{F}}(\phi),\phi\rangle,$ and a rather
general and straightforward argument gives, assuming differentiability of
$\tilde{\mathcal{F}},$

$D\mathcal{F}(\phi)\psi=\langle\tilde{\mathcal{F}}(\phi),\psi\rangle+\langle\
D\tilde{\mathcal{F}}(\phi)\psi,\phi\rangle.$ (A.2)

In particular, for $\phi=0,$ we have
$D\mathcal{F}(0)\psi=\langle\tilde{\mathcal{F}}(0),\psi\rangle.$

Next we define three intermediate spaces of real valued continuous functions:

$Y=\left\\{P\in C(T):||P||_{Y}:=\sup_{(\tau,a)\in T}e^{-\rho
a}|P(\tau,a)|<\infty\right\\}$

where $T=\\{(\tau,a)\in\mathbb{R}^{2}:0\leq\tau\leq a<\infty\\},$

$Z=\left\\{v\in C(T):||v||_{Z}:=\sup_{(\tau,a)\in
T}e^{-\rho_{1}a}|v(\tau,a)|<\infty\right\\}$

with $\rho_{1}>0$ to be chosen later,

$W=\left\\{S\in
C([0,\infty)):||S||_{W}:=\sup_{a\in[0,\infty)}e^{-\rho_{2}a}|S(a)|<\infty\right\\}$

with $\rho_{2}>0$ to be chosen later; and four maps:

$\mathcal{L}_{1}:\mathcal{X}\rightarrow Y\,\,\text{defined
by}\,\,(\mathcal{L}_{1}\phi)(\tau,a)=e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\phi(-s)\mbox{d}s,$ $\mathcal{G}:Y\rightarrow Z\,\,\text{defined
by}\,\,\mathcal{G}(P)=g\circ P,$ $\mathcal{L}_{2}:Z\rightarrow
W\,\,\text{defined
by}\,\,(\mathcal{L}_{2}v)(a)=x_{m}+\int_{0}^{a}v(\tau,a)\mbox{d}{\color[rgb]{0,0,1}\tau}$

and

$\mathcal{B}:W_{+}\rightarrow\mathcal{X}^{\prime}\,\,\text{defined
by}\,\,\mathcal{B}(S)(a)=e^{-\mu a}\,(\beta\circ S)(a),$

($W_{+}$ meaning the positive cone of $W$) in such a way that (at least
formally)
$\tilde{\mathcal{F}}=\mathcal{B}\circ\mathcal{L}_{2}\circ\mathcal{G}\circ\mathcal{L}_{1}.$
Then the claim will follow from the chain rule provided we prove that the four
maps are well defined and continuously differentiable with bounded derivative.

Step 1. $\mathcal{L}_{1}$ is bounded linear provided that $\rho\leq\mu$.

We have

$\sup_{(\tau,a)\in T}e^{-\rho a}\left|e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\phi(-s)ds\right|=\sup_{(\tau,a)\in T}e^{-\rho
a}\left|e^{-\mu(\tau-a)}\int_{-\infty}^{-a}e^{\mu s}\phi(s)ds\right|$

$\leq\sup_{a\geq 0}\int_{-\infty}^{-a}e^{(\mu-\rho)(a+s)}e^{\rho
s}|\phi(s)|ds\leq\int_{-\infty}^{0}e^{\rho s}|\phi(s)|ds,$

since $(\mu-\rho)(a+s)\leq 0$ in the last but one integral. Thus,
$||\mathcal{L}_{1}\phi||_{Y}\leq||\phi||_{\mathcal{X}}.$

Step 2. $\mathcal{G}$ is continuously differentiable with bounded derivative
provided that $2\rho\leq\rho_{1}$.
$\mathcal{G}$ is well defined because $g$ is bounded and continuous.
Let $P\in Y$ and $Q\in Y$ such that $||Q||_{Y}=1$, which implies
$|Q(\tau,a)|\leq e^{\rho a}.$
We start by proving that $Q\rightarrow g^{\prime}(P(\cdot))Q(\cdot)$ defines a
bounded linear map $\mathcal{Y}\rightarrow\mathcal{Z}$ with norm bounded
independently of $P$:

$\displaystyle\sup_{\left|\left|Q\right|\right|_{Y}=1}\left|\left|g^{\prime}(P(\cdot))Q(\cdot)\right|\right|_{Z}$
$\displaystyle=\sup_{\left|\left|Q\right|\right|_{Y}=1}\sup_{z\in
T}e^{-\rho_{1}a}\left|g^{\prime}(P(z))Q(z)\right|$
$\displaystyle\leq\left|\left|g^{\prime}\right|\right|_{\infty}\sup_{z\in
T}e^{-\rho a}\left|Q(z)\right|=\left|\left|g^{\prime}\right|\right|_{\infty}.$

Moreover, we can write, setting $z=(\tau,a),$ and using (A.1),

$\displaystyle e^{-\rho_{1}a}|g(P(z)+\varepsilon
Q(z))-g(P(z))-g^{\prime}(P(z))\varepsilon Q(z)|$ $\displaystyle\leq
C\varepsilon^{2}e^{-\rho_{1}a}|Q(z)|^{2}$ $\displaystyle\leq
C\varepsilon^{2}e^{(2\rho-\rho_{1})a}\leq C\varepsilon^{2},$

i.e.,

$\left|\left|\mathcal{G}(P+\varepsilon
Q)-\mathcal{G}(P)-g^{\prime}(P(\cdot))\varepsilon
Q(\cdot)\right|\right|_{Z}\leq C\varepsilon^{2}.$

Therefore, $(D\mathcal{G}(P)Q)(z):=g^{\prime}(P(z))Q(z)$ is the Fréchet
derivative of $\mathcal{G}$ at the point $P$, its norm is uniformly bounded by
$\left|\left|g^{\prime}\right|\right|_{\infty}$; and it is (uniformly)
continuous: for $Q\in Y$ with norm $1$ we have

$\displaystyle\left|\left|D\mathcal{G}(P_{1})Q-D\mathcal{G}(P_{2})Q\right|\right|_{Z}=\sup_{z\in
T}e^{-\rho_{1}a}\left|\big{(}g^{\prime}(P_{1}(z))-g^{\prime}(P_{2}(z))\big{)}Q(z)\right|$
$\displaystyle\leq 2C\sup_{z\in
T}e^{-\rho_{1}a}|P_{1}(z)-P_{2}(z)|\,|Q(z)|\leq 2C\sup_{a\geq
0}e^{(-\rho_{1}+2\rho)a}||P_{1}-P_{2}||_{Y}\leq 2C||P_{1}-P_{2}||_{Y}.$

Step 3. $\mathcal{L}_{2}$ is a positive continuous affine map provided that
$0<\rho_{1}<\rho_{2}$.
It suffices to see,

$\begin{array}[]{ll}\left|\left|\mathcal{L}_{2}v-x_{m}\right|\right|_{W}=\sup_{a\in[0,\infty)}e^{-\rho_{2}a}\left|\int_{0}^{a}v(\tau,a)d\tau\right|\\\
\leq\sup_{a\in[0,\infty)}e^{-\rho_{2}a}\int_{0}^{a}e^{\rho_{1}a}||v||_{Z}\,\mbox{d}\tau=\sup_{a\in[0,\infty)}ae^{(-\rho_{2}+\rho_{1})a}||v||_{Z}\leq\frac{1}{e(\rho_{2}-\rho_{1})}||v||_{Z}.\end{array}$

Step 4.
$\mathcal{B}$ is continuously differentiable with bounded derivative provided
that $\rho+2\rho_{2}\leq\mu$.
First notice that the assumptions on $\beta$ imply that there exist positive
constants $C_{1}$ and $C_{2}$ such that $\beta(s)\leq C_{1}+C_{2}s.$ Thus
$\mathcal{B}$ is well defined: for $S\in W_{+},$ since $|S(a)|\leq
e^{\rho_{2}a}||S||_{W},$

$\displaystyle e^{\rho a}|e^{-\mu a}\beta(S(a))|$ $\displaystyle\leq
e^{(\rho-\mu)a}(C_{1}+C_{2}|S(a)|)$ $\displaystyle\leq
e^{(\rho-\mu)a}(C_{1}+C_{2}e^{2\rho_{2}a}||S||_{W})\leq C_{1}+C_{2}||S||_{W}.$

As in Step 2, let us prove that $R\rightarrow
e^{-\mu\cdot}\beta^{\prime}(S(\cdot))R(\cdot)$ defines a bounded linear map
$W\rightarrow\mathcal{X^{\prime}}$ with norm bounded independently of $S$:

$\displaystyle\sup_{\left|\left|R\right|\right|_{W}=1}\left|\left|e^{-\mu\cdot}\beta^{\prime}(S(\cdot))R(\cdot)\right|\right|_{\cal{X}^{\prime}}=$
$\displaystyle\sup_{\left|\left|R\right|\right|_{W}=1}\sup_{a\geq 0}e^{\rho
a}\left|e^{-\mu a}\beta^{\prime}(S(a))R(a)\right|$ $\displaystyle\leq$
$\displaystyle\sup_{a\geq
0}e^{(\rho+\rho_{2}-\mu)a}\left|\left|\beta^{\prime}\right|\right|_{\infty}\leq\left|\left|\beta^{\prime}\right|\right|_{\infty}.$

Let us now proceed to show that $\mathcal{B}$ is differentiable: Let $S\in
W_{+}$ and $R\in W$ with norm equal to 1, which implies
$|R(a)|<e^{\rho_{2}a}$. Then, for $\epsilon$ small enough, $S+\epsilon R\in
W_{+}$. Then we can write, using the $\beta$-variant of (A.1),

$\begin{array}[]{ll}e^{\rho a}\left|e^{-\mu a}\beta(S(a)+\varepsilon
R(a))-e^{-\mu a}\beta(S(a))-e^{-\mu a}\beta^{\prime}(S(a))\varepsilon
R(a)\right|\\\ \leq C\varepsilon^{2}e^{(\rho-\mu+2\rho_{2})a}\leq
C\varepsilon^{2},\end{array}$

proving that $\big{(}D\mathcal{B}(S)R\big{)}(a):=e^{-\mu
a}\beta^{\prime}(S(a))R(a)$ is the Fréchet derivative of $\mathcal{B}$ at the
point $S,$ with norm uniformly bounded by
$\left|\left|\beta^{\prime}\right|\right|_{\infty}$
We also show that the derivative is continuous as in Step 2. Let $R\in W$ with
norm equal to 1. We have

$\displaystyle\left|\left|D\mathcal{B}(S_{1})R-D\mathcal{B}(S_{2})R\right|\right|_{\mathcal{X}^{\prime}}=$
$\displaystyle\sup_{a\in[0,\infty)}e^{(\rho-\mu)a}\left|\left(\beta^{\prime}(S_{1}(a))-\beta^{\prime}(S_{2}(a))\right)R(a)\right|$
$\displaystyle\leq$ $\displaystyle
2C\sup_{a\in[0,\infty)}e^{(\rho-\mu)a}|S_{1}(a)-S_{2}(a)|\,|R(a)|$
$\displaystyle\leq$ $\displaystyle 2C\sup_{a\geq
0}e^{(\rho-\mu+2\rho_{2})a}||S_{1}-S_{2}||_{W}\leq 2C||S_{1}-S_{2}||_{W}.$

Finally, given any $\rho\in(0,\frac{\mu}{5})$ we can take $\rho_{1}=2\rho$
(fulfilling the assumption of Step 2) and $\rho_{2}=\frac{\mu-\rho}{2}>2\rho$
(fulfilling the assumption of Step 3 and that of Step 4 since then
$2\rho_{2}+\rho=\mu)$ to conclude the proof.

As a consequence, the chain rule gives, taking into account that
$\mathcal{L}_{1}$ is linear and $\mathcal{L}_{2}$ is affine,

$D\tilde{\mathcal{F}}(\phi)\psi=D(\mathcal{B}\circ\mathcal{L}_{2}\circ\mathcal{G}\circ\mathcal{L}_{1})(\phi)\psi=D\mathcal{B}(\mathcal{L}_{2}\mathcal{G}(\mathcal{L}_{1}\phi))\,(\mathcal{L}_{2}-x_{m})\,D\mathcal{G}(\mathcal{L}_{1}\phi)\,\mathcal{L}_{1}\psi.$

Since we are interested in linearisation around steady states, we can restrict
to evaluation of the differential on constant functions $\bar{B}.$ So, we
compute, sequentially:

$\displaystyle\mathcal{L}_{1}\psi\,(\tau,a)=$ $\displaystyle
e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu s}\psi(-s)\mbox{d}s,$
$\displaystyle\mathcal{L}_{1}\bar{B}\,(\tau,a)=$
$\displaystyle\bar{B}e^{-\mu\tau}/\mu,$ $\displaystyle
D\mathcal{G}(\mathcal{L}_{1}\bar{B})\mathcal{L}_{1}\psi\,(\tau,a)=$
$\displaystyle
g^{\prime}(\bar{B}e^{-\mu\tau}/\mu)e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\psi(-s)\mbox{d}s,$
$\displaystyle(\mathcal{L}_{2}-x_{m})\,D\mathcal{G}(\mathcal{L}_{1}\bar{B})\mathcal{L}_{1}\psi\,(\tau,a)=$
$\displaystyle\int_{0}^{a}g^{\prime}(\bar{B}e^{-\mu\tau}/\mu)e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\psi(-s)\mbox{d}s\,\mbox{d}\tau$ $\displaystyle=$ $\displaystyle:h(a),$

and, also,

$\mathcal{L}_{2}\,\mathcal{G}(\mathcal{L}_{1}\bar{B})=x_{m}+\int_{0}^{a}g(\bar{B}e^{-\mu\tau}/\mu)\mbox{d}\tau(=\bar{S}(a)),$

where, in the last equality we assumed, furthermore, that $\bar{B}$ is not
only a constant function, but a steady state (see (5.4)). Therefore,

$\displaystyle D\tilde{\mathcal{F}}(\bar{B})\psi=$ $\displaystyle
D\mathcal{B}(\mathcal{L}_{2}\mathcal{G}(\mathcal{L}_{1}\bar{B}))\,(\mathcal{L}_{2}-x_{m})\,D\mathcal{G}(\mathcal{L}_{1}\bar{B})\,\mathcal{L}_{1}\,\psi(a)$
$\displaystyle=$ $\displaystyle
D\mathcal{B}(\mathcal{L}_{2}\,\mathcal{G}(\mathcal{L}_{1}\bar{B}))h(a)=e^{-\mu
a}\beta^{\prime}(\bar{S}(a))h(a)$ $\displaystyle=$ $\displaystyle e^{-\mu
a}\beta^{\prime}(\bar{S}(a))\,\int_{0}^{a}g^{\prime}(\bar{B}e^{-\mu\tau}/\mu)e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\psi(-s)\mbox{d}s\,\mbox{d}\tau.$

Finally, we will have,

$\langle\
D\tilde{\mathcal{F}}(\bar{B})\psi,\bar{B}\rangle=\int_{0}^{\infty}e^{-\mu
a}\beta^{\prime}(\bar{S}(a))\,\int_{0}^{a}g^{\prime}(\bar{B}e^{-\mu\tau}/\mu)e^{-\mu(\tau-a)}\int_{a}^{\infty}e^{-\mu
s}\psi(-s)\mbox{d}s\,\mbox{d}\tau\bar{B}\mbox{d}a,$

which, together with (A.2), gives (5.6).

##  Appendix B The PDE formulation

The classical formulation derived by imposing a conservation law leads to the
(non-local, quasi-linear and first-order) partial differential equation

$\displaystyle\frac{\partial}{\partial t}u(x,t)+\frac{\partial}{\partial
x}\left(g(E(x,t))u(x,t)\right)+\mu u(x,t)$ $\displaystyle=0,$ (B.1)
$\displaystyle g(E(x_{m},t))u(x_{m},t)$
$\displaystyle=\int_{x_{m}}^{\infty}\beta(y)u(y,t)\,\mathrm{d}y,$
$\displaystyle E(x,t)$ $\displaystyle=\int_{x}^{\infty}u(y,t)\,\mathrm{d}y.$

Here the second equation stands for the flux of newborns, offspring of
individuals of any size $y$ which have a size specific per capita fertility
(obviously nonnegative) $\beta(y)$. Notice that the fertility is indeed
indirectly affected by negative density dependence since a larger value of the
environmental variable leads to a smaller size achieved by the individuals.
From a dynamical point of view the solutions of (B.1) can be seen as orbits
$t\mapsto u(\cdot,t)$ in the space of integrable functions with respect to
height, i.e. $L^{1}(x_{m},\infty)$.

The slightly more general model with environmental interaction variable

$E(x,t)=\alpha\int_{0}^{x}u(y,t)\,\mathrm{d}y+\int_{x}^{M}u(y,t)\,\mathrm{d}y,\quad\alpha\in[0,1],$

(but with finite maximal size $M$) was studied for example in [5, 1, 2, 13],
and a very general model incorporating distributed recruitment in [4]. In [16]
the well posedness of the above problem was proven by rewriting the system in
terms of characteristic coordinates.

In this appendix we include a series of results showing that the PDE
formulation is tightly related to the delay formulation (as it should be since
both models are built from a description of the same biological processes). In
subsection B.1 we show that one can solve the PDE problem by solving a scalar
RE (with integration from 0 to $t$) for the population birth rate $B$ and that
the large time limiting form of this equation is exactly (2.9). In subsection
B.2 we show that the condition characterising the existence of non-trivial
steady states of (B.1) coincides with (4.1) (in addition a formula for the
non-trivial stationary population size-density is given). Finally in
subsection B.3 we show that the formal linearisation of system (B.1) leads to
the characteristic equation (5.8).

###  B.1 Solution of the PDE in terms of a renewal equation

The solution of (B.1) can be written as the sum of two terms: the first
considers the individuals born between $0$ and $t$ and the second considers
the individuals that already exist at time 0, i.e. those reflected in the
initial population density $u_{0}(x)$.

First notice that at time 0,

$\bar{E}(\xi)=\int_{\xi}^{\infty}u_{0}(\eta)d\eta$

gives the number of individuals with size larger than $\xi$, while at time
$\tau$

$\tilde{E}(\tau)=\left(\int_{0}^{\tau}B(\sigma)e^{\mu\sigma}d\sigma+\int_{0}^{\infty}u_{0}(\eta)d\eta\right)e^{-\mu\tau}$

gives the number of individuals with size larger than $x_{m}$. Since the
mortality rate is constant, these numbers decrease exponentially with rate
$\mu$ as time increases. As a consequence, the size at time $t$ of an
individual with size $\xi$ at time 0 is

$X(t,0,\xi)=\xi+\int_{0}^{t}g(\bar{E}(\xi)e^{-\mu\sigma})d\sigma$

and the size at time $t$ of an individual born at time $\tau>0$ with
$0<\tau<t$ is

$X(t,\tau,x_{m})=x_{m}+\int_{0}^{t-\tau}g\left(\tilde{E}(\tau)e^{-\mu\sigma}\right)d\sigma.$

So the birth rate has to satisfy the renewal equation

$B(t)=B_{\text{dsc}}(t)+B_{\text{fnd}}(t)$

where

$B_{\text{dsc}}(t)=\int_{0}^{t}\beta(X(t,\tau,x_{m}))B(\tau)e^{-\mu(t-\tau)}d\tau$

is the birth rate associated to the descendants of the founder population and

$B_{\text{fnd}}(t)=\int_{0}^{t}\beta(X(t,0,\xi))u_{0}(\xi)d\xi\,e^{-\mu t}$

is the known birth rate associated to the founder population. Once we solve
the renewal equation constructively, we can obtain an explicit expression for
the (weak) solution of the PDE by integrating along characteristics.

Note that $B_{\text{fnd}}(t)$ tends to 0 exponentially as
$t\rightarrow\infty$. By changing $\tau$ to $a$ with $t-\tau=a$ we can rewrite

$B_{\text{dsc}}(t)=\int_{0}^{t}\beta(X(t,t-a,x_{m}))B(t-a)e^{-\mu a}da.$

Now note that

$X(t,t-a,x_{m})=x_{m}+\int_{0}^{a}g(\tilde{E}(t-a)e^{-\mu\tau})d\tau$

and

$\tilde{E}(t-a)=\int_{0}^{t-a}B(\eta)e^{\mu\eta}d\eta\,e^{-\mu(t-a)}+\int_{0}^{\infty}u_{0}(\eta)d\eta\,e^{-\mu(t-a)}$

where the second summand at the right hand side tends exponentially to 0 as
$t\rightarrow\infty$. Since this term represents the founder population that
remains at time $t$, let us refer to it as $P_{\text{fnd}}(t)$. Next, by using
the transformation $\eta=t-s$ we have

$\tilde{E}(t-a)=\int_{a}^{t}B(t-s)e^{-\mu(s-a)}ds+P_{\text{fnd}}(t)=e^{\mu
a}\int_{a}^{t}B(t-s)e^{-\mu s}ds+P_{\text{fnd}}(t).$

Now note that, by ignoring $B_{\text{fnd}}(t)$ and $P_{\text{fnd}}(t)$ and by
replacing the upper integration boundary $t$ in the last integral by $\infty$,
we obtain (2.11).

###  B.2 Existence and characterization of non-trivial steady states

To establish criteria for the existence of non-trivial steady states $\bar{u}$
in the PDE formulation is apparently more complex than what we had to do for
the delay formulation in Section 4.

Let us first concentrate on the ordinary differential equation which arises
from the first and the third equations in (B.1) when one assumes that
$\bar{u}$ only depends on $x$. This leads to the following second order ode
for $E(x):=\int_{x}^{\infty}\bar{u}(s)ds$,

$\frac{d}{dx}\left(g(E(x))E^{\prime}(x)\right)+\mu E^{\prime}(x)=0,$

or, equivalently, to

$g(E(x))E^{\prime}(x)+\mu E(x)=C,$

for some constant $C.$ Since $E(x)$ tends to $0$ when $x$ tends to $\infty,$
$C$ has to coincide with (minus) the flux of individuals leaving the system at
infinity:
$C=\lim_{x\rightarrow\infty}g(E(x))E^{\prime}(x)=-g(0)\lim_{x\rightarrow\infty}\bar{u}(x)$
and so it has to be $0$ (since otherwise
$\lim_{x\rightarrow\infty}\bar{u}(x)=-C/g(0)\neq 0$ and $\bar{u}$ would not be
integrable).
Therefore we look for solutions of the differential equation

$\frac{dE}{dx}(x)=-\mu\frac{E(x)}{g(E(x))}$

with initial condition $E(x_{m})=N$ (the total population) and such that
$\lim_{x\rightarrow\infty}E(x)=0.$ Equivalently,

$\int_{E(x)}^{N}\frac{g(z)}{z}\mathrm{d}z=\mu(x-x_{m}).$

If $G$ is a primitive of $g(z)/z,$ the previous equation reads

$G(N)-G(E(x))=\mu(x-x_{m}),$

which, can be rewritten as

$E(x)=G^{-1}\left(G(N)-\mu(x-x_{m})\right).$

It follows that

$\displaystyle\bar{u}(x)=$
$\displaystyle-\frac{d}{dx}\left(G^{-1}\left(G(N)-\mu(x-x_{m})\right)\right)$
(B.2) $\displaystyle=$
$\displaystyle\frac{\mu}{G^{\prime}\left(G^{-1}\left(G(N)-\mu(x-x_{m})\right)\right)}=\mu\frac{G^{-1}\left(G(N)-\mu(x-x_{m})\right)}{g\left(G^{-1}\left(G(N)-\mu(x-x_{m})\right)\right)}.$

Since $\bar{u}(x_{m})=\frac{\mu N}{g(N)}$ we have

$g(E(x_{m}))\bar{u}(x_{m})=g(N)\bar{u}(x_{m})=\mu N.$

Therefore, using the boundary condition, a non-trivial steady state (given by
(B.2)) does exist if and only if a positive number $N$ exists such that

$N=\int_{x_{m}}^{\infty}\beta(x)\frac{G^{-1}\left(G(N)-\mu(x-x_{m})\right)}{g\left(G^{-1}\left(G(N)-\mu(x-x_{m})\right)\right)}\mathrm{d}{x}.$
(B.3)

This turns out to be equivalent to (4.1) with $N=B/\mu$. Indeed, we can write

$\displaystyle R(B)=$
$\displaystyle\int_{0}^{\infty}\beta\left(x_{m}+\int_{0}^{a}g\left(B\frac{e^{-\mu\tau}}{\mu}\right)\,\mbox{d}\tau\right)\,\,e^{-\mu
a}\,\mbox{d}a$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\beta\left(x_{m}+\int_{Ne^{-\mu
a}}^{N}\frac{g\left(z\right)}{\mu z}\,\mbox{d}z\right)\,\,e^{-\mu
a}\,\mbox{d}a$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\beta\left(x_{m}+\frac{G(N)-G(Ne^{-\mu
a})}{\mu}\right)e^{-\mu a}\,\mbox{d}a$ $\displaystyle=$
$\displaystyle\frac{1}{N}\int_{x_{m}}^{\infty}\beta(x)\frac{G^{-1}(G(N)-\mu(x-x_{m}))}{g\left(G^{-1}(G(N)-\mu(x-x_{m})\right)}\mbox{d}x$

where in the second equality we performed the change of variables
$z=B\frac{e^{-\mu\tau}}{\mu}$, and in the fourth one, the change of variables
$x=x_{m}+\frac{G(N)-G(Ne^{-\mu a})}{\mu}.$ See Section 5.2 where a particular
case is developed and where an explicit expression for a primitive $G$ is
available.

###  B.3 Linearisation in the PDE formulation

The (formal) linearisation of the PDE (B.1) around the steady state $u_{*}$ is
very economical as it simply reads (note that $g^{\prime}$ below stands for
the derivative of $g$ with respect to its argument $E$)

$\displaystyle
v_{t}(x,t)+\left(g(E_{*}(x))v(x,t)+g^{\prime}(E_{*}(x))u_{*}(x)\int_{x}^{\infty}v(y,t)\,\mathrm{d}y\right)_{x}=$
$\displaystyle-\mu v(x,t),$ (B.4) $\displaystyle
g(E_{*}(x_{m}))v(x_{m},t)+g^{\prime}(E_{*}(x_{m}))u_{*}(x_{m})\int_{x_{m}}^{\infty}v(x,t)\,\mathrm{d}x=$
$\displaystyle\int_{x_{m}}^{\infty}\beta(x)v(x,t)\,\mathrm{d}x.$

Substituting $v(x,t)=e^{\lambda t}V(x)$ into (B.4) we have

$\displaystyle\left(g(E_{*}(x))V(x)+g^{\prime}(E_{*}(x))u_{*}(x)\int_{x}^{\infty}V(y)\,\mathrm{d}y\right)_{x}=$
$\displaystyle-(\lambda+\mu)V(x),$ (B.5) $\displaystyle
g(E_{*}(x_{m}))V(x_{m})+g^{\prime}(E_{*}(x_{m}))u_{*}(x_{m})\int_{x_{m}}^{\infty}V(x)\,\mathrm{d}x=$
$\displaystyle\int_{x_{m}}^{\infty}\beta(x)V(x)\,\mathrm{d}x.$

Therefore, $\lambda\in\mathbb{C}$ is an eigenvalue, if and only (B.5) admits a
solution $V\not\equiv 0$. We also note that although the size domain is
unbounded, it can be shown that the part of the spectrum of the semigroup
generator in the half plane
$\left\\{z\in\mathbb{C}\,|\,\text{Re}(z)>-\mu\right\\}$ contains only
eigenvalues , see e.g. [13, Sect.4.] for more details, and therefore (linear)
stability can indeed be characterized by the leading eigenvalue of the
semigroup generator.

For the trivial steady state $u_{*}\equiv 0$ the left hand side of (B.5) has
only local terms and therefore easily leads to the characteristic equation

$g(0)=\int_{x_{m}}^{\infty}\beta(x)e^{-\frac{\lambda+\mu}{g(0)}(x-x_{m})}dx,$

which is exactly what one gets by inserting $y(t)=e^{\lambda t}$ into (5.1)
and making the change of variables $x=x_{m}+g(0)a$. Therefore, the stability
of $u_{*}\equiv 0$ is characterized by the net reproduction number ($R$
evaluated at the zero steady state, or the virgin environment as we previously
referred to), as expected. That is, if

$R(0)=(R_{0}=)\displaystyle\int_{x_{m}}^{\infty}\frac{\beta(x)}{g(0)}e^{-\frac{\mu}{g(0)}(x-x_{m})}\,\mathrm{d}x>1,$

then $u_{*}\equiv 0$ is unstable; while $R(0)<1$ implies that the trivial
steady state is asymptotically stable.

To deduce the characteristic equation we integrate the first equation of (B.5)
from $x$ to $\infty$, to obtain

$g(E_{*}(x))V(x)+g^{\prime}(E_{*}(x))u_{*}(x)\int_{x}^{\infty}V(y)\,\mathrm{d}y=(\lambda+\mu)\int_{x}^{\infty}V(y)\,\mathrm{d}y.$
(B.6)

Substituting $x=x_{m}$ into (B.6) and combining it with the second equation in
(B.5) yields

$(\lambda+\mu)\int_{x_{m}}^{\infty}V(x)\,\mathrm{d}x=\int_{x_{m}}^{\infty}\beta(x)V(x)\,\mathrm{d}x.$
(B.7)

Note that $\lambda\in\mathbb{C}$ is an eigenvalue if and only if (B.6)-(B.7)
admits a solution $V\not\equiv 0$. To see for which $\lambda$ this is possible
let us introduce

$H(x):=\int_{x}^{\infty}V(y)\,\mathrm{d}y$ (B.8)

as unknown so that (B.6) boils down to differential equation

$-g(E_{*}(x))H^{\prime}(x)+g^{\prime}(E_{*}(x))u_{*}(x)H(x)=(\lambda+\mu)H(x),$

whose solution is

$H(x)=H(x_{m})\exp\left(\int_{x_{m}}^{x}\frac{g^{\prime}(E_{*}(r))u_{*}(r)-(\lambda+\mu)}{g(E_{*}(r))}\,\mathrm{d}r\right)=H(x_{m})\pi(x,\lambda),$
(B.9)

where we defined

$\pi(x,\lambda):=\exp\left(\int_{x_{m}}^{x}\frac{g^{\prime}(E_{*}(r))u_{*}(r)-(\lambda+\mu)}{g(E_{*}(r))}\,\mathrm{d}r\right),\quad\lambda\in\mathbb{C},\,\,x\in[x_{m},\infty).$
(B.10)

Then, substitution of (B.9) into (B.7) via (B.8) yields the characteristic
equation

$(\lambda+\mu)=-\int_{x_{m}}^{\infty}\beta(x)\frac{\partial}{\partial\,x}\pi(x,\lambda)\,\mathrm{d}x.$
(B.11)

This equation can be rewritten (assuming that $\beta$ is differentiable) as

$(\lambda+\mu)=-\beta(\infty)\pi(\infty,\lambda)+\beta(x_{m})+\int_{x_{m}}^{\infty}\beta^{\prime}(x)\pi(x,\lambda)\,\mathrm{d}x.$
(B.12)

Next note that since $E_{*}^{\prime}(x)=-u_{*}(x)$ we have

$\displaystyle\pi(x,\lambda)=$
$\displaystyle\exp\left(\int_{x_{m}}^{x}\frac{g^{\prime}(E_{*}(r))u_{*}(r)-(\lambda+\mu)}{g(E_{*}(r))}\,\mathrm{d}r\right)$
$\displaystyle=$
$\displaystyle\exp\left(-\int_{x_{m}}^{x}\frac{\lambda+\mu}{g(E_{*}(r)}\,\mathrm{d}r\right)\exp\left(-\int_{x_{m}}^{x}\frac{\frac{\mathrm{d}}{\mathrm{d}\,r}(g(E_{*}(r))}{g(E_{*}(r))}\,\mathrm{d}r\right)$
$\displaystyle=$
$\displaystyle\exp\left(-\int_{x_{m}}^{x}\frac{\lambda+\mu}{g(E_{*}(r)}\,\mathrm{d}r\right)\frac{g(E_{*}(x_{m}))}{g(E_{*}(x))}.$

Then, if $\mu>\displaystyle\sup_{x\geq
x_{m}}\left\\{g^{\prime}(E_{*}(x))u_{*}(x)\right\\}$ (for example if
$g^{\prime}\leq 0$) one has $\pi(\infty,\lambda)=0$ for every
$\lambda\in\mathbb{C}$, and the characteristic equation reduces to

$\lambda+\mu=\beta(x_{m})+\int_{x_{m}}^{\infty}\beta^{\prime}(x)\exp\left(-\int_{x_{m}}^{x}\frac{\lambda+\mu}{g(E_{*}(r)}\,\mathrm{d}r\right)\frac{g(E_{*}(x_{m}))}{g(E_{*}(x))}\,\mathrm{d}x.$
(B.13)

Now let us rewrite equation (B.13) such that the integration variable is age
$a$, which will show that it is the characteristic equation (5.8) in disguise.
Using that we have

$\frac{\mathrm{d}a}{\mathrm{d}x}(x)=\frac{1}{g(E_{*}(x))},\,\,\bar{S}(a)=\int_{0}^{a}g\left(\bar{B}\frac{e^{-\mu\tau}}{\mu}\right)\,\mathrm{d}\tau=\Gamma^{-1}(a)=x,\,\,\Gamma(x):=\int_{0}^{x}\frac{1}{g(E_{*}(r))}\,\mathrm{d}r,\,\,E_{*}(x)=\bar{B}\frac{e^{-\mu\Gamma(x)}}{\mu},$

equation (B.13) can be rewritten as

$\lambda+\mu=\beta(x_{m})+g\left(\frac{\bar{B}}{\mu}\right)\int_{0}^{\infty}\beta^{\prime}(\bar{S}(a))e^{-(\lambda+\mu)a}\,\mathrm{d}a,$

which is identical to the characteristic equation (5.8) that was deduced from
the delay formulation.

## Acknowledgements

This work was partially supported by the research projects MT2017-84214C2-2-P
and PID2021-123733NB-I00. We also thank the International Centre for
Mathematical Sciences for financial support we received from the Research in
Groups program during our stay at Edinburgh in July 2017. The first ideas for
the present manuscript arose there and then.

## References

  * [1] AS Ackleh and K Deng.  A quasilinear hierarchical size-structured model: well-posedness and approximation.  Applied Mathematics and Optimization 2005; 51:35-59.
  * [2] AS Ackleh and K Ito.  Measure-valued solutions for a hierarchically size-structured population.  Journal of Differential Equations 2005; 217:431-455.
  * [3] C Barril, A Calsina, O Diekmann and JZ Farkas.  On the Formulation of Size-Structured Consumer Resource Models (with Special Attention for the Principle of Linearised Stability)  _Mathematical Models and Methods in Applied Sciences_ , 32 (2022), 1141–1191.
  * [4] A Calsina and J Saldaña.  Basic theory for a class of models of hierarchically structured population dynamics with distributed states in the recruitment.  Mathematical Models and Methods in Applied Sciences 16 (2006), 1695-1722.
  * [5] JM Cushing.  The dynamics of hierarchical age-structured populations.  Journal of Mathematical Biology 32 (1994), 705-729.
  * [6] DA Dewsbury.  Dominance rank, copulatory behavior, and differential reproduction  _The Quarterly Review of Biology_ , 57 (1982), 135–159.
  * [7] O Diekmann, M Gyllenberg, H Huang, M Kirkilionis, JAJ Metz and HR Thieme.  On the formulation and analysis of general deterministic structured population models II. Nonlinear theory  _Journal of Mathematical Biology_ , 43 (2001), 157–189.
  * [8] O Diekmann, SA van Gils, SM Verduyn Lunel and HO Walther.  Delay equations: functional-, complex-, and nonlinear analysis. Springer-Verlag, New York, 1995.
  * [9] O Diekmann, M Gyllenberg, JAJ Metz, H. Thieme.  On the formulation and analysis of general deterministic structured population models I. Linear theory  _Journal of Mathematical Biology_ , 36 (4) (1998), 349-388.
  * [10] O Diekmann, M Gyllenberg.  Abstract delay equations inspired by population dynamics.  _Functional analysis and evolution equations_ , 187-200, Birkhäuser, Basel, 2007.
  * [11] O Diekmann, M Gyllenberg.  Equations with infinite delay: blending the abstract and the concrete  _Journal of Differential Equations_ , 252 (2012), no. 2, 819-851.
  * [12] J Dieudonne.  Foundations of modern analysis. Academic Press, New York, 1960.
  * [13] JZ Farkas, T Hagen.  Hierarchical size-structured populations: the linearized semigroup approach.  Dynamics of Continuous, Discrete & Impulsive Systems. Series A. Mathematical Analyis 2010; 17:639-657.
  * [14] G Gripenberg, S-O Londen, O Staffans.  Volterra integral and functional equations.  _Encyclopedia of Mathematics and its Applications_ , Cambridge U P, 1990\.
  * [15] D Henry.  Geometric theory of semilinear parabolic equations. Lecture Notes in Mathematics, 840, Springer-Verlag, 1981.
  * [16] EA Kraev.  Existence and uniqueness for height structured hierarchical population models.  Natural Resource Modeling 2001; 14:45–70.
  * [17] P Magal and S Ruan.  _Structured Population Models in Biology and Epidemiology_ ,  Lecture Notes in Mathematics, 1936, Springer, 2008.
  * [18] P Magal and Z Zhang.  Competition for light in forest population dynamics: From computer simulator to mathematical model  _Journal of Theoretical Biology_ , 419 (2012), 290-304.
  * [19] JAJ Metz and O Diekmann.  The dynamics of physiologically structured populations.  _Lecture Notes in Biomathematics_ , 68, 1986.
  * [20] AM de Roos and L Persson.  Population and Community Ecology of Ontogenetic Development.  Monographs in Population Biology 51, Princeton University Press, 2013.
  * [21] A Pazy.  Semigroups of linear operators and applications to partial differential equations.  Springer-Verlag, New York, 1983.
  * [22] GA Polis.  The evolution and dynamics of intraspecific predation.  _Annual Review of Ecology and Systematics_ , 12 (1981), 225–251.
  * [23] GF Webb.  Theory of nonlinear age-dependent population dynamics.  Monographs and textbooks in pure and applied mathematics 89, 1985.
  * [24] MA Zavala, O Angulo, RB de la Parra and JC López-Marcos.  An analytical model of stand dynamics as a function of tree growth, mortality and recruitment: The shade tolerance-stand structure hypothesis revisited  _Journal of Theoretical Biology_ , 244 (2007), 440-450.