A Nonlinear Population Model

1. Introduction

Our paper studies a population model with age structure. This model and its variants were considered in many works such as [2,3,6,10,11,14,16,17,20], to cite just a few. There is a vast literature on this subject by now.

The seminal work of [14] considers the age structure into the dynamics of one sex population model assuming that the female population dynamics can be modelled as a function of two variables, namely age and time. This model takes as inputs an age specific mortality intensity and an age specific fertility function. The number of individuals at a given time that have age less than a certain age is given by an integral of a function $p(a,t)$ of two variables (age and time). The size of population (the total number of individuals) can also be obtained by integrating $p(a,t)$. It appears only natural to find $p(a,t)$ which is described by a system of integral and differential equations. This system is shown to reduce to a Volterra integral equation. The later has a unique solution established by a fixed point theorem, but this solution is not available in closed form, and its numerical computation based on fixed point iterations seems complicated. The next class of population models are the nonlinear models obtained in situations when the mortality and fertility are functions of age and population size (see [12] for more on this).

Let us mention two seminal works in the paradigm of nonlinear population models. The first work to incorporate population size dependent mortality and fertility rates (rendering the model nonlinear) is [10]. They also characterize the equilibrium population and established a condition for the equilibrium to be locally asymptotically stable. A special class of models, the separable ones, are considered in their paper. [20] shows that the stability classification, depends in many cases on the marginal birth and death rates as measures the sensitivities of the fertility and mortality; on other cases more information is required to determine the stability.

Let us turn now to the contributions of our paper. Our setting is the nonlinear population model, in the special case when mortality and fertility are separable functions of age and population size. The age part of the fertility is assumed as in [18] which is a model with non constant (in age) fertility rate. In this paradigm we establish the existence of the steady state for the nonlinear equations characterizing the population dynamics. Moreover, these are shown to be equivalent to a nonlinear system of differential equations. The late is shown to have an equilibrium solution which we find explicitly. The equilibrium stability analysis can also be established.

The reminder of this paper is organized as follows: Section 2 presents the populations models, linear, nonlinear. Section 3 presents our nonlinear model and the main results of the paper. Section 4 contain some results that give new directions of study. The work end with some remarks in Section 5.

2. Population Models

Let us first introduce the linear population model. The exposition here follows [12]. The dynamics of population is expressed in terms of the density of the population of age a, at time t density denoted $p(a,t)$. The total population at t, denoted by $P\left(t\right)$ can be obtained by integrating its density, i.e.,

$$P\left(t\right)={\int }_0^{\infty }p(a,t)da.$$

Next let us introduce fertility and mortality. The number of offspring, borne by individuals during the infinitesimal time-interval $[t,t+dt]$, and the infinitesimal age-interval $[a,a+da]$ is

$$\beta (a,t),$$

also referred to as the age specific fertility. The number of offspring during the infinitesimal interval $[t,t+dt]$ is then

$$B\left(t\right)={\int }_0^{\infty }\beta (a,t)p(a,t)da.$$

The number of deaths of individuals during the infinitesimal time-interval $[t,t+dt]$, and the infinitesimal age-interval $[a,a+da]$ is

$$\mu (a,t).$$

The number of deaths during the infinitesimal interval $[t,t+dt]$ is then

$$D\left(t\right)={\int }_0^{\infty }\mu (a,t)p(a,t)da.$$

The probability that an individual of age $a-x$ at the time $t-x$ will survive up to time t (with age a) is given by

$$\pi (a,t,x)=e^{-{\int }_0^x\mu (a-\sigma ,t-\sigma )d\sigma }.$$

In the case of time independent mortality

$$\pi \left(a\right)=e^{-{\int }_0^a\mu \left(\sigma \right)d\sigma },$$

is the probability for a newborn to survive to age a, also known as the survival probability.

In the following we will derive the linear Lotka-McKendrick Equation. The fertility and mortality rates $\beta \left(a\right),$ and $\mu \left(a\right)$ are assume time independent, as they only depend on the age $a.$ The number of individuals with age less than a at time $t,$ denoted by $N(a,t)$ is given by

$$N(a,t)={\int }_0^{a}p(\sigma ,t)d\sigma .$$

Next, let us look at the number of individuals with age less than $a+h$ at time $t+h,$ i.e., $N(a+h,t+h)$. This number will comprise $N(a,t),$ and all newborn in the time interval $[t,t+h]$ (their age will be less than $a+h$), which is

$${\int }_t^{t+h}B\left(s\right)ds.$$

One needs to adjust then for the deaths from the newborns, through the time interval $[t,t+h],$ and the deaths on $[t,t+h]$ of individuals older than $a,$ and these number of deaths is

