Stability of nonlinear population systems with individual scale and migration

1.   Introduction

In recent years, the status of protecting the environment and population diversity has been improved, and environmental pollution has become the top priority of environmental issues. Our primary concern is the interaction between environmental pollution and biological populations. We need to consider the interdependence between biological populations as well as the interdependence between populations and the environment. In the dynamical model of a biological population, the equilibrium of the model corresponds to the static equilibrium of the biological system, and the stability of equilibrium describes the stability of ecological balance and represents the anti-interference ability of the ecosystem. Therefore, it is paramount to study the stability of the population dynamics model.

So far, scholars at home and abroad have established a large number of mathematical models in the process of exploring the evolution of biological populations. Age is an important parameter, affecting the size and behavior of the population. Sharpe and Lotka ^[1] established a linear model system with an age structure. Since the birth and death rate of individuals are always affected by the number of populations, Curtin and Maccamy ^[2] established the age-structured population model with density constraint, proved the existence and uniqueness of the model, and obtained that the equilibrium is locally asymptotically stable. Considering the environmental constraints on population growth, a nonlinear age-structured population model was studied^[3]. For the population model with age structure, many scholars have studied its optimal control problem^[4,5,6,7]. In^[8], the nonlinear age-structured population model with nonlocal diffusion and boundary conditions was studied. Lopez et al. ^[9] studied the discrete nonlinear age-structured population model of environmental changes and obtained conditions for the existence and asymptotic stability of positive equilibrium. Kilic et al. ^[10] used the forward Euler algorithm to study the discrete-time prey-predator population model with migration, and obtained the asymptotic stability of the unique positive fixed point.

Studies by ecologists show that the scale structure is a continuous variable describing the individual characteristics of the population (volume, diameter, length, maturity, etc.), it can better reflect the relationship between groups and individuals in the evolution of biological populations. So it has attracted extensive attention from ecologists and mathematicians ^[11,12,13]. Dercole et al. ^[14] explored a nonlinear plant growth model in which the scale variable represented the base diameter of the tree, and obtained the stable properties of the equilibrium. Calsina et al. ^[15,16] considered environmental factors, they studied the scale structure models when individual growth rate, death rate, and birth rate were related to the environment and the total population. Farkas ^[17,18,19] successively studied a class of nonlinear and migration models based on scale structure, in which mortality rate and birth rate were related to the total population, and obtained the stability conditions of the equilibrium. Kang et al.^[20] studied a nonlinear double physiologically structured population models with two internal variables and obtained the stability of the steady state. Ademosu et al. ^[21] analyzed the stability and optimal measure for controlling eco-epidemiological dynamics of the prey-predator model. Singh and Emerick ^[22] studied generalized stability conditions for host-parasitoid population dynamics and identified implications for biological control. Molla et al.^[23] analyzed a mathematical model on predator-prey interactions incorporating prey refuge and additive Allee effect on the prey species and obtained the various dynamical behaviors. Wu^[24] studied an optimal harvesting problem of a size-structured population model:

where $p(x, t)$ represented the number of individuals with scale $x$ at time $t,$ $g(x)$ represented the growth rate, $\mu, \beta$ represented individual natural mortality rate and reproductive rate respectively, $p_{0}(x)$ represented the initial distribution of the population, $q_{i}(x), i = 1, 2$ were different weight functions used to reflect the dependence of the lifetime parameters on density, $u(x)$ represented the harvest intensity of humans to the individual population.

With the development of industry and the intensification of human activities, a large number of toxic substances have been discharged into the ecological environment, destroying the habitat and community structure of species, and seriously threatening the survival of the species exposed to the environment (humans, fish, plankton, etc.). Liu et al. ^[25] studied the persistence and extinction of a single-species population system with random disturbance and pulse toxicity in a polluted environment. Qi et al.^[26] studied the dynamic behavior of a stage-structured population model with transient and non-transient pulse effects in a polluted environment. Wei et al.^[27] studied the psychological effect on single-species population model in a polluted environment and obtained the conditions to ensure local extinction and average persistence of the population. Ma et al.^[28] studied the stationary distribution and optimal control of a stochastic population model in a polluted environment.

