On the Question of Self-Organization of Population Dynamics on Earth

Abstract

A new mathematical model based on the predator–prey interactions has been proposed. Strictly analytical solution has been found for a system of nonlinear differential equations describing accurately the dynamics of population concentration on the Earth’s land surface, taking into account the main cut-off factors influencing its spatio-temporal evolution. Using the proposed model it is shown that the population of the Earth tends to increase, and then, due to the system’s permeability, after reaching the maximum, the population size begins to decline. This observation is confirmed by rigorous mathematical and numerical calculations.

SUPPLEMENT

In order to establish the correctness of equations (1), we will use the general methodology for deriving the basic physical equations, while taking some specific features of the evolutionary development of biological systems into account. For this purpose, it is convenient to use a technique based on the law of conservation of the total power of the system, which has been successfully applied to a number of problems (see, for example, [32–37]). Similarly to problems of mechanics, we will consider the problem of the predator–prey type using the language of energy and dissipative functions, although in the context of biophysical problems they should be understood not in the literal but in a conventional sense. It is known that the kinetic energy of a hydrodynamic flow can be represented as $\varepsilon = \frac{1}{2}\int_V {\rho {{{v}}^{2}}} dV,$ where ρ is the density of a liquid and v is its velocity [38]. In our model, the notion of energy cannot be introduced for populations, but we can write the prototypical energy, as in hydrodynamics, by introducing some quadratic form by the concentration of individuals n and the cut-off factor m. As in physical problems, we denote it by E and can obviously write

$$E = \frac{1}{2}\int\limits_\Sigma {\left[ {{{C}_{1}}(m){{n}^{2}} + {{C}_{2}}(n){{m}^{2}}} \right]d\sigma ,} $$

(S1)

where C₁(n) and C₂(m) are two essentially positive functions on the corresponding arguments, Σ is the surface of Earth, and dσ is an element of the surface.

The dissipative function $\dot {Q}$ has the sense of entropy production per unit of time [39] and is essentially positive in any physical process.

In the case of biological systems, this principle should be applied with caution, considering the objective fact that changes in population size represent a chaotic process; as a consequence, for open systems, entropy production per unit of time cannot remain essentially positive all the time, since during certain intervals it would have to obey two mutually exclusive conditions: $\dot {S} > 0$ and $\dot {S} < 0$.

Therefore, for the formal description of any scalar chaos, we are entitled to represent the dissipative function as the following conditionally positive combination (see [34]):

$$\dot {S}\sim \dot {Q} = \int\limits_\Sigma {\left[ {{{\Gamma }_{1}}(n,m){{m}^{2}} - {{\Gamma }_{2}}(n,m){{n}^{2}} + {{\Gamma }_{3}}{{n}^{3}} + {{\Gamma }_{4}}n{{m}^{2}} + {{\Gamma }_{5}}m{{n}^{2}}} \right]} d\sigma ,$$

(S2)

where Γ_1, 2(n, m) are functions of the concentrations n and m, which have the meaning of generalized decrement and increment, respectively, and Γ_3, 4, 5(n, m) in the general case are also functions of n and m, but they can have any sign. Furthermore, n > 0 and m > 0 are always observed. The fact that definitions (S1) and (S2) include a considerable number of free parameters should not be too surprising, since all of them in the limiting special cases are either expressed through each other or are equal to zero.

According to the general principle of obtaining any equations of dynamics, we are entitled to proceed from the condition $\dot {E} + \dot {Q} = 0$ (as exemplified, for instance, by the Navier−Stokes equations and any equations of mechanics in which the dissipative function is essentially positive), which obviously indicates that the extremum condition is also automatically satisfied for the functions written in the form (S1) and (S2) for open systems. It is clear that as a result we should obtain equations that describe the extreme distribution of both parameters of interest (see below).

In fact, by differentiating (S1) by time, we find

$$\dot {E} = \frac{1}{2}\int\limits_\Sigma {\left[ {C_{1}^{'}\dot {m}{{n}^{2}} + 2{{C}_{1}}n\dot {n} + C_{2}^{'}\dot {n}{{m}^{2}} + 2{{C}_{2}}m\dot {m}} \right]d\sigma ,} $$

(S3)

where $C_{1}^{'} = \frac{{d{{C}_{1}}}}{{dm}}$ and $C_{2}^{'} = \frac{{d{{C}_{2}}}}{{dn}}.$

By adding (S2) and (S3), we have

$$\int\limits_\Sigma {\left[ {\frac{1}{2}C_{1}^{'}\dot {m}{{n}^{2}} + {{C}_{1}}n\dot {n} + \frac{1}{2}C_{2}^{'}\dot {n}{{m}^{2}} + {{C}_{2}}m\dot {m} + {{\Gamma }_{1}}(n,m){{m}^{2}} - {{\Gamma }_{2}}(n,m){{n}^{2}} + {{\Gamma }_{3}}{{n}^{3}} + {{\Gamma }_{4}}n{{m}^{2}} + {{\Gamma }_{5}}m{{n}^{2}}\,...} \right]d\sigma = 0.} $$

By grouping the terms with n and m, we present the integrand function as the sum of the two expressions:

