Boundary perturbations and steady states of structured populations

TL;DR

This paper develops boundary perturbation techniques to analyze steady states in structured population models, enabling control of spectral bounds and recovering known reproductive numbers for age-structured models.

Contribution

It introduces new boundary perturbation results that help control spectral bounds in nonlinear structured population models, especially those with fertility-driven dynamics.

Findings

Established conditions for positive steady states.

Controlled spectral bounds via boundary perturbations.

Reproduced classical reproductive numbers using new methods.

Abstract

In this work we establish conditions which guarantee the existence of (strictly) positive steady states of a nonlinear structured population model. In our framework the steady state formulation amounts to recasting the nonlinear problem as a family of eigenvalue problems combined with a fixed point problem. Amongst other things our formulation requires us to control the growth behaviour of the spectral bound of a family of linear operators along positive rays. For the specific class of model we consider here, this presents a considerable challenge. We are going to show that the spectral bound of the family of operators, arising from the steady state formulation, can be controlled by perturbations in the domain of the generators (only). These new boundary perturbation results are particularly important for models exhibiting fertility controlled dynamics. As an important by-product of the application of the boundary perturbation results we employ here we recover (using a recent theorem by H.R. Thieme) the familiar net reproduction number (or function) for models with single state at birth, which include for example the classic McKendrick (linear) and Gurtin-McCamy (non-linear) age-structured models.

Equations

$$\frac{du}{dt}={\mathcal{A}}_{\mathbf{u}}u,D({\mathcal{A}}_{\mathbf{u}})\subseteq \mathcal{X},u(0)=u_0.$$

$$E:\mathcal{X}\to \mathcal{Y},E(u)=\mathbf{u}.$$

$${\mathcal{A}}_{\mathbf{u}}u=0u,\Psi (E(u))=\mathbf{u};0\not\equiv u\in {\mathcal{X}}_{+},$$

$$\mu \left(a,\frac{{\int }_a^{m}u(r,t)\mathrm{d}r}{{\int }_0^{m}u(r,t)\mathrm{d}r}\right),\beta \left(a,{\int }_0^{m}u(r,t)\mathrm{d}r\right).$$

$$\mu \left(\cdot ,\frac{{\int }_{\cdot }^{m}u(r,t)\mathrm{d}r}{{\int }_0^{m}u(r,t)\mathrm{d}r}\right)\equiv \mu \left(\cdot ,\frac{{\int }_{\cdot }^m\alpha u(r,t)\mathrm{d}r}{{\int }_0^m\alpha u(r,t)\mathrm{d}r}\right),\forall \alpha >0.$$

$$p_t(a,t)+p_a(a,t)=$$

$$p(0,t)=$$

$$F_1(p)(\cdot )=\frac{{\int }_{\cdot }^{m}p(r)\mathrm{d}r}{{\int }_0^{m}p(r)\mathrm{d}r},p\in {\mathcal{X}}_{+}\setminus \{0\},F_2(p)={\int }_0^{m}p(r)\mathrm{d}r,p\in {\mathcal{X}}_{+}.$$

$${\mathcal{B}}_{\mathbf{u}}p=$$

$$D({\mathcal{B}}_{\mathbf{u}})=$$

$${\Phi }_{\mathbf{u}}p=$$

$$\beta (m,\cdot )\ne 0.$$

$$S=\{\mathbf{u}\in {\mathcal{Y}}_{+}^0|s({\mathcal{B}}_{\mathbf{u}})=0\},$$

$$||x|{|}_{-1}:={\left|\left|{\mathcal{A}}^{-1}x\right|\right|}_{\mathcal{X}}.$$

$${\mathcal{A}}_{|}z:=\mathcal{A}z,D({\mathcal{A}}_{|}):=\left\{z\in D(\mathcal{A})\cap \mathcal{Z}|\mathcal{A}z\in \mathcal{Z}\right\}.$$

$${\hat{\mathcal{B}}}_{\mathbf{u}}p=-p^{\prime }-\mu \left(\cdot ,F_1(\mathbf{u})(\cdot )\right)p,D({\hat{\mathcal{B}}}_{\mathbf{u}})=\left\{p\in W^{1,1}(0,m)|p(0)=0\right\},$$

$${\hat{\mathcal{C}}}_{\mathbf{u}}p=-{\Phi }_{\mathbf{u}}(p){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}\left(\exp \left\{-{\int }_0^{\cdot }\mu (r,F_1(\mathbf{u})(r))\mathrm{d}r\right\}\right),$$

$$\pi (a):=\pi \left(a,F_1(\mathbf{u})\right):=\exp \left\{-{\int }_0^a\mu (r,F_1(\mathbf{u})(r))\mathrm{d}r\right\},a\in [0,m].$$

$$D_1:=\left\{y\in \mathcal{X}|\left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y\in \mathcal{X}\right\}=\left\{y\in W^{1,1}(0,m)|y(0)={\Phi }_{\mathbf{u}}(y)\right\}=:D_2.$$

$$y_n(x):=$$

$${\pi }_n(x):=$$

