A Nonlinear Population Model

Preprint

Article

A Nonlinear Population Model

Altmetrics

Downloads

Views

Comments

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

26 December 2023

Posted:

27 December 2023

You are already at the latest version

Alerts

Abstract

This paper considers a nonlinear model for population dynamics with age structure. The fertility rate with respect to age is non constant and has the form proposed by . Moreover, its multiplicative structure and the multiplicative structure of mortality makes the model separable. In this setting it is shown that the number of births in unit time is given by a system of nonlinear ordinary differential equations. The steady state solution together with the equilibrium solution is found explicitly.

Keywords:

Subject: Computer Science and Mathematics - Applied Mathematics

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 (t)

can be obtained by integrating its density, i.e.,

P (t) = \int_{0}^{\infty} p (a, t) d a .

Next let us introduce fertility and mortality. The number of offspring, borne by individuals during the infinitesimal time-interval

[t, t + d t]

, and the infinitesimal age-interval

[a, a + d a]

β (a, t),

also referred to as the age specific fertility. The number of offspring during the infinitesimal interval

[t, t + d t]

is then

B (t) = \int_{0}^{\infty} β (a, t) p (a, t) d a .

The number of deaths of individuals during the infinitesimal time-interval

[t, t + d t]

, and the infinitesimal age-interval

[a, a + d a]

μ (a, t) .

The number of deaths during the infinitesimal interval

[t, t + d t]

is then

D (t) = \int_{0}^{\infty} μ (a, t) p (a, t) d a .

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

π (a, t, x) = e^{- \int_{0}^{x} μ (a - σ, t - σ) d σ} .

In the case of time independent mortality

π (a) = e^{- \int_{0}^{a} μ (σ) d σ},

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

β (a),

and

μ (a)

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 (σ, t) d σ .

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 (s) d s .

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} μ (σ) p (σ, t + s) d σ d s .

This is the case because

\int_{0}^{a + s} μ (σ) p (σ, t + s) d σ,

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 (s) d s - \int_{0}^{h} \int_{0}^{a + s} μ (σ) p (σ, t + s) d σ d s .

Let us differentiate this with respect to h, and then take

h = 0

to get

p (a, t) + \int_{0}^{a} p_{t} (σ, t) d σ = B (t) - \int_{0}^{a} μ (σ) p (σ, t) d σ .

Next, let us differentiate this with respect to

a,

to get

p (a, t) + p_{a} (a, t) + μ (a) p (a, t) = 0 .

Also by setting

a = 0

yields

p (0, t) = B (t),

but on the other hand

B (t) = \int_{0}^{\infty} β (σ) p (σ, t) d σ .

As such we obtained the following system