$${\int }_0^h{\int }_0^{a+s}\mu \left(\sigma \right)p(\sigma ,t+s)d\sigma ds.$$

This is the case because

$${\int }_0^{a+s}\mu \left(\sigma \right)p(\sigma ,t+s)d\sigma ,$$

gives the number of individuals who die at $t+s$ younger than $a+s$. Therefore

$$N(a+h,t+h)=N(a,t)+{\int }_t^{t+h}B\left(s\right)ds-{\int }_0^h{\int }_0^{a+s}\mu \left(\sigma \right)p(\sigma ,t+s)d\sigma ds.$$

Let us differentiate this with respect to h, and then take $h=0$ to get

$$p(a,t)+{\int }_0^a{p}_t(\sigma ,t)d\sigma =B\left(t\right)-{\int }_0^a\mu \left(\sigma \right)p(\sigma ,t)d\sigma .$$

Next, let us differentiate this with respect to $a,$ to get

$$p(a,t)+p_a(a,t)+\mu \left(a\right)p(a,t)=0.$$

Also by setting $a=0$ yields

$$p(0,t)=B\left(t\right),$$

but on the other hand

$$B\left(t\right)={\int }_0^{\infty }\beta \left(\sigma \right)p(\sigma ,t)d\sigma .$$

As such we obtained the following system