$${\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{y}_n\to {\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}y,\text{and}-{\Phi }_{\mathbf{u}}(y_n){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{\pi }_n\to -{\Phi }_{\mathbf{u}}(y){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}\pi ={\hat{\mathcal{C}}}_{\mathbf{u}}y,$$

$${\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{y}_n-{\Phi }_{\mathbf{u}}(y_n){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{\pi }_n\to \left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y,\text{in}{\mathcal{X}}_{-1}.$$

$${\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{y}_n-{\Phi }_{\mathbf{u}}(y_n){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{\pi }_n\to \left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y=-y^{\prime }-\mu y,\text{in}\mathcal{X}.$$

$${\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{y}_n-{\Phi }_{\mathbf{u}}(y_n){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{\pi }_n=-y_n^{\prime }-\mu y_n+{\Phi }_{\mathbf{u}}(y_n)\left({\pi }_n^{\prime }+\mu {\pi }_n\right)=-y^{\prime }-\mu y,$$

$${\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{y}_n$$

$$=$$

$${\left|\left|{\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{y}_n-{\Phi }_{\mathbf{u}}(y_n){\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}{\pi }_n-(-y^{\prime }-\mu y)\right|\right|}_{\mathcal{X}}\le {\int }_0^{\frac{1}{n}}|y^{\prime }+\mu y|\mathrm{d}x$$

$$+{\int }_0^{\frac{1}{n}}\left|{\Phi }_{\mathbf{u}}(y_n)n\pi \left(\frac{1}{n}\right)-ny\left(\frac{1}{n}\right)\right|\mathrm{d}x+{\int }_0^{\frac{1}{n}}\left|nx\left({\Phi }_{\mathbf{u}}(y_n)\mu \pi \left(\frac{1}{n}\right)-\mu y\left(\frac{1}{n}\right)\right)\right|\mathrm{d}x$$

$$\le {\int }_0^{\frac{1}{n}}|y^{\prime }+\mu y|\mathrm{d}x+\left|{\Phi }_{\mathbf{u}}(y_n)\pi \left(\frac{1}{n}\right)-y\left(\frac{1}{n}\right)\right|+\frac{\sup (\mu )}{2n}\left|{\Phi }_{\mathbf{u}}(y_n)\pi \left(\frac{1}{n}\right)-y\left(\frac{1}{n}\right)\right|$$

$$\underset{n\to \infty }{\to }0,$$

$$\left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y_n\underset{n\to \infty }{\to }\left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y=z\in \mathcal{X},$$

$${\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}^{-1}\left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y_n$$

$$=y_n-{\Phi }_{\mathbf{u}}(y_n)\pi \to {\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}^{-1}z=:w\in D({\hat{\mathcal{B}}}_{\mathbf{u}}),\text{in}\mathcal{X};$$

$$\left({\left({\hat{\mathcal{B}}}_{\mathbf{u}}\right)}_{-1}+{\hat{\mathcal{C}}}_{\mathbf{u}}\right)y$$

Full text

1. Recap and motivation

In recent years we have developed a general framework to treat steady state problems of nonlinear partial differential and partial integro-differential equations, see [6, 7, 8]. We have been mainly interested in studying models of physiologically structured populations (for reference on structured population models see for example the monographs [14, 27, 31, 37]). The types of nonlinearities naturally arising in these models present considerable challenges. For example a basic size-structured population model can be formulated as a quasi-linear first order hyperbolic partial differential equation with a nonlinear and non-local boundary condition. On the other hand, the method we have developed to treat the steady state problem capitalises on the fact that the models we have considered describe the dynamics of populations, and consequently in our framework we have to study families of positive linear operators, admitting desirable spectral properties. In particular, the operators that naturally arise generate positive strongly continuous semigroups, often with further desirable regularity properties, such as eventual compactness and irreducibility. The spectral theory of these operators is well developed (at least on abstract Lebesgue spaces), see for example the monographs [3, 10, 20, 33].

To start let us briefly recall from [7] how we can (re)formulate the positive steady state problem for a nonlinear evolution equation, in general. We are then going to show in the subsequent section how to apply this general framework to a specific class of model exhibiting fertility controlled dynamics. Let $\mathcal{X}$ and $\mathcal{Y}$ be Banach lattices, and let us consider the (parametrised) family of abstract Cauchy problems.

[TABLE]

Above, for every ${\bf u}\in\mathcal{Y}$, $\mathcal{A}_{\bf u}$ is a linear operator, typically, the generator of a strongly continuous semigroup on the state space $\mathcal{X}$. For many of the concrete applications we have considered (see e.g. [6, 7, 22]) $\mathcal{X}$ can be naturally chosen as the Lebesgue space of integrable functions; while in the simplest case (as for example in case of the model we have considered in [22]) we have $\mathcal{Y}=\mathbb{R}$.

The relationship between ${\bf u}\in\mathcal{Y}$-parameter space, and $u\in\mathcal{X}$-state space, is determined in general by the so-called environmental operator:

[TABLE]

The latter condition we call the environmental condition, or environmental feedback. From the biological point of view, $E$ determines how the standing population affects (via population level interaction) individual development, hence the terminology (see e.g. [18, 19] for more details). From the mathematical point of view, the parametrisation can often be chosen in such a way that $E$ is a positive linear operator, for example an integral operator. In particular, if problem (1.1)-(1.2) can be set up such that the range of $E$ is contained in $\mathbb{R}^{n}$ for some $n\in\mathbb{N}$, then we say that problem (1.1) incorporates $n$-dimensional nonlinearity. Otherwise we say that our model incorporates an infinite dimensional nonlinearity. We note that for models with infinite dimensional nonlinearities (e.g. the one we have considered in [6]) it is often desirable (from the mathematical point of view) to set $\mathcal{Y}=\mathcal{X}$. Also note that (in particular also in the case when $\mathcal{Y}=\mathcal{X}$), $E$ is not necessarily surjective.

The advantage of the abstract formulation above is that the (positive) steady state problem can be formally written as:

[TABLE]

where $\Psi:\mathcal{Y}_{+}\to\mathcal{Y}_{+}$ is an appropriately defined non-linear map. At the same time, (at least) in case of (structured) population models, the generators $\mathcal{A}_{\bf u}$ naturally determine well-posed positive linear evolution problems. There is a well-developed mathematical framework to analyse these linear evolution problems, for example using the theory of positive strongly continuous semigroups. Solving problem (1.3) (that is, establishing the existence of a positive steady state of our nonlinear model) then amounts to study the spectral behaviour of the family of operators $\{\mathcal{A}_{\bf u}\,|\,{\bf u}\in\mathcal{Y}\}$ (note that $\mathcal{Y}$ is an infinite dimensional vector space, in general), and to prove the existence of a fixed point of a (possibly set-valued and non-convex) nonlinear map, which is related to the feedback condition $\Psi(E(u))={\bf u}$ in (1.3). Both of these problems pose different types of challenges depending on the type of nonlinearity, and the structure of the model equations, we are dealing with. In the present paper we are going to tackle one of these important challenges through a concrete example, and therefore further develop, and at the same time also demonstrate, the applicability of the general theory outlined in [7].

It is evident that there are natural models of biological populations which give rise to equations incorporating different types of nonlinearities, in particular both finite and infinite dimensional ones, simultaneously. A case in point is a population exhibiting a hierarchical structure, in which the hierarchy affects individual mortality. More specifically we have a population in mind in which older (more experienced) individuals possess an advantage when competing for resources, for example food. At the same time we assume that individual fertility rate is modulated by the total population size, due to scramble competition effects.

As a motivating example consider the following age-specific mortality rate $\mu$, and fertility rate $\beta$.

[TABLE]

Above $a\in[0,m]$ denotes age, while $u(a,t)$ denotes the density of individuals of age $a$ at time $t$. Note that while the mortality of an individual of age $a$ depends on the proportion of older individuals in the population, the fertility of an individual of age $a$ depends on the total number of individuals in the population. In particular, the corresponding age-structured partial differential equation model (see the next section) describes the time evolution of a hierarchic age-structured population, where individuals in the population compete with older individuals for resources, e.g. food, affecting their mortality. Similar hierarchic age-structured population models were in fact already introduced and studied for example by Cushing in [13]. Importantly, here we assume that individual mortality is affected by the proportion of older individuals, while individual fertility (or the success of mating) is affected by the total population size. This in fact can be viewed as mortality is affected by pure contest competition (mass independent), while fertility is affected by scramble (mass dependent) competition. Realistic ecological situations may include populations exhibiting cannibalistic behaviour, where the death rate of individuals is modulated by the proportion of older (larger) individuals in the population. In general, the effects of frequency dependent survival probabilities have been studied in the case of specific real world (e.g. fish) populations, see e.g. [32].

It turns out however that incorporating certain types of nonlinearities (such as the ones above in (1.4)) in a structured population model, poses some challenges, when trying to apply the abstract framework (1.1)-(1.3) to study the existence of positive steady states. In particular, our steady state problem formulation in (1.3) requires, amongst other things, to control the spectral bound $s(\mathcal{A}_{\bf u})$ on positive rays in the parameter space $\mathcal{Y}$. But as we can see $\mu$ is in fact constant along such positive rays, since clearly

[TABLE]

Hence the information how the spectral bound of the generator $\mathcal{A}_{\bf u}$ changes with respect to perturbations of ${\bf u}$ along rays is encoded in the fertility function $\beta$. In general for distributed states at birth models (such as the one we considered in [22]), since the recruitment operator appears as a bounded perturbation of the main part of the generator, this does not pose additional problems. However, for single state at birth models, such as the one we consider here, the function $\beta$ only plays a role in the definition of the domain of the unbounded generator $\mathcal{A}_{\bf u}$, and as a consequence any useful information pertaining the growth (or decay) of its spectral bound can only be extracted from its domain $D(\mathcal{A}_{\bf u})$. Our main aim in this work is to show (in the next section) that we can deal with such domain perturbation problems by using the framework of Sobolev towers, as we have already indicated in [5].

As a by-product, using the Sobolev tower construction and the boundary perturbation results contained in [20], together with a recent result of H. R. Thieme from [36], we recover the familiar net reproduction number/function. In particular, we rigorously justify that the biologically plausible definition of the net reproduction number (see e.g. [27]) - for the linear McKendrick model -, and function (see e.g. [23]) - for the nonlinear Gurtin-MacCamy age-structured model -, is indeed the spectral radius of an appropriate operator, which arises when we apply the natural splitting to the generator $\mathcal{A}_{\bf u}$. Note that the idea (which we proposed earlier in [5], see in particular the Concluding Remarks section) to lift the age-structured problem into the extrapolated space is crucial, as this allows us to apply the necessary and natural splitting of the generator $\mathcal{A}_{\bf u}$ which is required to apply Thieme’s theorem and to arrive at the biologically plausible net reproduction number/function (see also [21] for more details in the direction of net reproduction functions).

2. Boundary perturbations - through a concrete example

Our general goal in this section is to show how the spectral bound of the family of operators $\{\mathcal{A}_{\bf u}\,|\,{\bf u}\in\mathcal{Y}_{+}\}$ can be controlled (along positive rays in $\mathcal{Y}$) by perturbations in the domain of $\mathcal{A}_{\bf u}$ only. We are going to demonstrate this through a concrete example, but the approach we present here is quite general and can be applied to a variety of models with non-local boundary conditions (single state at birth). Here as an example we consider a model of a hierarchic age-structured population with mortality and fertility rates as described in (1.4). In particular, with the specific choice of the vital rates in (1.4) the familiar (single state at birth) Gurtin-McCamy-type age-structured population model reads as follows.

[TABLE]

We set the state space as $\mathcal{X}=L^{1}(0,m)$, and we note that our model (2.5)-(2.6) is highly nonlinear. In particular the nonlinearities are induced by two essentially different type of competition effects (scramble and contest), which are incorporated in the model via the two operators $F_{1}$ and $F_{2}$, which are defined as follows.

[TABLE]

The results obtained for this model will therefore also illustrate that the general framework we have outlined in [7] can be applied to models incorporating multiple (different) types of competition effects (inducing different types of nonlinearities), simultaneously.

To formulate the (positive) steady state problem for model (2.5)-(2.6) in the form of (1.3), we set the parameter space as $\mathcal{Y}:=L^{1}(0,m)=\mathcal{X}$, and we define for every ${\bf u}\in\mathcal{Y}_{+}^{0}:=\mathcal{Y}_{+}\setminus\{{\bf 0}\}$ the linear operators $\mathcal{B}_{\bf u}$ and $\Phi_{\bf u}$ as follows.

[TABLE]

In (2.9) above $W^{1,1}(0,m)$ denotes the usual Sobolev space of integrable functions with generalized first derivative in $L^{1}$. For every ${\bf u}\in\mathcal{Y}_{+}^{0}$ it can be shown that $\mathcal{B}_{\bf u}$ generates a positive and eventually compact semigroup of operators, under some very mild smoothness assumptions on the model ingredients $\beta$ and $\mu$ (see e.g. [23] for more details). For example, we may, and we will, in the rest of the paper assume that $\mu$ and $\beta$ are non-negative, continuous and bounded, and that $\mu\not\equiv 0,\beta\not\equiv 0$.

Moreover, the semigroup generated by $\mathcal{B}_{\bf u}$ can be shown irreducible if we impose the natural assumption that

[TABLE]

To establish the existence of a positive steady state of model (2.5)-(2.6), using the general framework from [7], we need to assure (amongst other things) the existence of a ‘desirable’ set in the parameter space $\mathcal{Y}$, where the spectral bound of the operators $\mathcal{B}_{\bf u}$ vanish. In particular, our goal is to show that this level set, defined specifically as

[TABLE]

intersects every positive ray in $\mathcal{Y}$ in a unique point (element), and that it is contained in a ball $B\subset\mathcal{Y}$ of a sufficiently large radius. To this end, it is sufficient to show that the spectral bound is strictly monotone decreasing (or, alternatively, increasing) along positive rays, that is $\forall\,{\bf u}\in\mathcal{Y}^{0}_{+},$ we have $s\left(\mathcal{B}_{\bf u}\right)>s\left(\mathcal{B}_{\alpha\,{\bf u}}\right),\,\forall\,\alpha>1$; and in particular that it changes sign along every ray in the set $B\cap\mathcal{Y}_{+}$. But, as we noted before, $\mu$ is constant along positive rays (since $F_{1}$ is constant along rays in $\mathcal{Y}$), hence to establish that the spectral bound is strictly monotone along rays in the positive cone we need to control the spectral bound via the boundary operator $\Phi_{\bf u}$, which only appears in the domain of the generator $\mathcal{B}_{\bf u}$.

It is the main purpose of this paper to show how this can be done using the framework of Sobolev towers, in particular using the extrapolated space of $\mathcal{X}$. To the best of our knowledge the notion of extrapolated spaces was originally introduced to establish existence of solutions of some nonlinear evolution equations, see e.g. [16]. Significantly, later they were also employed to deal with boundary perturbations of generators of semigroups, see e.g. [17, 25] for more details. We also mention that there is a well-developed and fruitful theory of integrated semigroups, in particular with applications to age-structured population models. We refer the interested reader for example to [29, 30] for more details and further references.

Recall from [20, Ch.II] that the extrapolated space $\mathcal{X}_{-1}$ corresponding to the generator $\mathcal{A}$ (assuming it is invertible) is defined as the completion of $\mathcal{X}$ with respect to the following norm:

[TABLE]

Then, $\mathcal{X}$ is dense in $\mathcal{X}_{-1}$ and the generator $\mathcal{A}$ of the strongly continuous semigroup $\mathcal{T}$ has a unique continuous extension $\mathcal{A}_{-1}$ with $D(\mathcal{A}_{-1})=\mathcal{X}$, which is the generator of the so-called extended semigroup $\mathcal{T}_{-1}$ on $\mathcal{X}_{-1}$. Recall also from [20, Ch.II] that in general $\mathcal{A}_{|}$ is called the part of $\mathcal{A}$ in $\mathcal{Z}\subseteq\mathcal{X}$ if

[TABLE]

Note that $\mathcal{A}$ is the part of $\mathcal{A}_{-1}$ in $\mathcal{X}$.

In our setting we consider the modified generator

[TABLE]

and introduce the operator

[TABLE]

where $\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}$ is the unique continuous extension of $\hat{\mathcal{B}}_{\bf u}$ with domain $D\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}\right)=\mathcal{X}$. Note that $\hat{\mathcal{C}}_{\bf u}$ is a bounded linear operator from $\mathcal{X}$ to $\mathcal{X}_{-1}$, and in fact it is of rank $1$. For convenience we also introduce the notation

[TABLE]

Proposition 2.1**.**

*For every ${\bf u}\in\mathcal{Y}^{0}_{+}$, the operator $\mathcal{B}_{\bf u}$ is the part of the operator

$\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}+\hat{\mathcal{C}}_{\bf u}\right)$ in $\mathcal{X}$.*

Proof.

First we show that the domain of the part of $\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}+\hat{\mathcal{C}}_{\bf u}\right)$ in $\mathcal{X}$ coincides with the domain of $\mathcal{B}_{\bf u}$. That is, we need to show that

[TABLE]

First assume that $y\in D_{2}$. We define the following sequences of functions

[TABLE]

Note that $y_{n}(0)=\pi_{n}(0)=0,\,\forall\,n\in\mathbb{N}$, and $y_{n}\to y,\,\pi_{n}\to\pi$ in $\mathcal{X}$, moreover $\forall\,n\in\mathbb{N},\,\,y_{n},\pi_{n}\in D(\hat{\mathcal{B}}_{\bf u})$.

Since

[TABLE]

in $\mathcal{X}_{-1}$; we have that

[TABLE]

But we are going to show that in fact

[TABLE]

To this end note that for $x>\frac{1}{n}$ we have (below for simplicity we are suppressing the arguments in the functions $y$, $\mu$ and $\pi$)

[TABLE]

since $\pi^{\prime}_{n}(x)+\mu\pi_{n}(x)=0,\,\,\forall\,x>\frac{1}{n}$.

While for $x\leq\frac{1}{n}$ we have

[TABLE]

Therefore we obtain

[TABLE]

since $y$ and $y^{\prime}$ belong to $L^{1}$, $\Phi_{\bf u}(y_{n})\rightarrow\Phi_{\bf u}(y)=y(0)$ ($y_{n}$ tends to $y$ in $\mathcal{X}$), $y\left(\frac{1}{n}\right)\rightarrow y(0)$ (note that $y$ is a continuous function), and $\pi\left(\frac{1}{n}\right)\rightarrow 1$, as $n\to\infty$. Hence we conclude that $y\in D_{1}$.

Next assume that $y\in D_{1}$ and take an arbitrary sequence of smooth functions $y_{n}$, such that $y_{n}\in W^{1,1}(0,m)$ and $y_{n}(0)=0$ for every $n\in\mathbb{N}$, and such that $y_{n}\to y$ in $\mathcal{X}.$ Then we have

[TABLE]

with convergence in $\mathcal{X}_{-1}$. Using the definition of the $\mathcal{X}_{-1}$ norm it is shown that (2.28) is equivalent to

[TABLE]

i.e. $w\in W^{1,1}(0,m)$, and $w(0)=0.$

Since $y_{n}\xrightarrow[\mathcal{X}]{}y$ and $\Phi_{\bf u}(y_{n})\to\Phi_{\bf u}(y),$ we have $w=y-\Phi_{\bf u}(y)\pi$, which implies that $y\in W^{1,1}(0,m)$, and $y(0)=w(0)+\Phi_{\bf u}(y)\pi(0)=\Phi_{\bf u}(y)$; that is $y\in D_{2}.$

Now, for $y\in D_{1},$ we have

[TABLE]

where we used that $(y-y(0)\pi)\in D\left(\hat{\mathcal{B}}_{\bf u}\right)$, and that $y\in D_{2}=D_{1}.$ Hence the proof is completed.

∎

From the point of view of applying our steady state framework to model (2.5)-(2.6), the significance of recovering the generator $\mathcal{B}_{\bf u}$ as part of $\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}+\hat{\mathcal{C}}_{\bf u}\right)$ is that, since the semigroups $\mathcal{T}(t)$ and $\mathcal{T}_{-1}(t)$ are similar, that is

[TABLE]

holds; their spectra coincide, see [20, Ch.II] for more details.

This in particular means that we can study the behaviour of $s\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}+\hat{\mathcal{C}}_{\bf u}\right)$ as a function of ${\bf u}$, instead of studying directly the behaviour of $s(\mathcal{B}_{\bf u})$ as a function of the parameter ${\bf u}$, and vice versa.