\{\begin{matrix} p_{a} (a, t) + p_{t} (a, t) + μ (a) p (a, t) = 0 & for & a, t \geq 0, \\ p (0, t) = \int_{0}^{\infty} β (σ) ρ (σ, t) d σ & for & t > 0, \\ p (a, 0) = p_{0} (a) & for & a > 0 . \end{matrix}

(1)

2.1. A Special Linear Population Model

Inspired by [1,7] and [8], the recent work [4] considers a model with survival probability

π (a)

pseudo exponential, i.e.

π (a) = \sum_{i = 1}^{n} c_{i} e^{- μ_{i} a},

(2)

for positive constants

c_{i}, μ_{i}

with

\sum_{i = 1}^{n} c_{i}^{2} \neq 0, \sum_{i = 1}^{n} μ_{i}^{2} \neq 0 and \sum_{i = 1}^{n} c_{i} = 1 .

Moreover, we present the case of

π (a) = (\sum_{i = 0}^{n} c_{i} a^{i}) e^{- μ_{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

μ (a, p),

and the fertility rate

β (a, p),

are functions of age and population size. The age density function at time

t \geq 0,

p (a, t)

satisfies the following nonlinear equations

\{\begin{matrix} p_{t} (a, t) + p_{a} (a, t) + μ (a, P (t)) p (a, t) = 0 \\ p (0, t) = \int_{0}^{\infty} β (σ, P (t)) p (σ, t) d σ \\ p (a, 0) = p_{0} (a) \\ P (t) = \int_{0}^{\infty} p (σ, t) d σ, \end{matrix}

where

P (t)

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

π (a, t, x, P) = e^{- \int_{0}^{x} μ (a - σ, P (t - σ)) d σ}

and

p (a, t) = \{\begin{matrix} p_{0} (t - a) π (a, t, t, P) & f o r & t \leq a \\ B (t - a) π (a, t, a, P) & f o r & t > a . \end{matrix}

(4)

Here the function B solves the following system of equations:

\{\begin{matrix} B (t) = \int_{0}^{t} β (σ, P (t)) π (σ, t, σ, P) B (σ) d σ + F (t, P) \\ P (t) = \int_{0}^{t} π (σ, t, σ, P) B (σ) d σ + G (t, P) \end{matrix}

(5)

where

\{\begin{matrix} F (t, P) = \int_{t}^{\infty} β (σ, P (t)) π (σ, t, t, P) p_{0} (σ - t) d σ \\ G (t, P) = \int_{t}^{\infty} π (σ, t, t, P) p_{0} (σ - t) d σ . \end{matrix}

(6)

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

\{\begin{matrix} B^{k + 1} (t) = \int_{0}^{t} β (σ, P^{k} (t)) π (σ, t, σ, P^{k}) B^{k} (σ) d σ + F (t, P^{k}) \\ P^{k + 1} (t) = \int_{0}^{t} π (σ, t, σ, P^{k}) B^{k} (σ) d σ + G (t, P^{k}), \end{matrix}

(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

β (a, p) = R_{0} β_{0} Φ (p) e^{- ρ a}, μ (a, p) = μ_{0} + Ψ (p) .

By plugging this in the nonlinear system one gets

\{\begin{matrix} p_{t} (a, t) + p_{a} (a, t) + (μ_{0} + Ψ (P (t))) (a, t) = 0 \\ p (0, t) = \int_{0}^{\infty} R_{0} Φ (P (t)) e^{- ρ a} p (σ, t) d σ \\ p (a, 0) = p_{0} (a) \\ P (t) = \int_{0}^{\infty} p (σ, t) d σ . \end{matrix}

These nonlinear equations can be reduced to the following ODE system

\{\begin{matrix} P^{'} (t) = - (μ_{0} + Ψ (P (t))) P (t) + R_{0} β_{0} Φ (P (t)) Q (t) \\ Q^{'} (t) = (R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t))) Q (t) \end{matrix}

(8)

where

Q (t) = \int_{0}^{\infty} e^{- a ρ} p_{t} (a, t) d a,

see [12] for more details.

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

R (x) = R_{0} Φ (x) \int_{0}^{\infty} β_{0} e^{- (ρ + μ_{0} + Ψ (x)) a} d a .

This quantity was introduced first by [10]. According with [10] the quantity

R (x)

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 (P^{*}) = 1,

which in this setting reads

R_{0} β_{0} Φ (P^{*}) = ρ + μ_{0} + Ψ (P^{*}) .

3. A Special Nonlinear Population Model

We consider the following logistic system

\{\begin{matrix} p_{t} (a, t) + p_{a} (a, t) + μ_{0} p (a, t) + Ψ (P (t)) p (a, t) = 0 \\ p (0, t) = R_{0} Φ (P (t)) \int_{0}^{\infty} \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} σ^{i} e^{- ρ σ} p (σ, t) d σ \\ p (a, 0) = p_{0} (a) \\ P (t) = \int_{0}^{\infty} p (σ, t) d σ \end{matrix}

(9)

where

p (a, t)

is the age density function at time

t \geq 0

P (t)

is the total population at time t,

p_{0} (a)

is the initial population at time

t = 0

a \in [0, \infty)

μ_{0} > 0

is the intrinsic mortality term,

R_{0} > 0

is the birth modulus and

ρ

σ

are prescribed positive parameters. In real world, the age profile is given by

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

The fertility is given by

β (a, p) = R_{0} Φ (p) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} a^{i} e^{- ρ a} with β_{0}, β_{1}, . . ., β_{n - 1} \in (0, \infty),

and the mortality (death rate) by

μ (a, P) = μ_{0} + Ψ (P) .

The couple

(β (a), μ_{0})

can be interpreted as an intrinsic birth-death process that is age dependent while

Ψ (P (t))

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

Φ

and

Ψ

are continuous on

[0, \infty)

continuously differentiable on

(0, \infty)

and

\begin{matrix} Φ (x) & \geq & 0, Φ^{'} (x) < 0, Φ (0) = 1, Φ (+ \infty) = 0, \end{matrix}

(10)

\begin{matrix} Ψ (x) & \geq & 0, Ψ^{'} (x) > 0, Ψ (0) = 0, Ψ (+ \infty) = + \infty . \end{matrix}

(11)

We also adopt the normalization condition

\int_{0}^{\infty} \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} a^{i} e^{- (ρ + μ_{0}) a} d a = 1,

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

\overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0})}^{i + 1}} = 1,

(12)

so that the parameter

R_{0}

in (9) has the role of an intrinsic basic reproduction number denoted by

R (x)

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{matrix} R (x) & = & R_{0} Φ (x) \int_{0}^{\infty} \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} a^{i} e^{- (ρ + μ_{0} + Ψ (x)) a} d a \\ = & R_{0} Φ (x) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} \int_{0}^{\infty} a^{i} e^{- (ρ + μ_{0} + Ψ (x)) a} d a \\ = & R_{0} Φ (x) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} \int_{0}^{\infty} {(\frac{t}{ρ + μ_{0} + Ψ (x)})}^{i} e^{- t} \frac{d t}{ρ + μ_{0} + Ψ (x)} \\ = & R_{0} Φ (x) \overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i}}{{(ρ + μ_{0} + Ψ (x))}^{i + 1}} \int_{0}^{\infty} t^{i} e^{- t} d t \\ = & R_{0} Φ (x) \overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i} Γ (i + 1)}{{(ρ + μ_{0} + Ψ (x))}^{i + 1}} \\ = & R_{0} Φ (x) \overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (x))}^{i + 1}} \end{matrix}