However, there is little literature on the stability of scale-structured populations in a polluted environment. Considering the biological background of the population model, this paper focuses on the external living environment of species and the impact of migration on species survival, so that the bio-mathematics model is more realistic. Therefore, based on reference ^[24], we study the following population model with scale structure and migration in a polluted environment. Generally speaking, the increase in population has different effects on the birth and death rates of individuals. We consider the density restriction in the process of population evolution with scale structure, and adopt different individual weighting functions to represent the different ways in which life parameters depend on density. The model is as follows:

where the fixed constant $m > 0,$ which represents the largest scale of the individual, $0 < T\leq\infty,$ $T$ represents the maximum time of the evolution. The variables of the model are described in Table 1.

Note 1: In system (1.1), $l$ represents the individual scale, here denotes as $s.$ In system (1.2), $l$ becomes $m,$ and $m$ represents the largest scale of the individual.

Note 2: We consider the external migration of individuals, which is assumed as $f(s)\geq0.$ In a certain environment, due to limited resources, the birth rate will be affected by the proportion of females in the population density $\gamma(s)$, that is, the birth rate will decrease with the increase of the population density, and the mortality rate will be affected by the weight function of total population $\delta(s),$ that is, the mortality rate increases with the increase of population density. This is reasonable in the actual external environment.

We make the following assumptions:

$(H_1)$ $f(s)$ is a non-negative continuous function;

$(H_2)\; \; \tilde{\mu}(s, t): = \mu(s, c_0(t),$ $J(t))\in L_{loc}^{1}((0, m)\times(0, T)),$ $0\leq\mu(s, c_{0}(t), J(t))\leq \mu^{0}, \; \forall x\in(0, +\infty), t\in (0, T),$ $\mu(\cdot, c_{0}(t), x)$ is monotonically increasing with respect to $x$ and locally Lipschitz continuous $|\mu(s, c_{0}(t), x_{1})-\mu(s, c_{0}(t), x_{2})|\leq L(n_{1})|x_{1}-x_{2}|, L(n_{1})$ is the Lipschitz constant;

$(H_3)\; \tilde{\beta}(s, t): = \beta(s, c_0(t), R(t))\in L_{loc}^{1}((0, m)\times(0, T)),$ $0\leq \beta(s, c_{0}(t), R(t))\leq \beta^{0},$ $\forall x\in(0, +\infty), t\in(0, T),$ $\beta(\cdot, c_{0}(t), x)$ is monotonically decreasing with respect to $x$ and locally Lipschitz continuous $|\beta(s, c_{0}(t), x_{1})-\beta(s, c_{0}(t), x_{2})|\leq L(n_{2})|x_{1}-x_{2}|, L(n_{2})$ is the Lipschitz constant;

$(H_4)\; \; g\in C^{1}([0, m])$ and $g(s)$ is continuously differentiable, $\frac{1}{g(\cdot)}\in L^{1}[0, m],$ $\frac{\delta(\cdot)}{g(\cdot)}\in L^{1}[0, m],$ $\frac{\gamma(\cdot)}{g(\cdot)}\in L^{1}[0, m],$ $\frac{\beta(\cdot, 0, 0)}{g(\cdot)}\in L^{1}[0, m];$

$(H_5)\; \; 0\leq p_{0}(s)\leq M_{1},$ $\forall s\in [0, m); \lim\limits_{s\rightarrow m}p_{0}(s) = 0;$

$(H_6)\; \; g\leq k\leq g+m, v(t)\leq h.$

This paper mainly studies the stability of equilibrium. Over time, the population either persists or becomes extinct. The extinction of a population means that the total number of individuals in the population tends to zero in a finite or infinite time interval. However, with the continuous development of the population, does the population have one or more equilibria? If so, under what conditions are these equilibria stable or unstable? These problems have guiding significance for the sustainability of ecological balance.

Note 3: In ^[29], we have established the well-posedness of the system, the system owns a unique solution, which is non-negative, bounded and continuous with respect to the initial distributions. The method of proving the existence and uniqueness of model (1.2) in this paper is similar to that in ^[29].

2.   The existence and uniqueness of the equilibrium solution

When discussing stability, it is assumed that there is no human interference, namely $u\equiv0.$ Since the second and the third equations in system (1.2) are two first order ordinary differential equations, which could be solved by method of variation of constant, we get

For system (1.2), when it satisfies the condition $0 < k_{1} < g_{1}+n, \sup v(t)\leq h,$ for $t\in(0, T),$ $0 \leq c_{0}(t)\leq 1, 0\leq c_{e}(t)\leq1$ is established.