Note that if we define the positive cone of $\mathcal{X}_{-1}$ as the completion of $\mathcal{X}_{+}$ with respect to the norm $||\cdot||_{-1}$, then it is shown that $\mathcal{X}_{-1}$, with the usual partial ordering, is a Banach lattice itself.

We now formulate conditions which guarantee the existence of a (strictly) positive steady state of the nonlinear model (2.5)-(2.6).

Theorem 2.2**.**

Assume that $\beta$ is a strictly monotone decreasing function with respect to its second variable, and that

1. (1)

$\mu>\mu_{0}$, for some $\mu_{0}>0$. 2. (2)

$\exists\,K>0\quad\text{such that}\quad\displaystyle\max_{a\in[0,m]}\,\{\beta(a,K)\}<\mu_{0}$. 3. (3)

$\beta(m,x)>0,\,\,\forall\,x\in\mathbb{R}_{+}$.

Then, if there exists an $r>0$, such that $s(\mathcal{B}_{\bf u})>0$ holds for all ${\bf u}\in\mathcal{Y}^{0}_{+}$ , $||{\bf u}||\leq r$, then model (2.5)-(2.6) admits a strictly positive steady state.

Proof.

Similarly as we have applied our framework previously e.g. in [6] for a model with infinite dimensional nonlinearity, but with distributed recruitment process; our goal is to transform the steady state problem into a fixed point problem in the parameter space. This requires us to establish some desirable properties of the level set $S$ (defined in (2.12)). In particular, we are going to establish the following properties.