where

Γ

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

lim_{x \to \infty} R (x) = 0 and R (x) is a decreasing function (i . e . R^{'} (x) < 0) .

(13)

Indeed

R (x)

is a decreasing function because

R^{'} (x) = \overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i} i! [R_{0} Φ^{'} (x) (ρ + μ_{0} + Ψ (x)) - R_{0} Φ (x) Ψ^{'} (x)]}{{(ρ + μ_{0} + Ψ (x))}^{i + 2}} - \overset{n - 1}{\underset{i = 0}{Σ}} \frac{i β_{i} i! R_{0} Φ (x) Ψ^{'} (x)}{{(ρ + μ_{0} + Ψ (x))}^{i + 2}} < 0,

for all

x \geq 0

. Moreover,

lim_{x \to \infty} R (x) = 0,

is satisfied in light of the asymptotic conditions on

Ψ, Φ .

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 (P^{*}) = 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} Φ (P^{*})}{ρ + μ_{0} + Ψ (P^{*})} \overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{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

\overset{n - 1}{\underset{i = 0}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}} = \frac{ρ + μ_{0} + Ψ (P^{*})}{R_{0} Φ (P^{*})} .

(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 (t) .

We will denote by

P_{i} (t) = \int_{0}^{\infty} σ^{i - 1} e^{- ρ σ} p (σ, t) d σ,

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{matrix} p (0, t) & = & R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} \int_{0}^{\infty} σ^{i} e^{- ρ σ} p (σ, t) d σ \\ = & R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1} (t) . \end{matrix}

More that, the calculation of the first derivative of

P (t)

gives