$$\left\{\begin{array}{ccc}{p}_a(a,t)+p_t(a,t)+\mu \left(a\right)p(a,t)=0\hfill & \mathrm{for}\hfill & a,t\ge 0,\hfill \\ p(0,t)=\underset{0}{\overset{\infty }{\int }}\beta \left(\sigma \right)\rho (\sigma ,t)d\sigma \hfill & \mathrm{for}\hfill & t>0,\hfill \\ p(a,0)=p_0\left(a\right)\hfill & \mathrm{for}\hfill & a>0.\hfill \end{array}\right.$$

(1)

2.1. A Special Linear Population Model

Inspired by [1,7] and [8], the recent work [4] considers a model with survival probability $\pi \left(a\right)$ pseudo exponential, i.e.

$$\pi \left(a\right)=\sum _{i=1}^n{c}_i{e}^{-{\mu }_{i}a},$$

(2)

for positive constants $c_i,{\mu }_i$ with

$$\sum _{i=1}^n{c}_i^2\ne 0,\sum _{i=1}^n{\mu }_i^2\ne 0\mathrm{and}\sum _{i=1}^n{c}_i=1.$$

Moreover, we present the case of

$$\pi \left(a\right)=\left(\sum _{i=0}^n{c}_i{a}^i\right)e^{-{\mu }_{1}a}.$$

(3)

In this setting finding $p(a,t)$ is reduced to a linear ODE system. In special situations (2)-(3) a closed form solution is obtained by means of Laplace transform.

2.2. Nonlinear Population Models

These are models in which death rate $\mu (a,p),$ and the fertility rate $\beta (a,p),$ are functions of age and population size. The age density function at time $t\ge 0,$$p\left(a,t\right)$ satisfies the following nonlinear equations

$$\left\{\begin{array}{c}{p}_t\left(a,t\right)+p_a\left(a,t\right)+\mu (a,P\left(t\right))p\left(a,t\right)=0\hfill \\ p\left(0,t\right)={\int }_0^{\infty }\beta (\sigma ,P\left(t\right))p\left(\sigma ,t\right)d\sigma \hfill \\ p\left(a,0\right)=p_0\left(a\right)\hfill \\ P\left(t\right)={\int }_0^{\infty }p\left(\sigma ,t\right)d\sigma ,\hfill \end{array}\right.$$

where $P\left(t\right)$ is the total population at t. This system of nonlinear equations can be reduced to integral equations, as it is shown in (see [12] for more details). Indeed, let

$$\pi (a,t,x,P)=e^{-{\int }_0^x\mu (a-\sigma ,P(t-\sigma ))d\sigma }$$

and

$$p(a,t)=\left\{\begin{array}{ccc}{p}_0(t-a)\pi (a,t,t,P)\hfill & for& \hfill t\le a\\ B(t-a)\pi (a,t,a,P)\hfill & for& \hfill t>a.\end{array}\right.$$

(4)

Here the function B solves the following system of equations:

$$\left\{\begin{array}{ccc}B\left(t\right)={\int }_0^t\beta (\sigma ,P\left(t\right))\pi (\sigma ,t,\sigma ,P)B\left(\sigma \right)d\sigma +F(t,P)\hfill & & \\ P\left(t\right)={\int }_0^t\pi (\sigma ,t,\sigma ,P)B\left(\sigma \right)d\sigma +G(t,P)\hfill & \end{array}\right.$$

(5)

where

$$\left\{\begin{array}{ccc}F(t,P)={\int }_t^{\infty }\beta (\sigma ,P\left(t\right))\pi (\sigma ,t,t,P)p_0(\sigma -t)d\sigma \hfill & & \\ G(t,P)={\int }_t^{\infty }\pi (\sigma ,t,t,P)p_0(\sigma -t)d\sigma .\hfill & \end{array}\right.$$

(6)

This system (5) can be solved through the following iterative method

$$\left\{\begin{array}{ccc}{B}^{k+1}\left(t\right)={\int }_0^t\beta (\sigma ,P^k\left(t\right))\pi (\sigma ,t,\sigma ,P^k)B^k\left(\sigma \right)d\sigma +F(t,P^k)\hfill & & \\ P^{k+1}\left(t\right)={\int }_0^t\pi (\sigma ,t,\sigma ,P^k)B^k\left(\sigma \right)d\sigma +G(t,P^k),\hfill & \end{array}\right.$$

(7)

see [12] for more details on this. A special class of nonlinear models are the separable population models presented in the next subsection.

2.3. Separable Population Models

Now let us specialize the nonlinear model with the following choice of fertility and mortality

$$\beta (a,p)=R_0{\beta }_0\Phi \left(p\right)e^{-\rho a},\mu (a,p)={\mu }_0+\Psi \left(p\right).$$

By plugging this in the nonlinear system one gets

$$\left\{\begin{array}{c}{p}_t\left(a,t\right)+p_a\left(a,t\right)+({\mu }_0+\Psi \left(P\left(t\right)\right))\left(a,t\right)=0\hfill \\ p\left(0,t\right)={\int }_0^{\infty }{R}_0\Phi \left(P\left(t\right)\right)e^{-\rho a}p\left(\sigma ,t\right)d\sigma \hfill \\ p\left(a,0\right)=p_0\left(a\right)\hfill \\ P\left(t\right)={\int }_0^{\infty }p\left(\sigma ,t\right)d\sigma .\hfill \end{array}\right.$$

These nonlinear equations can be reduced to the following ODE system

$$\left\{\begin{array}{c}{P}^{\prime }\left(t\right)=-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P\left(t\right)+R_0{\beta }_0\Phi \left(P\left(t\right)\right)Q\left(t\right)\hfill \\ Q^{\prime }\left(t\right)=\left(R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\right)Q\left(t\right)\hfill \end{array}\right.$$

(8)

where

$$Q\left(t\right)={\int }_0^{\infty }{e}^{-a\rho }{p}_t\left(a,t\right)da,$$

see [12] for more details.

Let us turn now to the existence of steady states. The net reproduction number is given by

$$R\left(x\right)=R_0\Phi \left(x\right){\int }_0^{\infty }{\beta }_0{e}^{-\left(\rho +{\mu }_0+\Psi \left(x\right)\right)a}da.$$

This quantity was introduced first by [10]. According with [10] the quantity $R\left(x\right)$ is the number of children expected to be born to an individual when the population is x. The steady state $P^{*}$ is given by the following equation

$$R\left(P^{*}\right)=1,$$

which in this setting reads

$$R_0{\beta }_0\Phi \left(P^{*}\right)=\rho +{\mu }_0+\Psi \left(P^{*}\right).$$

3. A Special Nonlinear Population Model

We consider the following logistic system

$$\left\{\begin{array}{c}{p}_t\left(a,t\right)+p_a\left(a,t\right)+{\mu }_{0}p\left(a,t\right)+\Psi \left(P\left(t\right)\right)p\left(a,t\right)=0\hfill \\ p\left(0,t\right)=R_0\Phi \left(P\left(t\right)\right){\int }_0^{\infty }\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{\sigma }^i{e}^{-\rho \sigma }p\left(\sigma ,t\right)d\sigma \hfill \\ p\left(a,0\right)=p_0\left(a\right)\hfill \\ P\left(t\right)={\int }_0^{\infty }p\left(\sigma ,t\right)d\sigma \hfill \end{array}\right.$$

(9)

where $p\left(a,t\right)$ is the age density function at time $t\ge 0$, $P\left(t\right)$ is the total population at time t, $p_0\left(a\right)$ is the initial population at time $t=0$, $a\in \left[0,\infty \right)$, ${\mu }_0>0$ is the intrinsic mortality term, $R_0>0$ is the birth modulus and $\rho$, $\sigma$ are prescribed positive parameters. In real world, the age profile is given by

$$\frac{p\left(a,t\right)}{P\left(t\right)}.$$

The fertility is given by

$$\beta \left(a,p\right)=R_0\Phi \left(p\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{a}^i{e}^{-\rho a}\mathrm{with}{\beta }_0,{\beta }_1,...,{\beta }_{n-1}\in \left(0,\infty \right),$$

and the mortality (death rate) by

$$\mu (a,P)={\mu }_0+\Psi \left(P\right).$$

The couple $\left(\beta \left(a\right),{\mu }_0\right)$ can be interpreted as an intrinsic birth-death process that is age dependent while $\Psi \left(P\left(t\right)\right)$ models an external mortality that is the same for all ages and just depends on the weighted sizes.

Since we are interested in the existence of steady states solutions for the system (9), i.e. for solutions that are constant in time, we assume for the start that $\Phi$ and $\Psi$ are continuous on $\left[0,\infty \right)$ continuously differentiable on $\left(0,\infty \right)$ and

$$\begin{array}{ccc}\hfill \Phi \left(x\right)& \ge & 0,{\Phi }^{\prime }\left(x\right)<0,\Phi \left(0\right)=1,\Phi \left(+\infty \right)=0,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \Psi \left(x\right)& \ge & 0,{\Psi }^{\prime }\left(x\right)>0,\Psi \left(0\right)=0,\Psi \left(+\infty \right)=+\infty .\hfill \end{array}$$

(11)

We also adopt the normalization condition

$${\int }_0^{\infty }\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{a}^i{e}^{-\left(\rho +{\mu }_0\right)a}da=1,$$

from where, for example with the use of Gamma integral, we obtain

$$\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0\right)}^{i+1}}=1,$$

(12)

so that the parameter $R_0$ in (9) has the role of an intrinsic basic reproduction number denoted by $R\left(x\right)$ in the next.

Concerning the existence of steady states, the net reproduction number at size x (see page 154 in [12] or page 288 in [10]) takes the following form within our setting

$$\begin{array}{ccc}\hfill R\left(x\right)& =& R_0\Phi \left(x\right){\int }_0^{\infty }\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{a}^i{e}^{-\left(\rho +{\mu }_0+\Psi \left(x\right)\right)a}da\hfill \\ & =& R_0\Phi \left(x\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{\int }_0^{\infty }{a}^i{e}^{-\left(\rho +{\mu }_0+\Psi \left(x\right)\right)a}da\hfill \\ & =& R_0\Phi \left(x\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{\int }_0^{\infty }{\left(\frac{t}{\rho +{\mu }_0+\Psi \left(x\right)}\right)}^i{e}^{-t}\frac{dt}{\rho +{\mu }_0+\Psi \left(x\right)}\hfill \\ & =& R_0\Phi \left(x\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_i}{{\left(\rho +{\mu }_0+\Psi \left(x\right)\right)}^{i+1}}{\int }_0^{\infty }{t}^i{e}^{-t}dt\hfill \\ & =& R_0\Phi \left(x\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_i\Gamma \left(i+1\right)}{{\left(\rho +{\mu }_0+\Psi \left(x\right)\right)}^{i+1}}\hfill \\ & =& R_0\Phi \left(x\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(x\right)\right)}^{i+1}}\hfill \end{array}$$

where $\Gamma$ represent the Gamma integral. A first observation regarding the set of assumptions (10)-(11) is that

$$\underset{x\to \infty }{lim}R\left(x\right)=0\mathrm{and}R\left(x\right)\mathrm{is}\mathrm{a}\mathrm{decreasing}\mathrm{function}(\mathrm{i}.\mathrm{e}.R^{\prime }\left(x\right)<0).$$

(13)

Indeed $R\left(x\right)$ is a decreasing function because

$$R^{\prime }\left(x\right)=\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_{i}i!\left[R_0{\Phi }^{\prime }\left(x\right)\left(\rho +{\mu }_0+\Psi \left(x\right)\right)-R_0\Phi \left(x\right){\Psi }^{\prime }\left(x\right)\right]}{{\left(\rho +{\mu }_0+\Psi \left(x\right)\right)}^{i+2}}-\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{i{\beta }_{i}i!R_0\Phi \left(x\right){\Psi }^{\prime }\left(x\right)}{{\left(\rho +{\mu }_0+\Psi \left(x\right)\right)}^{i+2}}<0,$$

for all $x\ge 0$. Moreover,

$$\underset{x\to \infty }{lim}R\left(x\right)=0,$$

is satisfied in light of the asymptotic conditions on $\Psi ,\Phi .$

As is well known (see page 154 in [12] or [10, Theorem 6, pages 288-289]), since we have a single weighted size P, we use that non-trivial stationary sizes $P^{*}$ must satisfy

$$R\left(P^{*}\right)=1,$$

(14)

and this is a necessary and sufficient condition for a non-trivial stationary sizes to exist with total population $P^{*}$. In this case, (14) becomes

$$\frac{R_0\Phi \left(P^{*}\right)}{\rho +{\mu }_0+\Psi \left(P^{*}\right)}\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}=1.$$

(15)

Due to (13) this equation (15) has one, and only one, non-trivial solution if, and only if, $R_0>1$. The fact that the condition $R_0>1$ is necessary and sufficient for the existence of a non-trivial equilibrium means that $R_0$ acts as a bifurcation parameter: under the assumptions of the model it is clear that we have a forward bifurcation at the point $R_0$.

One can rewrite (15) for the non-trivial stationary sizes $P^{*}$ in the form

$$\stackrel{n-1}{\underset{i=0}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}=\frac{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{R_0\Phi \left(P^{*}\right)}.$$

(16)

Let us summarize the results in the following Lemma.

Lemma 1. 

Given assumptions (10), (11) our model has a unique steady state $P^{*}$ given by equation (15) if and only if $R_0>1.$

Now let us turn to the problem of finding $P\left(t\right).$ We will denote by

$$P_i\left(t\right)={\int }_0^{\infty }{\sigma }^{i-1}{e}^{-\rho \sigma }p\left(\sigma ,t\right)d\sigma ,$$

for $i=1,2,..n$. The next step, is to observe that the renewal condition or the total birth rate or fertility rate, at the time t can be can be written in the new notations such

$$\begin{array}{ccc}\hfill p\left(0,t\right)& =& R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{\int }_0^{\infty }{\sigma }^i{e}^{-\rho \sigma }p\left(\sigma ,t\right)d\sigma \hfill \\ & =& R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}\left(t\right).\hfill \end{array}$$

More that, the calculation of the first derivative of $P\left(t\right)$ gives

$$\begin{array}{ccc}\hfill P^{\prime }\left(t\right)& =& {\int }_0^{\infty }{p}_t\left(a,t\right)da\hfill \\ & =& -{\int }_0^{\infty }{p}_a\left(a,t\right)da-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right){\int }_0^{\infty }p\left(a,t\right)da\hfill \\ & =& p\left(0,t\right)-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P\left(t\right)\hfill \\ & =& -\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P\left(t\right)+R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}\left(t\right)\hfill \end{array}$$

and, similarly

$$\begin{array}{ccc}\hfill P_1^{\prime }\left(t\right)& =& {\int }_0^{\infty }{e}^{-a\rho }{p}_t\left(a,t\right)da\hfill \\ & =& -{\int }_0^{\infty }{e}^{-a\rho }{p}_a\left(a,t\right)da-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right){\int }_0^{\infty }{e}^{-a\rho }p\left(a,t\right)da\hfill \\ & =& p\left(0,t\right)-\rho {\int }_0^{\infty }{e}^{-a\rho }{p}_a\left(a,t\right)da-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P_1\left(t\right)\hfill \\ & =& R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}-\rho P_1\left(t\right)-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P_1\left(t\right)\hfill \\ & =& \left(R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\right)P_1\left(t\right)+R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=1}{\Sigma }}{\beta }_i{P}_{i+1}.\hfill \end{array}$$

In the same way, the calculation of the first derivative of $P_{i+1}\left(t\right)$ ($i=1,...,n-1$) gives

$$\begin{array}{ccc}\hfill P_{i+1}^{\prime }\left(t\right)& =& {\int }_0^{\infty }{a}^i{e}^{-\rho a}{p}_t\left(\sigma ,t\right)d\sigma \hfill \\ & =& -{\int }_0^{\infty }{a}^i{e}^{-a\rho }{p}_a\left(a,t\right)da-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right){\int }_0^{\infty }{a}^i{e}^{-a\rho }p\left(a,t\right)da\hfill \\ & =& {\int }_0^{\infty }\left(i{a}^{i-1}{e}^{-a\rho }-a^i\rho e^{-a\rho }\right)p\left(a,t\right)da-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right){\int }_0^{\infty }{a}^i{e}^{-a\rho }p\left(a,t\right)da\hfill \\ & =& i{\int }_0^{\infty }{a}^{i-1}{e}^{-a\rho }p\left(a,t\right)da-{\int }_0^{\infty }{a}^i\rho e^{-a\rho }p\left(a,t\right)da-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right){\int }_0^{\infty }{a}^i{e}^{-a\rho }p\left(a,t\right)da\hfill \\ & =& i{P}_i\left(t\right)-\rho P_{i+1}\left(t\right)-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P_{i+1}\left(t\right)\hfill \\ & =& i{P}_i\left(t\right)-\left(\rho +{\mu }_0+\Psi \left(P\left(t\right)\right)\right)P_{i+1}\left(t\right).\hfill \end{array}$$