• (i)

$s\left(\mathcal{B}_{\bf u}\right)$ is a continuous function of the parameter ${\bf u}$.

• (ii)

$s\left(\mathcal{B}_{\bf u}\right)$ is strictly monotone decreasing (as a function of ${\bf u}$) along rays in $\mathcal{Y}_{+}^{0}$.

• (iii)

The zero level set $S\subset\mathcal{Y}_{+}^{0}$ defined in (2.12) is contained in a bounded ball.

To establish the continuity of $s\left(\mathcal{B}_{\bf u}\right)$ as a function of the parameter ${\bf u}$, previously (see [7]) we used the notion of generalised convergence of operators and perturbation results from [28]. Here, in contrast, we employ a more direct approach, noting that $\mathcal{B}_{\bf u}$ generates an eventually compact positive semigroup, therefore its spectral bound $s\left(\mathcal{B}_{\bf u}\right)$ is an isolated eigenvalue, for every ${\bf u}\in\mathcal{Y}_{+}^{0}$. In fact $s\left(\mathcal{B}_{\bf u}\right)$ is the dominant real solution of the characteristic equation

[TABLE]

We now set

[TABLE]

and apply the Implicit Function Theorem. In particular we use Theorem A from Appendix A in [12]. Note that if we naturally extend the functions $\beta$ and $\mu$ as $\beta(\cdot,x)=\beta(\cdot,0),\,\forall x<0$, and $\mu(\cdot,x):=\mu(\cdot,0),\,\forall\,x<0$, then $G$ as defined above in (2.35) (on an open set) is a continuous map.