\begin{matrix} P^{'} (t) & = & \int_{0}^{\infty} p_{t} (a, t) d a \\ = & - \int_{0}^{\infty} p_{a} (a, t) d a - (μ_{0} + Ψ (P (t))) \int_{0}^{\infty} p (a, t) d a \\ = & p (0, t) - (μ_{0} + Ψ (P (t))) P (t) \\ = & - (μ_{0} + Ψ (P (t))) P (t) + R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1} (t) \end{matrix}

and, similarly

\begin{matrix} P_{1}^{'} (t) & = & \int_{0}^{\infty} e^{- a ρ} p_{t} (a, t) d a \\ = & - \int_{0}^{\infty} e^{- a ρ} p_{a} (a, t) d a - (μ_{0} + Ψ (P (t))) \int_{0}^{\infty} e^{- a ρ} p (a, t) d a \\ = & p (0, t) - ρ \int_{0}^{\infty} e^{- a ρ} p_{a} (a, t) d a - (μ_{0} + Ψ (P (t))) P_{1} (t) \\ = & R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1} - ρ P_{1} (t) - (μ_{0} + Ψ (P (t))) P_{1} (t) \\ = & (R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t))) P_{1} (t) + R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 1}{Σ}} β_{i} P_{i + 1} . \end{matrix}

In the same way, the calculation of the first derivative of

P_{i + 1} (t)

(

i = 1, . . ., n - 1

) gives

\begin{matrix} P_{i + 1}^{'} (t) & = & \int_{0}^{\infty} a^{i} e^{- ρ a} p_{t} (σ, t) d σ \\ = & - \int_{0}^{\infty} a^{i} e^{- a ρ} p_{a} (a, t) d a - (μ_{0} + Ψ (P (t))) \int_{0}^{\infty} a^{i} e^{- a ρ} p (a, t) d a \\ = & \int_{0}^{\infty} (i a^{i - 1} e^{- a ρ} - a^{i} ρ e^{- a ρ}) p (a, t) d a - (μ_{0} + Ψ (P (t))) \int_{0}^{\infty} a^{i} e^{- a ρ} p (a, t) d a \\ = & i \int_{0}^{\infty} a^{i - 1} e^{- a ρ} p (a, t) d a - \int_{0}^{\infty} a^{i} ρ e^{- a ρ} p (a, t) d a - (μ_{0} + Ψ (P (t))) \int_{0}^{\infty} a^{i} e^{- a ρ} p (a, t) d a \\ = & i P_{i} (t) - ρ P_{i + 1} (t) - (μ_{0} + Ψ (P (t))) P_{i + 1} (t) \\ = & i P_{i} (t) - (ρ + μ_{0} + Ψ (P (t))) P_{i + 1} (t) . \end{matrix}

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

\{\begin{matrix} P^{'} (t) = - (μ_{0} + Ψ (P (t))) P (t) + R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1} \\ P_{1}^{'} (t) = (R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t))) P_{1} (t) + R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 1}{Σ}} β_{i} P_{i + 1} \\ P_{i + 1}^{'} (t) = i P_{i} (t) - (ρ + μ_{0} + Ψ (P (t))) P_{i + 1} (t), i = 1, . . ., n - 1 and n \geq 2 \end{matrix}

(17)

coupled with the initial conditions

P (0) = P_{0}, P_{i} (0) = P_{0}^{i}, i = 1, 2, . . ., n

(18)

where

P (0) = \int_{0}^{\infty} p_{0} (a) d a, P_{i} (0) = \int_{0}^{\infty} a^{i - 1} e^{- ρ a} p_{0} (a) d a, 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 (t), P_{1} (t), . . ., P_{n} (t))

solves (17), then by setting

B (t) = p (0, t) = R_{0} Φ (P (t)) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1} (t)

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