Finally, to achieve our goal of obtaining the existence of solutions to the model (9) we are led to the system of differential equations of first order

$$\left\{\begin{array}{c}{P}^{\prime }\left(t\right)=-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P\left(t\right)+R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}\hfill \\ P_1^{\prime }\left(t\right)=\left(R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\right)P_1\left(t\right)+R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=1}{\Sigma }}{\beta }_i{P}_{i+1}\hfill \\ P_{i+1}^{\prime }\left(t\right)=i{P}_i\left(t\right)-\left(\rho +{\mu }_0+\Psi \left(P\left(t\right)\right)\right)P_{i+1}\left(t\right),i=1,...,n-1\mathrm{and}n\ge 2\hfill \end{array}\right.$$

(17)

coupled with the initial conditions

$$P\left(0\right)=P_0,P_i\left(0\right)=P_0^i,i=1,2,...,n$$

(18)

where

$$P\left(0\right)={\int }_0^{\infty }{p}_0\left(a\right)da,P_i\left(0\right)={\int }_0^{\infty }{a}^{i-1}{e}^{-\rho a}{p}_0\left(a\right)da,i=1,2,..n.$$

The study of the existence of solutions for the system (17) is equivalent to the study of existence of solutions to (9) because, if the pair

$$(P\left(t\right),P_1\left(t\right),...,P_n\left(t\right))$$