Also note that for any $\left({\bf u}_{*},s(\mathcal{B}_{{\bf u}_{*}})\right)\in\left\{\mathcal{Y}\setminus\left\{{\bf u}\in\mathcal{Y}:\int_{0}^{m}{\bf u}(r)\,\mathrm{d}r=0\right\}\right\}\times\mathbb{R}$ we have

[TABLE]

and the map $\lambda\to G({\bf u},\lambda)$ is continuously differentiable for every parameter value ${\bf u}\in\mathcal{Y}\setminus\left\{{\bf u}\in\mathcal{Y}:\int_{0}^{m}{\bf u}(r)\,\mathrm{d}r=0\right\}$.

In fact, since

[TABLE]

we have that the map $x\to\frac{\partial G}{\partial\lambda}\left({\bf u}_{*},s(\mathcal{B}_{{\bf u}_{*}})\right)x$ is a linear homeomorphism from $\mathbb{R}$ onto $\mathbb{R}$.

Then, the Implicit Function Theorem implies that $s(\mathcal{B}_{{\bf u}_{*}})$ is a continuous function of ${\bf u}_{*}$ for all ${\bf u}_{*}\in\mathcal{Y}\setminus\left\{{\bf u}\in\mathcal{Y}:\int_{0}^{m}{\bf u}(r)\,\mathrm{d}r=0\right\}$.