p (a, t) = \{\begin{matrix} p_{0} (a - t) e^{- \int_{0}^{t} (μ_{0} + Ψ (P (σ))) d σ} & a \geq t, \\ B (t - a) e^{- \int_{t - a}^{t} (μ_{0} + Ψ (P (σ))) d σ} & a < t . \end{matrix}

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

(P^{*}, P_{1}^{*}, . . ., P_{n}^{*})

of (17) is given by

\{\begin{matrix} 0 = - (μ_{0} + Ψ (P^{*})) P^{*} + R_{0} Φ (P^{*}) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1}^{*} \\ 0 = (R_{0} β_{0} Φ (P^{*}) - ρ - μ_{0} - Ψ (P^{*})) P_{1}^{*} + R_{0} Φ (P^{*}) \overset{n - 1}{\underset{i = 1}{Σ}} β_{i} P_{i + 1}^{*} \\ 0 = i P_{i}^{*} - (ρ + μ_{0} + Ψ (P^{*})) P_{i + 1}^{*}, i = 1, . . ., n - 1 . \end{matrix}

(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{matrix} P_{i + 1}^{*} & = & \frac{i}{ρ + μ_{0} + Ψ (P^{*})} P_{i}^{*} \\ = & \frac{i (i - 1)}{ρ + μ_{0} + Ψ (P^{*})} P_{i - 1}^{*} \\ = \\ . . . \\ = & \frac{i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}} P_{1}^{*} . \end{matrix}

The next step, is to replace the determined quantities

P_{i + 1}^{*} = \frac{i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}} P_{1}^{*}, i = 1, . . ., n - 1 .

(20)

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

\begin{matrix} 0 & = & - (μ_{0} + Ψ (P^{*})) P^{*} + R_{0} Φ (P^{*}) β_{0} P_{1}^{*} + R_{0} Φ (P^{*}) \overset{n - 1}{\underset{i = 1}{Σ}} β_{i} P_{i + 1}^{*} \\ \Leftrightarrow \\ 0 & = & - (μ_{0} + Ψ (P^{*})) P^{*} + R_{0} Φ (P^{*}) β_{0} P_{1}^{*} + R_{0} Φ (P^{*}) \overset{n - 1}{\underset{i = 1}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}} P_{1}^{*} \\ \Leftrightarrow \\ 0 & = & - (μ_{0} + Ψ (P^{*})) P^{*} + [β_{0} + \overset{n - 1}{\underset{i = 1}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}}] R_{0} Φ (P^{*}) P_{1}^{*} . \end{matrix}

Finally, since the non-trivial stationary sizes

P^{*}

is given by (16), we obtain from equation above that

P_{1}^{*} = \frac{(μ_{0} + Ψ (P^{*}))}{[β_{0} + \overset{n - 1}{\underset{i = 1}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}}] R_{0} Φ (P^{*})} P^{*} = \frac{(μ_{0} + Ψ (P^{*}))}{\frac{ρ + μ_{0} + Ψ (P^{*})}{R_{0} Φ (P^{*})} R_{0} Φ (P^{*})} = \frac{μ_{0} + Ψ (P^{*})}{ρ + μ_{0} + Ψ (P^{*})} P^{*}

(21)

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