If there is an equilibrium solution $p_{*}(s)$ in model (1.2), the equilibrium solution satisfies the following equation:

From the first equation of (2.1), it can be determined

Since

let

solved by the above formula

Substituting the above formula into (2.1), it can be determined

The net regeneration number of this population is defined as

According to (2.3), $N(c_{0}(t), R, J)$ must be a normal number less than 1.

Lemma 1. If (2.3) has a unique positive solution $(p_{*}(0), c_{0*}, c_{e*}, J_{*}, R_{*}),$ then system (2.1) has a unique non-trivial equilibrium state $p_{*}(s).$

Next, we discuss the existence and uniqueness of the positive solutions of (2.3).

Lemma 2. If $N(0, 0, 0) < 1,$ $\exists\; M_{2}, M_{3}, M_{4},$ such that $0\leq p_{*}(0)\leq M_{2}, 0\leq c_{0*}\leq1, 0\leq c_{e*}\leq1, 0\leq J_{*}\leq M_{3}, 0\leq R_{*}\leq M_{4},$ where

Proof. Let

It can be known by $(H_{2}), (H_{3})$ and the monotonicity of $T(c_{0}(t), J, R)$, $F(s, c_{0}(t), J, R)$ that

Then, we obtain by the first formula of (2.3) that

therefore

From the fourth formula of (2.3), it can be determined that

From the fifth formula of (2.3), it can be determined that

Lemma 3. Set $\Omega = [0, M_{2}]\times[0, 1]\times[0, 1]\times[0, M_{3}]\times[0, M_{4}]\subset {\mathbb R}^5,$ $A = (a_{1}, a_{2}, a_{3}, a_{4}, a_{5})^{{\rm T}}\in\Omega,$ and define the mapping $\phi:\Omega\subset {\mathbb R}^5\rightarrow {\mathbb R}^5,$ by $\phi A = Y,$ where $Y = (y_{1}, y_{2}, y_{3}, y_{4}, y_{5})^{{\rm T}},$

Then, there is a constant $\overline{M},$ such that $\forall A^{1}, A^{2}\in \Omega,$ $\|\phi A^{1}-\phi A^{2}\|\leq\overline{M}\|A^{1}-A^{2}\|,$ where

Proof. By assumptions $(H_2), (H_3), (H_5),$ when $|x_{i}|\leq \xi, i = 1, 2, \xi = \min\{M_{2}, M_{3}, M_{4}\},$

Firstly, since $A\in\Omega,$ there is

since $M_{2} = \frac{\frac{1}{g(0)}\|f\|_{L^{1}}\|\frac{\beta(\cdot, 0, 0)}{g(\cdot)}\|_{L^{1}}}{1-\|\frac{\beta(\cdot, 0, 0)}{g(\cdot)}\|_{L^{1}}},$ therefore, $0\leq y_{1}\leq M_{2}.$

Since $0\leq c_{0*}\leq1, 0\leq c_{e*}\leq1,$ therefore, $0\leq y_{2}\leq1, 0\leq y_{3}\leq1.$

since $M_{3} = (g(0)M_{2}+\|f\|_{{L^{1}}})\|\frac{\delta(\cdot)}{g(\cdot)}\|_{{L^{1}}},$ therefore, $0\leq y_{4}\leq M_{3}.$

since $M_{4} = (g(0)M_{2}+\|f\|_{{L^{1}}})\|\frac{\gamma(\cdot)}{g(\cdot)}\|_{{L^{1}}},$ therefore, $0\leq y_{5}\leq M_{4}.$

That is $Y = (y_1, y_2, y_3, y_4, y_5)^{{\rm T}}\in\Omega.$ Secondly, it can be proved that $\phi$ maps $\Omega$ into itself, let $\phi A^{i} = Y^{i},$ where

It can be determined from the first formula of (2.6) that

It can be determined from the second formula of (2.6) that

It can be determined from the third formula of (2.6) that

It can be determined from the fourth formula of (2.6) that

It can be determined from the fifth formula of (2.6) that

It can be known by $(2.7)$–$(2.11)$ that

Since $\overline{M} = \max\{N_{1}, N_{2}, N_{3}, N_{4}, N_{5}\},$ it can be known by the definition of the norm on ${\mathbb R}^5$ that

Thus, lemma is proved.