Next we show that the function ${\bf u}\to s(\mathcal{B}_{\bf u})$ is strictly monotone decreasing along positive rays in $\mathcal{Y}$. Using the extended generator, we have for every ${\bf u}\in\mathcal{Y}_{+}^{0}$, $v\in\mathcal{X}_{+}$, and $0<\alpha_{1}<\alpha_{2}$ (note that ${\bf u}$ is a parameter, while $v$ is an element of the positive cone of the state space $\mathcal{X}$)

[TABLE]

a positive perturbation.

Indeed, using that $F_{1}$ is invariant along rays and the notation (2.17), we have

[TABLE]

This is because $F_{2}(\alpha_{1}{\bf u})-F_{2}(\alpha_{2}{\bf u})=(\alpha_{1}-\alpha_{2})\int_{0}^{m}{\bf u}(a)da<0$, and hence the integral above is positive because $\beta$ decreases with respect to its second variable and, moreover, we have $\left(-\left(\hat{\mathcal{B}}_{{\bf u}}\right)_{-1}\pi\right)\in\left(\mathcal{X}_{-1}\right)_{+}$, because

[TABLE]

where the limits above are in $\mathcal{X}_{-1}$, (see also the proof of Prop. 2.1).

Note that $\left(\hat{\mathcal{B}}_{\alpha_{1}{\bf u}}\right)_{-1}-\left(\hat{\mathcal{B}}_{\alpha_{2}{\b{u}}}\right)_{-1}=0$ follows from the fact that the parts of $\left(\hat{\mathcal{B}}_{\alpha_{1}{\bf u}}\right)_{-1}$ and $\left(\hat{\mathcal{B}}_{\alpha_{2}{\bf u}}\right)_{-1}$ in $\mathcal{X}$ satisfy $\hat{\mathcal{B}}_{\alpha_{1}{\bf u}}=\hat{\mathcal{B}}_{\alpha_{2}{\bf u}}$, and the only unique continuous extension of the zero operator from a dense subset to $\mathcal{X}$ is the zero operator.

As we have noted before, $\mathcal{B}_{\bf u}$ generates a positive and eventually compact semigroup of operators, which is also irreducible under our hypotheses, and therefore the spectral bound $s(\mathcal{B}_{\bf u})$ belongs to its spectrum and it is an eigenvalue of algebraic multiplicity one; and as a consequence the same holds for $s\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}+\hat{\mathcal{C}}_{\bf u}\right)$.

Also note that for every sufficiently large $\lambda$, we have the identity

[TABLE]

and therefore Proposition A.2 in [1] (see also Theorem 1.3 in [2]) implies that

[TABLE]

which indeed shows that the function ${\bf u}\to s(\mathcal{B}_{\bf u})$ is strictly monotone decreasing along positive rays in $\mathcal{Y}$.

Finally we are going to prove that the zero level set $S$ is contained in a bounded ball. To this end we let

[TABLE]

We are going to show that for all ${\bf u}\in\mathcal{Y}_{+}^{0}$ such that $||{\bf u}||>K$ holds, we have $\lambda_{\bf u}<0$, which implies the assertion that $S$ is contained in a bounded ball.

Note that for every ${\bf u}\in\mathcal{Y}_{+}^{0}$, $\lambda_{\bf u}$ is a dominant real eigenvalue of multiplicity one with a corresponding strictly positive eigenvector $p_{\bf u}$, satisfying the equation

[TABLE]

which is equivalent to

[TABLE]

Note that the definition of $\hat{\mathcal{B}}_{\bf u}$ (in particular the homogeneous boundary condition) implies that

[TABLE]

which means that for any $\lambda$ the operator $\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}-\lambda\,\mathcal{I}\right)$ is invertible. Hence from (2.41) we obtain

[TABLE]

which yields (since $p_{\bf u}$ is strictly positive)

[TABLE]

Note that $F$ is monotone decreasing for $\lambda_{\bf u}>0$, and at $\lambda_{\bf u}=0$ we have

[TABLE]

for every ${\bf u}\in\mathcal{Y}_{+}^{0}$ satisfying $||{\bf u}||>K$, which shows that the spectral bound indeed changes sign in a bounded ball along every ray in $\mathcal{Y}^{0}_{+}$, and hence $S$ itself is contained in a bounded ball.

We also note that the fact that $S$ is bounded away from the origin follows directly from our last assumption in Theorem 2.2.

To complete the proof let us denote by $B_{+}$ the unit sphere of $\mathcal{Y}$ intersected with the positive cone $\mathcal{Y}_{+}$, and define a map $\Theta$ from the closed and convex set $B_{+}$ into itself as:

[TABLE]