$$\begin{gathered} n\left[ {{{C}_{1}}\dot {n} + \frac{1}{2}C_{1}^{'}\dot {m}n + {{\Gamma }_{3}}{{n}^{2}} - {{\Gamma }_{2}}(n,m)n + \frac{1}{2}{{\Gamma }_{4}}{{m}^{2}} + \frac{1}{2}{{\Gamma }_{5}}mn\,\,...} \right] \\ + \,\,m\left[ {{{C}_{2}}\dot {m} + {{\Gamma }_{1}}(n,m)m + \frac{1}{2}C_{2}^{'}\dot {n}m + \frac{1}{2}{{\Gamma }_{4}}nm + \frac{1}{2}{{\Gamma }_{5}}{{n}^{2}}\,...} \right] = 0. \\ \end{gathered} $$

Using the independence of n and m, we separately set either of them equal to zero and find the system of equations we are interested in:

$$\left\{ \begin{gathered} {{C}_{1}}\dot {n} + \frac{1}{2}C_{1}^{'}\dot {m}n + {{\Gamma }_{3}}{{n}^{2}} - {{\Gamma }_{2}}n + \frac{1}{2}{{\Gamma }_{4}}{{m}^{2}} + \frac{1}{2}{{\Gamma }_{5}}mn\,\,... = 0, \hfill \\ {{C}_{2}}\dot {m} + \frac{1}{2}C_{2}^{'}\dot {n}m + {{\Gamma }_{1}}m + \frac{1}{2}{{\Gamma }_{4}}nm + \frac{1}{2}{{\Gamma }_{5}}{{n}^{2}}\,... = 0. \hfill \\ \end{gathered} \right.$$

(S4)

Let us rewrite these equations in the following form:

$$\left\{ \begin{gathered} \dot {n} + \frac{{C_{1}^{'}}}{{2{{C}_{1}}}}\dot {m}n = - \frac{1}{{{{C}_{1}}}}\left( {{{\Gamma }_{3}}{{n}^{2}} - {{\Gamma }_{2}}n + \frac{1}{2}{{\Gamma }_{4}}{{m}^{2}} + \frac{1}{2}{{\Gamma }_{5}}mn\,\,...} \right) = {{\varphi }_{1}}(n,m), \hfill \\ \dot {m} + \frac{1}{{2{{C}_{2}}}}C_{2}^{'}\dot {n}m = - \frac{1}{{{{C}_{2}}}}\left( {{{\Gamma }_{1}}m + \frac{1}{2}{{\Gamma }_{4}}nm + \frac{1}{2}{{\Gamma }_{5}}{{n}^{2}}\,...} \right) = {{\varphi }_{2}}(n,m). \hfill \\ \end{gathered} \right.$$

(S5)

By solving this system, which is linear with respect to the derivatives $\dot {n}$ and $\dot {m}$, we obtain

$$\left\{ \begin{gathered} \dot {m} + \frac{1}{{{{C}_{2}} - \frac{{C_{1}^{'}C_{2}^{'}}}{{4{{C}_{1}}}}mn}}\left[ {{{\Gamma }_{1}}m + \frac{1}{2}\left( {\frac{{C_{1}^{'}}}{{{{C}_{1}}}}{{\Gamma }_{2}} + {{\Gamma }_{4}}} \right)mn + \frac{{{{\Gamma }_{5}}}}{2}{{n}^{2}} - \frac{{C_{2}^{'}}}{{4{{C}_{1}}}}{{\Gamma }_{4}}{{m}^{3}} - \frac{{C_{2}^{'}}}{{2{{C}_{1}}}}{{\Gamma }_{3}}m{{n}^{2}} - \frac{{C_{2}^{'}}}{{4{{C}_{1}}}}{{\Gamma }_{5}}n{{m}^{2}}\,...} \right] = 0, \hfill \\ \dot {n} + \frac{1}{{{{C}_{1}} - \frac{{C_{1}^{'}C_{2}^{'}}}{{4{{C}_{2}}}}mn}}\left[ { - {{\Gamma }_{2}}n + {{\Gamma }_{3}}{{n}^{2}} + \frac{{{{\Gamma }_{4}}}}{2}{{m}^{2}} + \frac{1}{2}\left( {{{\Gamma }_{5}} - \frac{{C_{1}^{'}}}{{{{C}_{2}}}}{{\Gamma }_{1}}} \right)nm - \frac{{C_{1}^{'}{{\Gamma }_{4}}}}{{4{{C}_{2}}}}m{{n}^{2}} - \frac{{C_{1}^{'}{{\Gamma }_{5}}}}{{4{{C}_{2}}}}{{n}^{3}}\,...} \right] = 0. \hfill \\ \end{gathered} \right.$$

(S6)

From the general system of equations (S6), it can be seen immediately that for m = 0 and Γ₅ = 0 we obtain the Verhulst equation: $\dot {n} = \alpha n - \beta {{n}^{2}}$. For other proportions, we get a generalized system of Lotka–Volterra equations.

Assuming in (S6) that C₂(n, m) = const and omitting all cubic terms by n and m, as well as introducing explicit dependences on n and m for the sake of convenience, we have

$$\left\{ \begin{gathered} \dot {m} + \frac{1}{{{{C}_{2}}}}\left( {{{\Gamma }_{1}}(n,m)m + \frac{{{{\Gamma }_{4}}(n,m)}}{2}mn + \frac{{{{\Gamma }_{5}}(n,m)}}{2}{{n}^{2}}} \right) = 0, \hfill \\ \dot {n} + \frac{1}{{{{C}_{1}}(n,m)}}\left[ { - {{\Gamma }_{2}}(n,m)n + {{\Gamma }_{3}}(n,m){{n}^{2}} + \frac{{{{\Gamma }_{4}}(n,m)}}{2}{{m}^{2}} + \frac{1}{2}\left( {{{\Gamma }_{5}}(n,m) - \frac{{C_{1}^{'}}}{{{{C}_{2}}}}{{\Gamma }_{1}}(n,m)} \right)nm} \right] = 0. \hfill \\ \end{gathered} \right.$$

(S7)

The case considered in our study corresponds to the conditions C₁(n, m) = C₁m, Γ₁(n, m) = Γ₁m, Γ₂(n, m) = Γ₂m, Γ₃(n, m) = Γ₃ = const, Γ₄ = const, Γ₅ = 0.

As a result, it follows from equations (S7) that

$$\left\{ \begin{gathered} \dot {n} + \left[ { - \frac{{{{\Gamma }_{2}}}}{{{{C}_{1}}}}n + \frac{{{{\Gamma }_{3}}}}{{{{C}_{1}}m}}{{n}^{2}} + \frac{{{{\Gamma }_{4}}}}{{2{{C}_{1}}}}m - \frac{{{{C}_{1}}}}{{2{{C}_{2}}}}{{\Gamma }_{1}}(m)nm} \right] = 0, \hfill \\ \dot {m} + \frac{1}{{{{C}_{2}}}}\left( {{{\Gamma }_{1}}(m)m + \frac{{{{\Gamma }_{4}}}}{2}mn} \right) = 0. \hfill \\ \end{gathered} \right.$$

(20)

To relate to the parameters of equations (1), we assume that the coefficient Γ₁ is negative and

$$\alpha = \frac{{{{\Gamma }_{2}}}}{{{{C}_{1}}}},\,\,\,\,\beta (t) = \frac{{{{\Gamma }_{3}}(t)}}{{{{C}_{1}}}},\,\,\,\,\gamma = \frac{{{{\Gamma }_{4}}}}{{2{{C}_{2}}}},\,\,\,\,\frac{{{{\Gamma }_{1}}(m)}}{{{{C}_{2}}}} = - \lambda \varepsilon {{m}^{\lambda }}.$$

As a result, from (S8), it follows that:

$$\left\{ \begin{gathered} \dot {n} = \alpha n - \frac{{\beta (t)}}{m}{{n}^{2}} - \gamma \frac{{{{C}_{2}}}}{{{{C}_{1}}}}m - \frac{{\lambda \varepsilon }}{2}n{{m}^{{1 + \lambda }}}, \hfill \\ \dot {m} = - \gamma mn + \lambda \varepsilon {{m}^{{1 + \lambda }}}. \hfill \\ \end{gathered} \right.$$

(S9)

If we compare equations (S9) to equations (1), we see that the right-hand side of the upper equation (S9) contains two additional terms that are absent from equations (1). It would seem that this point could be a strong reason to criticize the entire work. However, we should not jump to conclusions. As it can be seen from (S9), these terms have the correct signs; if we assume here C₁ = C₂ and utilize the numerical values of the parameters provided above (α = 10^–2 kg/(km² day), β = 10^–5 kg/(km² day), γ = 10^–6 kg/(km² day), λ = 10^–2, ε = 10^–2 (SI, km, day)), system (S9) acquires the following explicit form:

$$\left\{ \begin{gathered} \dot {x} = - {{10}^{{ - 6}}}xy + {{10}^{{ - 4}}}{{x}^{{1.01}}}, \hfill \\ \dot {y} = {{10}^{{ - 2}}}y - \frac{{{{{10}}^{{ - 5}}}}}{x}{{y}^{2}} - {{10}^{{ - 6}}}x - 0.5 \times {{10}^{{ - 4}}}y{{x}^{{1.01}}}. \hfill \\ \end{gathered} \right.$$

(S10)

The results of numerical integration of equations (S10) are shown in Fig. 3. Comparison to Figs. 1 and 2 reveals their striking similarity.

Fig. 3.

The evolution of the population size x(y) according to equations (S10) under assumption that the probability of action of the cut-off factor β is constant. Here, simulation was performed for the initial conditions x(0) = 1, y(0) = 1 until the moment t = 8000. Neglecting the last two terms in the bottom equation of system (S10) does not have any effect on the behavior of dependence x(y) and results in a shortened system of equations (1), which was investigated in the present work.

Thus, we have obtained a correct and well-founded (although shortened in comparison to equations (S9)) system of equations (1) that models the evolution of the human population size taking the factor that cuts off its development into account.