Theorem 1. If $\overline{M} = \max\{N_{1}, N_{2}, N_{3}, N_{4}, N_{5}\},$ then system (2.1) has a unique positive solution $p_{*}(s).$

Proof. From the assumption $\overline{M} < 1,$ the mapping $\phi$ on $\Omega$ is a compression map. According to the compression mapping theorem, $\phi$ has a unique fixed-point $(a_{1}^{*}, a_{2}^{*}, a_{3}^{*}, a_{4}^{*}, a_{5}^{*}),$ which is the unique solution of (2.3). Therefore, (2.1) has a unique solution. It is known that the solution of system (2.3) is positive. The theorem is verified.

3.   Characteristic equation

First, linearizing system (2.1), performing a Taylor expansion of $\mu(s, c_{0}(t), J(t))$ with respect to the third variable $J(t),$ there is

It can be determined from the first formula of (1.2) that

Let $d(s, t) = p(s, t)-p_{*}(s),$ that is $p(s, t) = d(s, t)+p_{*}(s).$

Then, by substituting that into above formula, it can be determined that

Substituting the first formula of (2.1) and set

It can be determined that

Performing a Taylor expansion of $\beta(s, c_{0}(t), R(t))$ with respect to the third variable $R(t),$ there is

substituting the above formula into the fourth formula of (2.1)

That is

The characteristic equation of equilibrium is derived below.

Let

Assuming (3.3) has a formal solution $d(s, t) = D(s)e^{\lambda t},$ substituting (3.4) into (3.1) and (3.2), it can be determined that

So one gets

where

Substituting $D(s)$ into $\overline{D}_{1},$ it can be determined that

where

Substituting $D(s)$ into $\overline{D}_{2},$ it can be determined that

where

Substituting $D(s)$ into $\overline{D}(0),$ it can be determined that

where

It can be determined from $(3.7)$–$(3.9)$ that

where $(D(0), \bar{D}_1, \bar{D}_2)\neq (0, 0, 0)\Rightarrow D(s)\not \equiv 0.$

The second equation of (3.10) is multiplied by $A_{33}(\lambda)$ and used as input into the third equation,

Therefore (3.10) is transformed into

Let

therefore (3.11) is simplified to

There exists a nonzero solution to equation (3.12) if and only if the determinant of its coefficient is 0. From the above discussion, the following theorem can be obtained.

Theorem 2. Linear systems (3.1) and (3.2) have unique formal solutions $d(s, t) = D(s)e^{\lambda t},$ if and only if $\lambda$ is the root of the characteristic equation

4.   Equilibrium stability

The following theorems can be obtained from the linear stability theory.

For the reader convenience we recall that any solution of the first-order linear differential equation system $\frac{\mbox{d}X}{\mbox{d}t} = AX$ with constant coefficients, it can be expressed as $\sum\limits_{m = 0}^{l_{i}}c_{im}t^{m}e^{\lambda_{i}t}(1\leq i\leq n),$ $\lambda_{i}$ is a root of characteristic equation $\det(A-\lambda E) = 0$ ($E$ is the identity matrix), $l_{i}$ is either zero or a positive integer is determined by $\lambda_{i}.$ Moreover, if the characteristic equation $\det(A-\lambda E) = 0$ has roots with negative real parts, the solution is asymptotically stable. If the characteristic equation has roots with positive real parts, the solution is unstable.

Theorem 3. If a root of the characteristic equation $G(\lambda) = 0$ has a negative real part, the equilibrium $p_{*}(s)$ of (2.1) is asymptotically stable; if there is a root with the positive real part, the equilibrium $p_{*}(s)$ of (2.1) is unstable.

In order to obtain a specific discriminant condition, the death rate $\mu$ is considered to be independent of weighted size $J(t)$.

Corollary 1. Let $\mu$ be independent of $J(t),$ $g(0) = 1,$ then

(1) when $\Lambda < 0,$ the equilibrium of (2.1) $p_{*}(s)$ is asymptotically stable;

(2) when $\Lambda > 0,$ the equilibrium of (2.1) $p_{*}(s)$ is unstable,

where

Proof. Because $\mu$ is independent of $J(t),$ $g(0) = 1,$ then

So one gets

therefore

Let $K(\lambda) = 1-G(\lambda),$ so the characteristic equation has the following form

From the first formula of (2.3), it can be determined

thereby