P_{1}^{*} = \frac{μ_{0} + Ψ (P^{*})}{ρ + μ_{0} + Ψ (P^{*})} P^{*} and P_{i + 1}^{*} = \frac{i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{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{matrix} 0 & = & (R_{0} β_{0} Φ (P^{*}) - ρ - μ_{0} - Ψ (P^{*})) P_{1}^{*} + R_{0} Φ (P^{*}) \overset{n - 1}{\underset{i = 1}{Σ}} β_{i} P_{i + 1}^{*} \\ ⟺ \\ 0 & = & (R_{0} β_{0} Φ (P^{*}) - ρ - μ_{0} - Ψ (P^{*})) P_{1}^{*} - R_{0} Φ (P^{*}) β_{0} P_{1}^{*} + R_{0} Φ (P^{*}) \overset{n - 1}{\underset{i = 0}{Σ}} β_{i} P_{i + 1}^{*} \\ \Leftrightarrow \\ 0 & = & (- ρ - μ_{0} - Ψ (P^{*})) P_{1}^{*} + [β_{0} + \overset{n - 1}{\underset{i = 1}{Σ}} \frac{β_{i} i!}{{(ρ + μ_{0} + Ψ (P^{*}))}^{i}}] R_{0} Φ (P^{*}) P_{1}^{*} \\ \Leftrightarrow \\ 0 & = & (- ρ - μ_{0} - Ψ (P^{*})) P_{1}^{*} + \frac{ρ + μ_{0} + Ψ (P^{*})}{R_{0} Φ (P^{*})} R_{0} Φ (P^{*}) P_{1}^{*} \end{matrix}

and the last equality is true.

We give an alternative proof for the positivity of the equilibrium. Since

R (x)

is a bijection there exists his inverse denoted by

R^{- 1} (x)

. Since

R^{- 1} (x) = \frac{1}{R^{'} (R^{- 1} (x))}

and

R^{'} (t) < 0

for all

t \geq 0

we have that

R^{- 1} (x) < 0

for all

x \geq 0

, i.e.

R^{- 1}

is a decreasing function. Finally, let us observe that (14) implies

P^{*} = R^{- 1} (1) .

(23)

On the other hand from (12) we have

R (0) = R_{0} ⟹ 0 = R^{- 1} (R_{0}) .

(24)

Due to the fact that

R_{0}

> 1

and

R^{- 1}

is a decreasing function we have

0 = R^{- 1} (R_{0}) < R^{- 1} (1) = 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} (R_{0}) > R^{- 1} (1) = P^{*}

and so the only trivial equilibrium exists.

We summarize our result in the next Lemmas.

Lemma 2.

R_{0} > 1

, then a unique non-trivial equilibrium solution

(P^{*}, P_{1}^{*}, . . ., P_{n}^{*})

of (17) exists and is given by (15), (21), and (20).

Lemma 3.

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

\{\begin{matrix} P^{'} (t) = - (μ_{0} + Ψ (P (t))) P (t) + R_{0} β_{0} Φ (P (t)) P_{1} (t) \\ P_{1}^{'} (t) = (R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t))) P_{1} (t) \end{matrix}

(26)

Assuming

P (0) > P_{1} (0) > 0,

then it follows that

P (t) > P_{1} (t) > 0

according to Proposition 8.5, page 237 in [12].

We have the following results about population behavior at infinity.

Lemma 4.

R_{0} < 1

then

lim_{t \to \infty} P (t) = 0 .

(27)

Proof

In light of

P (t) > P_{1} (t) > 0

it follows that

P^{'} (t) < (R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t))) P (t) .

Moreover,

R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t)) < 0,

in light of

R (P (t)) \leq R_{0} < 1 .

Thus,

P (t)

is decreasing and let

lim_{t \to \infty} P (t) = L,

where

L \geq 0

. Next, we show that

L = 0

. The condition

R_{0} < 1

implies that

R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t)) \leq (R_{0} - 1) (ρ + μ_{0} + Ψ (P (t))) .

Moreover

ρ + μ_{0} + Ψ (P (t)) \geq ρ + μ_{0} + Ψ (L) .

Therefore

R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t)) \leq (R_{0} - 1) (ρ + μ_{0} + Ψ (L)) < 0 .

The result follows by a simple application of Gronwall inequality.

Lemma 5.

R (P (t)) > 1 for all t \geq 0 and R_{0} > 1

then

lim_{t \to \infty} P (t) \neq 0 .

(28)

Proof The condition

R (P (t)) > 1

implies that

R_{0} β_{0} Φ (P (t)) - ρ - μ_{0} - Ψ (P (t)) > 0,

and as such