Above, $h$, defined as $h({\bf u})=\frac{{\bf u}}{||{\bf u}||}$, is continuous with a continuous inverse (this is because $S$ is bounded away from the origin and it is contained in a bounded ball). The projection $p$ along positive rays in $L_{+}^{1}$ is also continuous and bounded. $k$ is the compact injection of $W^{1,1}(0,m)$ into $L^{1}(0,m)$. The map $e$, which assigns the strictly positive normalised eigenvector ${\bf V}_{\bf u}$ of $\mathcal{B}_{\bf u}$ to $\mathcal{B}_{\bf u}$, is analytic, (see e.g. [11, Lemma 1.3]). Finally, the map $g$ above, which assigns the linear operator $(\mathcal{B}_{\bf u},D(\mathcal{B}_{\bf u}))\in C(\mathcal{X})$ (the set of closed operators on $\mathcal{X}$) to the parameter value ${\bf u}\in S$, is shown to be continuous using the notion of generalised convergence, see [28, Sect.2 Ch. IV]. (Note that compared to e.g. [6, 8] the situation is more delicate, as both the operator and its domain depend on the parameter value ${\bf u}$.) In particular, continuity is established via [28, Th.2.25 Ch. IV.], which characterizes generalised convergence of operators via convergence of their resolvents. Since the resolvent of $\mathcal{B}_{\bf u}$ is explicitly given as ($f\in L^{1}(0,m)$)

[TABLE]

it is shown, by applying Lebesgue’s dominated convergence theorem, that for any sequence $\left\{{\bf u_{n}}\right\}\subset S$ tending to ${\bf\hat{u}}\in S$ there exists a $\lambda$ large enough, such that

[TABLE]

as $n\to\infty$.

By Schauder’s fixed point theorem the map $\Theta$ has a fixed point. Equivalently, the map $h^{-1}\circ\Theta\circ h$, which maps $S$ into itself, has a fixed point ${\bf u}_{*}$ which turns out to be a steady state of the problem:

[TABLE]

for some $0\neq\alpha({\bf V}_{{\bf u}_{*}})\in\mathbb{R}$; since ${\bf V}_{{\bf u}_{*}}$ belongs to the kernel of $\mathcal{B}_{{\bf u}_{*}}$. ∎

The following (counterpart) result can be established along similar lines.

Theorem 2.3**.**

Assume that $\beta$ is a strictly monotone increasing function with respect to its second variable, and that

1. (1)

$\beta(m,x)>0,\,\,\forall\,x\in\mathbb{R}_{+}$. 2. (2)

There exists an $r>0$, such that $s(\mathcal{B}_{\bf u})<0$ holds for all ${\bf u}\in\mathcal{Y}_{+}^{0}$ , $||{\bf u}||\leq r$. 3. (3)

There exists an $R>r>0$, such that $s(\mathcal{B}_{\bf u})>0$ holds for all ${\bf u}\in\mathcal{Y}_{+}^{0}$ , $||{\bf u}||\geq R$.

Then, model (2.5)-(2.6) admits a strictly positive steady state.

As we have already noted in [7], the crucial condition that the spectral bound is monotone (and changes sign) along positive rays of the parameter space is closely related to the biologically meaningful condition that an appropriately defined net reproduction function, or in our case a net reproduction functional, changes monotonously (and crosses the value $1$) along positive rays of the parameter space. We now further explore this connection, and in particular we recover the familiar net reproduction number/function for age-structured (single state at birth) models.

To this end, for every fixed ${\bf u}\in\mathcal{Y}_{+}^{0}$ we define the net reproduction number as

[TABLE]

(notice that $\mathcal{R}({\bf u})$ coincides with $F(0)$ as defined in (2.44)), and for ${\bf u}\equiv{\bf 0}$ (i.e. the extinction state/environment) we may naturally define the net reproduction number as

[TABLE]

Note that naturally $\mathcal{R}\,:\,\mathcal{Y}_{+}^{0}\to\mathbb{R}_{+}$ can be considered as a functional, and so we may call it the net reproduction functional, see also [21] for more details in this direction. Also note that for the specific model we considered here we had to formally define $\mathcal{R}$ at the extinction steady state (zero population density) separately, since $F_{1}$ is not defined at ${\bf 0}$ in general. Note that in general the value $\mathcal{R}({\bf 0})$ is often denoted by simply $\mathcal{R}$ or $\mathcal{R}_{0}$ is usually referred to as the net reproduction number or basic reproduction number in the literature. We would like to emphasize though that for nonlinear models one really has to keep in mind that these scalar values (for example $\mathcal{R}_{0}$) are simply the values of a net reproduction function/functional at fixed population densities (environments). Indeed, net reproduction numbers, functions and functionals play an important role when studying the existence and the local asymptotic stability of equilibria of physiologically structured population dynamical models, for more details see e.g. the recent papers [4, 5, 15, 19, 21, 22, 23, 35, 36].

We now rigorously establish the connection between the spectral bound $s(\mathcal{B}_{\bf u})$ and the net reproduction functional $\mathcal{R}({\bf u})$, as defined in (2.49).

Proposition 2.4**.**

For every ${\bf u}\in\mathcal{Y}_{+}^{0}$ we have

[TABLE]

and in particular

[TABLE]

Proof.

To establish the relationship (2.51) for ${\bf u}\in\mathcal{Y}_{+}^{0}$ we apply Theorem 3.5 from [36] to the extended generator $\left(\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}+\hat{\mathcal{C}}_{\bf u}\right)$ in $\mathcal{X}_{-1}$. Note that one can verify that the operators $\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1},\,\hat{\mathcal{C}}_{\bf u}$ satisfy the assumptions of Theorem 3.5 in [36]. Hence applying Theorem 3.5 from [36] we obtain

[TABLE]

It remains to show that

[TABLE]

To this end note that the equation

[TABLE]

shows that $\lambda=0$ is an eigenvalue with corresponding eigenvectors $p=\left(\hat{\mathcal{B}}_{\bf u}\right)_{-1}w$ (spanning an infinite dimensional subspace), such that $w\in\mathcal{X}$ and

[TABLE]

holds.

On the other hand if $\lambda\neq 0$ then we have that

[TABLE]

for some $\alpha\neq 0$, which together with

[TABLE]

yields

[TABLE]

Hence we have shown that (below $\sigma$ stands for the spectrum)

[TABLE]

which establishes the assertions of Proposition 2.4. ∎

**Remark 2.5  ** Note that by appropriately defining the operator $\mathcal{B}_{\bf 0}$ we can also establish that $s(\mathcal{B}_{\bf 0})\gtrless 0\iff R({\bf 0})\gtrless 1$.

**Remark 2.6  ** We would like to emphasize that in the proof of Proposition 2.4 we did not really use that $F_{1},F_{2}$ assume the particular form as in (2.7).

**Remark 2.7  ** We recall from [26] the classic nonlinear (Gurtin-McCamy) age-structured population dynamical model here.

[TABLE]

Note that the proof of Proposition 2.4 can be directly adapted for this model (with a one dimensional nonlinearity, i.e. $\mathcal{Y}=\mathbb{R}$, since the environment individuals are experiencing is simply determined by the total population size $P(t)=\int_{0}^{m}p(a,t)\,\mathrm{d}a$), and of course also for the basic linear age-structured McKendrick model (i.e. when $\beta$ and $\mu$ do not depend on the total population size $P(t)$).

Thus we have rigorously justified the biologically inspired definition of the net reproduction number

[TABLE]

(appearing frequently in the literature, see e.g. [27]) for the linear age-structured McKendrick model; and similarly the net reproduction function, which, for the Gurtin-McCamy (non-linear) age-structured model (2.55)-(2.56) reads as:

[TABLE]

Also note that analogously one arrives at the definition of the net reproduction function

[TABLE]

for the size-structured version of the Gurtin-McCamy model (2.55)-(2.56), see e.g. [23], where stability criteria for positive steady states were also established in terms of the derivative of the net reproduction function above.

3. Concluding remarks

In this paper we have treated the (positive) steady state problem of a nonlinear structured population model with a nonlinear and non-local boundary condition. In particular we have showed that the growth behaviour of the spectral bound of the parametrised family of linear operators arising in the steady state formulation can be controlled by boundary perturbations only. This way we have demonstrated how our general framework developed in [7] can be applied to treat classes of nonlinear models with single state at birth, and in particular with fertility/reproduction controlled dynamics.

We note that for certain classes of models, a more direct fixed point approach as employed previously e.g. in [24] can be also applied. The question is whether one formulates the fixed point map on the state space as e.g. in [24], or on the parameter (environment) space. Sometimes, in particular when it is not feasible to obtain an implicit formula for the steady state, it may not be possible to formulate a fixed point map on the state space. Another disadvantage of formulating the fixed point problem on the state space is that often fixed point theorems rely on contraction arguments, and so in the context of population models special care has to be taken to avoid the extinction steady state, as e.g. in [24], where a fixed point theorem on conical shells in Banach spaces was employed.

For these reasons in general we prefer to formulate a fixed point problem on the parameter space, and make use of the spectral properties of the family of operators arising from the parametrisation (using the environment as parameter), partly because then this has a natural biological intepretation, too. In fact, the goal of the present work was to show how the spectral bound of the family of operators can be controlled by boundary perturbations only, and the results here can be readily extended to nonlinear size-structured models with non-local boundary condition, see e.g. [7] or even [8] (including a diffusion operator, too) in this direction.

The general idea of lifting a problem into the extrapolated space, where the boundary condition becomes a bounded perturbation of the generator is very powerful and can be applied to other models, as well. To name another concrete example, let us consider here the following general selection-mutation model (see also [9] for more details).

[TABLE]

In the model above $u(l,a,t)$ denotes the density of individuals of age $a$ and maturation age $l$ at time $t$. $E_{1}$ and $E_{2}$ are the environmental operators, as described in Section 1. $\mu$ denotes the mortality rate, $\beta$ the fertility function and $b$ is the probability density function describing the mutation, that is $\int_{l_{1}}^{l_{2}}b(l,\hat{l})\,\mathrm{d}l$ is the probability that the offspring of an individual of age at maturity $\hat{l}$ has maturity age $l\in(l_{1},l_{2})$.

In [9] the existence and uniqueness of a positive steady state of model (3.60) was established, for the case when $E_{1}\equiv E_{2}$, and when the mortality $\mu$ is a strictly monotone increasing function of its first variable, while the fertility $\beta$ is a strictly monotone decreasing function of its first variable, as well. In contrast, in [7] we treated the case of $E_{1}\equiv E_{2}\equiv E[u(\cdot,\cdot,t)]=(P(t),Q(t))$ (i.e. two-dimensional nonlinearity), where

[TABLE]

(note the finite maximal age) denote the juvenile and adult populations respectively, but crucially without any further monotonicity assumptions on $\mu$ and $\beta$.

Here we point out (without elaborating the details) that following the lines presented in the previous section, the positive steady state problem for model (3.60) could be treated, for example in case of the following choice of $E_{1}$ and $E_{2}$:

[TABLE]

where $w$ is some non-negative weight function and $\alpha\in[0,1]$.

Finally we mention that in the future we intend to explore the possibility to use our framework in the dynamic setting, that is, we aim to investigate the possibility of using the (re)formulation (1.1) not just to prove the existence of positive steady states, but also to establish further qualitative results, for example results concerning the global (in)stability of the extinction steady state. We anticipate that lifting the problem into the extrapolated space will prove to be fruitful again. We also note that it will be worthwhile to explore the connection between our general approach and the theory of monotone dynamical systems and in particular persistence theory (see e.g. [34]), and to compare the advantages/disadvantages of the two approaches.

Acknowledgments

The authors were partially supported by the research project MTM2014-52402-C03-2P.