It can be determined that $\lim\limits_{\lambda\rightarrow +\infty}K(\lambda) = 0$ and $K^{'}(\lambda) < 0.$

If $\Lambda > 0,$ then $K(0) = \Lambda+1 > 1.$ The equation $K(\lambda) = 1$ has at least one positive real root, therefore, the equilibrium $p_{*}(s)$ of (2.1) is unstable.

On the other hand, if $\Lambda < 0,$ suppose $K(\lambda) = 1(Re(K(\lambda)) = 1)$ has a root of $\lambda = x+iy(x\geq0),$

for $x\geq0,$ we have $e^{-x\tau(s)}\leq1, \cos(y\tau(s))\leq1,$ obviously

Contradiction with $Re(K(\lambda)) = 1,$ therefore, when $\Lambda < 0,$ all the roots of the characteristic equation have negative real parts, so that the equilibrium of system (2.1) is asymptotically stable.

5.   Algorithms and examples

Given the degree of complexity of the characteristic equation $G(\lambda) = 0,$ it is difficult to analyze the distribution of roots. We use numerical analysis to study the stability of the equilibrium

The model is a class of first-order hyperbolic equations with nonlinear variable coefficients whose solutions have local dependence on the initial values. Therefore, a windward differential format, with good approximation accuracy, is chosen. The upwind differential format for system (1.2) can be summarized as follows:

where $a$ is the number of evolution time steps and $\tau$ is the time step. The scale interval is taken as $[0, 1]$ and equally divided into $N$ segments. The scale step size $h = \frac{1}{N},$ in order to make the difference format stable, according to the freezing coefficient method and the Fourier analysis, the following two steps can be determined

Example 1. Population persistence

The initial distribution function of the population is assumed as

The parameters for setting model (1.2) are as follows:

In this example, the scale step size is $h = 0.01,$ the time step is $\tau = 0.02,$ and the time iteration step is $a = 120.$ Figure 1 is obtained from the calculation results. It can be seen that as time goes by, the density function of the population gradually moves away from the level surface, the scale distribution interval of the individual population becomes longer and longer, and the peak value increases. That is, the sustainable viability of the population is enhanced.

Example 2. Population extinction

The initial distribution function of the population is assumed as

The parameters for setting model (1.2) are as follows:

In this example, the scale step size is $h = 0.025,$ the time step is $\tau = 0.01,$ and the time iteration step is $a = 50.$ Figure 2 is obtained from the calculation results. It can be seen that over time, the density of the population gradually approaches zero, that is, the population tends to become extinct.

6.   Conclusions

Many achievements have been made in an age-structured model. In the growth process of populations, individuals of different ages can exist in the same scale stage, and scale is a better evaluation of population growth than age. The population model with scale structure proposed in this paper still has many problems to be solved, for example, the continuous variation of individual scales directly leads to mathematical models in the form of partial differential-integral equations, usually with global feedback boundary conditions, which are not a standard mathematical object and cannot be directly applied to the existing research results of mathematical theory. When the scale structure is considered, the growth rate of the population is no longer synchronized with time, which greatly increases the difficulty of analyzing the model. Through careful analysis, we obtained the existence and uniqueness condition of the equilibrium, derived the characteristic equation, gave the stability rule for judging the equilibrium, and obtained the gradual change of the equilibrium solution of the system without considering the effect of density restriction on mortality. Considering the nonlinear characteristics of the state system, it cannot be solved accurately. Therefore, an approximation algorithm is given in this paper, and the evolutionary behavior of the population is shown through two examples.

From the above stability conditions, it can be clearly seen that when the population has continuous immigration from the outside world, the stability conditions of its equilibrium are different from the situation in ^[30] that does not consider migration. This has some reasonableness in the actual situation. For example, for endangered species, due to the continuous immigration, the population can maintain a stable distribution under condition $\Lambda < 0$, and condition $\Lambda < 0$ also provides a reference for protecting endangered species. When $\mu$ is related to $J(t)$, its stability analysis is difficult and needs further study. In this paper, stability conditions are used to describe the process of population development, which plays an important role in bio-diversity conservation and species management.

Acknowledgments

The research has obtained the support from the Priority Research and Projects for Ningxia in China (2020BEG03021), the Science and Technique Research Foundation of Ningxia Institutions of Higher Education (NGY2020010), the First-class Major Foundation of Ningxia Institutions of High Education in China (NXYLXK2021A03).

Conflict of interest

The authors declare there is no conflict of interest.