solves (17), then by setting

$$B\left(t\right)=p\left(0,t\right)=R_0\Phi \left(P\left(t\right)\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}\left(t\right)$$

we obtain the solution to (9) via the usual formula

$$p\left(a,t\right)=\left\{\begin{array}{cc}{p}_0\left(a-t\right)e^{-{\int }_0^t\left({\mu }_0+\Psi \left(P\left(\sigma \right)\right)\right)d\sigma }& a\ge t,\\ B\left(t-a\right)e^{-{\int }_{t-a}^t\left({\mu }_0+\Psi \left(P\left(\sigma \right)\right)\right)d\sigma }& a<t.\end{array}\right.$$

Let us summarize the results here.

Theorem 1. 

The system (9) is equivalent to the ordinary nonlinear system (17).

Thus, in order to determine the existence of solutions to the model (9) we can focus on the analysis of (17).

Clearly, the system (17) has at least the trivial solution and so the existence of stationary solutions to this problem (17) are of our interest that may not be unique for some values of the parameters and then may lead to complex bifurcations. In the model of the problem any stationary solution is called an equilibrium density function.

3.1. The Equilibrium Solution and their dynamic behavior

The equilibrium solution $\left(P^{*},P_1^{*},...,P_n^{*}\right)$ of (17) is given by

$$\left\{\begin{array}{c}0=-\left({\mu }_0+\Psi \left(P^{*}\right)\right)P^{*}+R_0\Phi \left(P^{*}\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}^{*}\hfill \\ 0=\left(R_0{\beta }_0\Phi \left(P^{*}\right)-\rho -{\mu }_0-\Psi \left(P^{*}\right)\right)P_1^{*}+R_0\Phi \left(P^{*}\right)\stackrel{n-1}{\underset{i=1}{\Sigma }}{\beta }_i{P}_{i+1}^{*}\hfill \\ 0=i{P}_i^{*}-\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)P_{i+1}^{*},i=1,...,n-1.\hfill \end{array}\right.$$

(19)

To solve the nonlinear algebraic system (19), in our attention is the third equation ($i=1,...,n-1$) from where we obtain successively

$$\begin{array}{ccc}\hfill P_{i+1}^{*}& =& \frac{i}{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{P}_i^{*}\hfill \\ & =& \frac{i\left(i-1\right)}{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{P}_{i-1}^{*}\hfill \\ & =& \\ & & ...\hfill \\ & =& \frac{i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}{P}_1^{*}.\hfill \end{array}$$

The next step, is to replace the determined quantities

$$P_{i+1}^{*}=\frac{i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}{P}_1^{*},i=1,...,n-1.$$

(20)

in the first equation of (19). By equivalence, we obtain

$$\begin{array}{ccc}\hfill 0& =& -\left({\mu }_0+\Psi \left(P^{*}\right)\right)P^{*}+R_0\Phi \left(P^{*}\right){\beta }_0{P}_1^{*}+R_0\Phi \left(P^{*}\right)\stackrel{n-1}{\underset{i=1}{\Sigma }}{\beta }_i{P}_{i+1}^{*}\hfill \\ & \iff & \hfill \\ \hfill 0& =& -\left({\mu }_0+\Psi \left(P^{*}\right)\right)P^{*}+R_0\Phi \left(P^{*}\right){\beta }_0{P}_1^{*}+R_0\Phi \left(P^{*}\right)\stackrel{n-1}{\underset{i=1}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}{P}_1^{*}\hfill \\ & \iff & \hfill \\ \hfill 0& =& -\left({\mu }_0+\Psi \left(P^{*}\right)\right)P^{*}+\left[{\beta }_0+\stackrel{n-1}{\underset{i=1}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}\right]R_0\Phi \left(P^{*}\right)P_1^{*}.\hfill \end{array}$$

Finally, since the non-trivial stationary sizes $P^{*}$ is given by (16), we obtain from equation above that

$$P_1^{*}=\frac{\left({\mu }_0+\Psi \left(P^{*}\right)\right)}{\left[{\beta }_0+\stackrel{n-1}{\underset{i=1}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}\right]R_0\Phi \left(P^{*}\right)}{P}^{*}=\frac{\left({\mu }_0+\Psi \left(P^{*}\right)\right)}{\frac{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{R_0\Phi \left(P^{*}\right)}{R}_0\Phi \left(P^{*}\right)}=\frac{{\mu }_0+\Psi \left(P^{*}\right)}{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{P}^{*}$$

(21)

The existence of a non-trivial stationary solution for the system (17), in the form

$$P_1^{*}=\frac{{\mu }_0+\Psi \left(P^{*}\right)}{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{P}^{*}\mathrm{and}{P}_{i+1}^{*}=\frac{i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}{P}_1^{*},i=1,...,n-1,$$

(22)

is proved if the second equation in (19) is checked by (22). But, this is a simple exercises by the following equivalence

$$\begin{array}{ccc}\hfill 0& =& \left(R_0{\beta }_0\Phi \left(P^{*}\right)-\rho -{\mu }_0-\Psi \left(P^{*}\right)\right)P_1^{*}+R_0\Phi \left(P^{*}\right)\stackrel{n-1}{\underset{i=1}{\Sigma }}{\beta }_i{P}_{i+1}^{*}\hfill \\ & ⟺& \\ \hfill 0& =& \left(R_0{\beta }_0\Phi \left(P^{*}\right)-\rho -{\mu }_0-\Psi \left(P^{*}\right)\right)P_1^{*}-R_0\Phi \left(P^{*}\right){\beta }_0{P}_1^{*}+R_0\Phi \left(P^{*}\right)\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{P}_{i+1}^{*}\hfill \\ & \iff & \\ \hfill 0& =& \left(-\rho -{\mu }_0-\Psi \left(P^{*}\right)\right)P_1^{*}+\left[{\beta }_0+\stackrel{n-1}{\underset{i=1}{\Sigma }}\frac{{\beta }_{i}i!}{{\left(\rho +{\mu }_0+\Psi \left(P^{*}\right)\right)}^i}\right]R_0\Phi \left(P^{*}\right)P_1^{*}\hfill \\ & \iff & \\ \hfill 0& =& \left(-\rho -{\mu }_0-\Psi \left(P^{*}\right)\right)P_1^{*}+\frac{\rho +{\mu }_0+\Psi \left(P^{*}\right)}{R_0\Phi \left(P^{*}\right)}{R}_0\Phi \left(P^{*}\right)P_1^{*}\hfill \end{array}$$

and the last equality is true.

We give an alternative proof for the positivity of the equilibrium. Since $R\left(x\right)$ is a bijection there exists his inverse denoted by $R^{-1}\left(x\right)$. Since

$$R^{-1}\left(x\right)=\frac{1}{R^{\prime }\left(R^{-1}\left(x\right)\right)}$$

and $R^{\prime }\left(t\right)<0$ for all $t\ge 0$ we have that $R^{-1}\left(x\right)<0$ for all $x\ge 0$, i.e. $R^{-1}$ is a decreasing function. Finally, let us observe that (14) implies

$$P^{*}=R^{-1}\left(1\right).$$

(23)

On the other hand from (12) we have

$$R\left(0\right)=R_0⟹0=R^{-1}\left(R_0\right).$$

(24)

Due to the fact that $R_0$$>1$ and $R^{-1}$ is a decreasing function we have

$$0=R^{-1}\left(R_0\right)<R^{-1}\left(1\right)=P^{*}$$

(25)

where we have used (23) and (24). Clearly $P^{*}>0$ imply that $P_i^{*}>0$ for all $i=1,...,n$. Also, (25) proved that: if $R_0$$<1$ then

$$0=R^{-1}\left(R_0\right)>R^{-1}\left(1\right)=P^{*}$$

and so the only trivial equilibrium exists.

We summarize our result in the next Lemmas.

Lemma 2. 

If $R_0>1$, then a unique non-trivial equilibrium solution $\left(P^{*},P_1^{*},...,P_n^{*}\right)$ of (17) exists and is given by (15), (21), and (20).

Lemma 3. 

If $R_0<1$ then the only trivial equilibrium exists.

Let us point that, the stability of the equilibrium point $(P^{*},P_1^{*},...,P_n^{*})$ is determined by the sign of trace of the Jacobian matrix of the system. Depending on the parameter combinations chosen, the model can show stability as well as instability of the non-trivial equilibrium. Also, existence of periodic solutions occurs when passing from one case to the other.

4. Asymptotic behavior at infinity of the solution for the system (17)

In the next section we consider the special case $n=1$. The ODE system becomes

$$\left\{\begin{array}{c}{P}^{\prime }\left(t\right)=-\left({\mu }_0+\Psi \left(P\left(t\right)\right)\right)P\left(t\right)+R_0{\beta }_0\Phi \left(P\left(t\right)\right)P_1\left(t\right)\hfill \\ P_1^{\prime }\left(t\right)=\left(R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\right)P_1\left(t\right)\hfill \end{array}\right.$$

(26)

Assuming $P\left(0\right)>P_1\left(0\right)>0,$ then it follows that $P\left(t\right)>P_1\left(t\right)>0$ according to Proposition 8.5, page 237 in [12].

We have the following results about population behavior at infinity.

Lemma 4. 

If $R_0<1$ then

$$\underset{t\to \infty }{lim}P\left(t\right)=0.$$

(27)

Proof

In light of $P\left(t\right)>P_1\left(t\right)>0$ it follows that

$$P^{\prime }\left(t\right)<\left(R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\right)P\left(t\right).$$

Moreover,

$$R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)<0,$$

in light of

$$R\left(P\left(t\right)\right)\le R_0<1.$$

Thus, $P\left(t\right)$ is decreasing and let

$$\underset{t\to \infty }{lim}P\left(t\right)=L,$$

where $L\ge 0$. Next, we show that $L=0$. The condition $R_0<1$ implies that

$$R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\le (R_0-1)(\rho +{\mu }_0+\Psi \left(P\left(t\right)\right)).$$

Moreover

$$\rho +{\mu }_0+\Psi \left(P\left(t\right)\right)\ge \rho +{\mu }_0+\Psi \left(L\right).$$

Therefore

$$R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)\le (R_0-1)(\rho +{\mu }_0+\Psi \left(L\right))<0.$$

The result follows by a simple application of Gronwall inequality.

Lemma 5. 

If

$$R\left(P\left(t\right)\right)>1\mathrm{for}\mathrm{all}t\ge 0\mathrm{and}{R}_0>1$$

then

$$\underset{t\to \infty }{lim}P\left(t\right)\ne 0.$$

(28)

Proof The condition $R\left(P\left(t\right)\right)>1$ implies that

$$R_0{\beta }_0\Phi \left(P\left(t\right)\right)-\rho -{\mu }_0-\Psi \left(P\left(t\right)\right)>0,$$

and as such $P_1^{\prime }\left(t\right)>0,$ whence $t\to P_1\left(t\right)$ is strictly increasing which combined with $P\left(t\right)>P_1\left(t\right)>0$ yields the result.

Lemma 6. 

If

$$R\left(P\left(t\right)\right)\le 1\mathrm{for}\mathrm{all}t\ge 0\mathrm{and}{R}_0>1$$

then

$$\underset{t\to \infty }{lim}P\left(t\right)\ne 0.$$

(29)

Proof The above condition implies that

$$P\left(t\right)>R^{-1}\left(1\right)=P^{*}>0\mathrm{for}\mathrm{all}t\ge 0,$$

by using the fact that $R^{-1}$ is a decreasing function. This is the case when ${lim}_{t\to \infty }P\left(t\right)\ne 0$ and then the population does not vanish at infinity.

5. Some remarks

About all of the results we have the next two remarks.

Remark 1. 

Our model analysis can be extended to more general mortality

$$\mu (a,p)={\mu }_0\left(t\right)+\Psi \left(p\right),$$

where ${\mu }_0\left(t\right)$ is a deterministic function, i.e., Gompertz function

$${\mu }_0\left(t\right)=a+b{e}^{-ct},$$

for some constants $a,b,$ and c. However in such a case the trivial equilibrium is the only equilibrium.

Remark 2. 

Our approach can be applied to more general fertility function

$$\beta \left(a,p\right)=R_0{\beta }_0\Phi \left(p\right)F\left(a\right),$$

for any continuous function $F\left(a\right).$ Indeed, this is the case since $F\left(a\right)$ can be approximated by

$$\stackrel{n-1}{\underset{i=0}{\Sigma }}{\beta }_i{a}^i{e}^{-\rho a}.$$