P_{1}^{'} (t) > 0,

whence

t \to P_{1} (t)

is strictly increasing which combined with

P (t) > P_{1} (t) > 0

yields the result.

Lemma 6.

R (P (t)) \leq 1 for all t \geq 0 and R_{0} > 1

then

lim_{t \to \infty} P (t) \neq 0 .

(29)

Proof The above condition implies that

P (t) > R^{- 1} (1) = P^{*} > 0 for all t \geq 0,

by using the fact that

R^{- 1}

is a decreasing function. This is the case when

{lim}_{t \to \infty} P (t) \neq 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

μ (a, p) = μ_{0} (t) + Ψ (p),

where

μ_{0} (t)

is a deterministic function, i.e., Gompertz function

μ_{0} (t) = a + b e^{- c t},

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

β (a, p) = R_{0} β_{0} Φ (p) F (a),

for any continuous function

F (a) .

Indeed, this is the case since

F (a)

can be approximated by

\overset{n - 1}{\underset{i = 0}{Σ}} β_{i} a^{i} e^{- ρ a} .

Funding

This research was supported by NSERC grant 5-36700 by Traian A. Pirvu and Horizon 2020-2017-RISE-777911 project by Catalin Sterbeti.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Beghriche, A.; Zeghdoudi, H.; Raman, V.; Chouia, S. New polynomial exponential distribution: properties and applications. Statistics in Transition new series. 2022, 23(3), 95–112. [Google Scholar] [CrossRef]
Brauer, F.; Castillo–Chavez, C. Mathematical Models in Population Biology and Epidemiology (Second Edition), Texts in Applied Mathematics. 40, Springer, New York, 2012.
Brauer, F.; Castillo-Chavez, C.; Feng, Z. Mathematical Models in Epidemiology (Second Edition), Texts in Applied Mathematics. 69, New York, 2019.
Covei, D. P.; Pirvu T. A..; Sterbeti C. A population model with pseudo exponential survival, 2023, Preprint.
Covei, D. P. Linear algebra elements, ASE Publishing House, p. 1-181, 2015.
Hoppensteadt, F.C. Mathematical Theories of Populations: Demographics, Genetics and Epidemics. SIAM, Philadelphia, 1975.
Dufresne, D. Fitting combinations of exponentials to probability distributions. Appl. Stochastic Models Bus. Ind.
Feller, W. On the integral equation of renewal theory. Ann. Math. Statist. 1941, 12(3), 243–267. [Google Scholar] [CrossRef]
Györi, I. Some mathematical aspects of modelling cell population dynamics. Computers Math. Applic. 1990, 20, 127–138. [Google Scholar] [CrossRef]
Gurtin, M.E.; Maccamy, R.C. Non-linear age-dependent population dynamics. Arch. Rational Mech. Anal. 1974, 54, 281–300. [Google Scholar] [CrossRef]
Iannelli, M. Mathematical Theory of Age-Structured Population Dynamics. Appl. Mathematical Monographs, C.N.R., 1995.
Iannelli, M.; Milner, F. The Basic Approach to Age-Structured Population Dynamics, Springer, Dordrecht, The Netherlands, 2017.
Keyfitz, B. L.; Keyfitz, N. The McKendrick Partial Differential Equation and Its Uses in Epidemiology and Population Study. Alathl. Comput. Modelling. 1997, 26(6), l–9. [Google Scholar] [CrossRef]
McKendrick, A.G. Applications of mathematics to medical problems. Proc. Edinburgh Math. Soc. 1926, 44, 98–130. [Google Scholar] [CrossRef]
Levine, D. Models og age-dependent predation and cannibalism via the McKendrick equation. Comp. & Math with Appls 1983, 9, 403–414. [Google Scholar]
Sharpe, F.R.; Lotka, A.J. A problem in age distribution, The London, Edinburgh and Dublin Philosophical Magazine and Journal of Science. 1911, 21(124), 435-438.
Sterbeti, C. On a model for population with age structure, ITM Web of Conferences, 2020, vol. 34, 02010.
Swart, J., H. , and Meijer A. , R., A Simplified Model for Age-Dependent Population Dynamics Mathematical Biosciences. 1994, 121, 15–36. [Google Scholar]
Webb, G. Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985.
Weinstock, E.; Rorres, C. , Local Stability of an Age-Structured Population with Density-Dependent Fertility and Mortality. SIAM J. APPL. MATH. 1987, 47, 589–604. [Google Scholar] [CrossRef]

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

A Nonlinear Population Model

Abstract

1. Introduction

2. Population Models

2.1. A Special Linear Population Model

2.2. Nonlinear Population Models

2.3. Separable Population Models

3. A Special Nonlinear Population Model

3.1. The Equilibrium Solution and their dynamic behavior

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

5. Some remarks

Funding

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe