Computer research of a nonlinear model of population dynamics taking into account trophic chains*

Abstract

We develop an approach to the computer study of a threedimensional nonlinear model of population dynamics, taking into account trophic chains and competition in prey populations. This approach is based on the application of stability research methods, numerical methods for solving ordinary differential equations, and modern methods of parametric optimization. A qualitative study of the system is carried out, the states of equilibrium are found, of the species population dynamics graphs and corresponding phase portraits are constructed. An assessment of the parameters influence on the model stability and on the permanent coexistence of populations is presented. A software package in the Python language is used as a research implementation tool. The problem of the optimal control of the model with phase constraints is formulated. A control quality criterion is proposed and a generalized algorithm for solving the optimal control problem is developed. The results obtained can be used in problems of modeling and forecasting multidimensional ecological systems, as well as optimal control problems.

Keywords

• Computer Modeling
• Nonlinear Model
• Stability
• Phase Portrait
• Optimal Control

-

Currently, the study of continuous nonlinear dynamical models described by multidimensional differential equations is an important aspect of many fundamental and applied scientific fields. The development of computer modeling tools in the field of qualitative and numerical analysis of dynamical behavior systems gives researchers a number of new opportunities in problems trajectories search with implementation of *Copyright © 2021 for this paper by its authors. Use permitted under Creative Commons License numerical methods. New opportunities in the field of construction and analysis of the population dynamics models based on trophic chains stability should be noted.

The population is constantly affected by various type factors: biotic and abiotic. The mathematical model allows us to identify the factors that have the maximum impact on the dynamics of a particular population. Population dynamics models are widely used in the study of interrelated ecological systems, as well as in socioeconomic interactions of society.

The organisms that make up populations are in complex relationships with each other, as a result of which there is a positive or negative influence of some species on others. This type of interaction is usually classified as [1]. Among all the diversity of relations between the two populations, the main ones are considered to be: predatorprey (the number of prey species grows more slowly than the number of predator species), competition (the number of each species grows more slowly in the presence of the other), and symbiosis(species contribute to each other's population growth). Classic predator-prey models and competitor-competitor is the subject of many papers (for example, [2–4]). In [5-6], various models are introduced in the presence of a beneficial effect of species living in biocenoses on each other.

When studying models of population dynamics, an urgent problem is the study of the nonlinear models stability [1–6]. The Lyapunov function method is one of the most important methods for studying stability [7-8]. In [9-10], a study of the population systems stability is carried out based on the construction of stochastic selfconsistent models. When studying models of ecological systems, the use of applied mathematical packages and programming languages is relevant [11-12]. A number of control problems for population dynamics models are studied in [13-14] and in other papers. In particular, some optimal control problems for distributed population dynamics models are considered in [13].

In this paper, we consider a three-dimensional nonlinear model of the population dynamics of interconnected communities, taking into account competition and the predator-prey interaction. To study the stability of this model, a computer study is conducted. Model parameters are estimated and phase portraits of the system and the population dynamics graphs are constructed using the Jupyter software package. The obtained effects are analyzed. 2

Models and methods

We study a nonlinear model in which there are two competing types. Competitor species are also prey species that interact with the predator population. It is assumed that the presence of a predator-prey interaction can have a positive effect on competitors, and both competing species can coexist. We denote as xi the preys dy population density, as y the predators population density, as y = , xi = dxi time dt dt derivatives, i = 1, 2. According to the ecological sense, we have y0  0, xi 0  0, i = 1, 2 . The population dynamics model is an autonomous system described by three ordinary nonlinear differential equations of the form [15]: x1  x1 (a1  11x1  12 x2  b1 y), 2  x2 (a2  21x1  22 x2  b2 y), x y  y(c  d1 x1  d 2 x2 y), where x1 is the population density of the first competitor, x2 is the population density of the second competitor, y is the density population predator, ai is the reproduction coefficient of population competitor in the absence of predator, bi is the coefficient of consumption by the predator population of the prey population, di is coefficient processing consumed predator biomass the prey in private biomass, i = 1, 2, εij are coefficients of intraspecific competition for i = j = 1 and i = j = 2, εij are coefficients of interspecific competition for i not equal to j, с is natural mortality of the predator, γ is intraspecific competition of the predator. According to the ecological sense, the coefficients have the form ij  0,   0, ai > 0, bi > 0, di > 0, с > 0, i, j = 1, 2. Next, we introduce the following notations:   b2 d 211  b1d1 22  b2 d112  b1d 2 21 ,    11 22  12 21 , D     x1*  x2*  (a1b2  a2 b1 )d 2  (b1 22  b2 12 )c  (a1 22  a2 12 ) ,

D (a2 b1  a1b2 )d1  (b2 11  b1 21 )c  (a2 11  a1 21 ) ,

D y *  (a1 22  a2 12 )d1  (a2 11  a1 21 )d 2  c .

D

The equilibrium states of the model ( 1 ) are obtained in a general form. These equilibrium states have the form: Е0 (0,0,0), Е1  a111 ,0,0  , Е2  0, a222 ,0 , Е3  0,0, c ,  Е4  0,   22  b2 d 2 a2   b2c , a2 d 2  c 22 , Е5  a111  bb11dc1 ,0, a1d1  c22 ,

 22  b2 d 2  11  b1d1  Е6  a1 22  a2 12 , a2 11  a1 21 ,0, Е7 x1* , x2* , y* 

    

State of equilibrium E7 is an internal equilibrium state for which we assume that the positivity condition is satisfied.

Next, we consider the following conditions for the model ( 1 ).

A1: a1 22  a2 12 , A2: a2 11  a1 21 , A3: a1 22  a2 12, A4: a2 11  a1 21 , A5: D> 0, ( 1 ) ( 2 ) ( 3 ) ( 4 ) A7: { 11 x1*   22 x2*  y *}{x1* x2*  ( 11  b1d1 ) x1* y *  ( 22  b2 d 2 ) x2* y*}  x1* x2* y* D , A8: { 11 x1*   22 x2*  y *}{x1* x2*  ( 11  b1d1 ) x1* y *  ( 22  b2 d 2 ) x2* y*}  x1* x2* y* D . The following statements are true.

1. Internal positive state of equilibrium E7 systems ( 1 ) exists if and only if it holds ( A1  A2 )  ( A3  A4 ) .

2. If it holds ( A1  A2 ) , then system ( 1 ) is asymptotically stable. 3. If it holds ( A3  A4 ) , then the internal equilibrium state E7 is a saddle point. 4. If it holds ( A5  A7 ) , then the internal equilibrium state E7 of system ( 1 ) is asymptotically stable.

5. If it holds ( A6  A8 ) , then the internal equilibrium state E7 of system ( 1 ) is unstable.

6. If it holds A5  ( A1  A2 ) , then the populations in system ( 1 ) permanently coexist.

In [15], theoretical studies of stability and permanent coexistence are carried out using conditions A1–A8 for the model ( 1 ). The above statements 1–6 are modifications of the results [15].

This paper is a continuation of the research conducted in [15–20]. In particular, in [19] for a special case of the model ( 1 ) the phase portraits are obtained and stability is studied. In [20] authors consider the population dynamics model «predator – two preys». In particular, in this article a deterministic stability of limit cycles of this three dimensional model in a period doubling bifurcations zone at the transition from an order to chaos is investigated.

In this paper, we study model number ( 1 ) by means of computational experiments on the basis of numerical solution of ordinary differential equations. We use methods of stability theory and qualitative theory of differential equations and modified Runge–Kutta methods 4 orders. In addition, we generalize model ( 1 ) to controlled case and consider the optimal control problem. 3

Computer experiments results

For the model ( 1 ) numerical experiments are carried out using a developed software package based on Python in the Jupyter development environment [21].

We consider the following set of parameters: (x1, x2, y) = (1.5, 2, 4), (a1, a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (16, 8, 2.5, 8, 4, 3.5, 1, 1, 1, 0.3, 0.3, 0). For this set of parameters, we obtain the following approximate equilibria states: (0, 2.5, 5.5), ( 0, 8, 0 ), ( 2, 0, 0 ), (0.33, 2.17, 4.67), and trivial equilibrium state (0, 0, 0). Let's denote as E( 1 ) point (x1*, x2* , y*)  (0.33, 2.17, 4.67) . Figure 1 shows the phase portrait of the 7 model ( 1 ) for a given set of parameters in the neighbourhood of a point E7( 1 ) .

Fig. 1. Phase portrait of the model ( 1 ) in the neighbourhood E7( 1 ) in three-dimensional space.

Figure 2 shows the phase portrait of the model ( 1 ) with respect to projections (х1, х2) in the neighbourhood of the equilibrium state E7( 1 ) .

Further, we consider the second set of parameters: (x1, x2, y) = ( 2, 1, 6 ), (x1, x2, y) = (0.4, 4, 5), (a1, a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (16, 8, 2.5, 8, 4, 4.125, 1, 1, 1, 1, 1, 0). The equilibrium states of the system ( 1 ) have the form: (0, 0, 0), (0, 2.5, 5.5), Е2( 0, 8, 0 ), Е3(1.88, 0.24, 0), Е4( 2, 0, 0 ), Е5(0.57, 1.9, 3.7). Denote by E ( 2 ) and E ( 2 ) points ( 2, 0, 0 ) and (0.57, 1.9, 3.7) respectively. Figure 3 shows the 1 7 phase portrait of the model ( 1 ) in the neighbourhood of equilibrium states E ( 2 ) and 1 E ( 2 ) in three-dimensional space. Figure 4 shows the phase portrait of the model ( 1 ) 7 relatively to the projections (х1, х2) in the neighbourhood of equilibrium states E ( 2 ) 1 and E7( 2 ) . It should be noted that two semi-orbits are observed in Fig. 3, 4. One tends states: E ( 2 ) under initial conditions (x1, x2, y) = ( 2, 1, 6 ) and

1 conditions (x1, x2, y) = (0.4, 4, 5). to the limit cycle, the other one tends to the asymptotically stable equilibrium point at the boundary. This fact is consistent with the results of [15].

According to figures 3 and 4, there are two asymptotically stable equilibrium E ( 2 ) under initial 7

Fig. 4. Phase portrait of the model ( 1 ) relatively to projections (х1, х2) in the

neighbourhood of the equilibrium state E1( 2 ) и E7( 2 ) .

Note that the choice of the parameters for experiments is consistent with the parameters considered in [15], and for these sets of parameters, the conditions of stability and permanent coexistence are met. Developed software package based on Python 3 in the Jupyter system allows you to specify the results of constructing phase portraits of the model ( 1 ).

Figures 5 and 6 show the dynamics of the population density under the two sets of initial conditions indicated above. According to the computational experiments of the model ( 1 ), periodic fluctuations in the population of species are established. The graphs confirm that competitive prey (and predator) populations are able to permanently coexist.

Figures 6 and 7 show the dynamics competing species and predator populations under different initial conditions. In this case, the populations in the system ( 1 ) with the specified set of parameters permanently coexist only for (x1, x2, y) = (0.4, 4, 5). Fig. 6. Dynamics of populations x1, x2, y under initial conditions: (x1, x2, y) = (0.4, 4, 5), (a1, a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (16, 8, 0.83, 8, 4, 4.125, 1, 1, 1, 1, 1, 0).

Fig. 7. Dynamics of populations x1, x2, y under initial conditions: (x1, x2, y) = ( 2, 1, 6 ), (a1, a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (16, 8, 0.83, 8, 4, 4.125, 1, 1, 1, 1, 1, 0).

Further, let us consider the influence of the system parameters correction ( 1 ) on stability. In particular, the following parameter, as the speed of competitors reproduction in the absence of the predator has essential meaning for the dynamics of the systems with trophic chains. We consider a set of parameters for a stable system ( 1 ) with permanent coexistence (similar to figure 5). Taking into account deviations a01 = a02 = 0.2 we get the change in population density shown in figures 8 and 9. Fig 8. Dynamics of populations x1, x2, y under initial conditions: (x1, x2, y) = (0.4, 4, 5), (a1, a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (15.8, 8, 0.83, 8, 4, 3.5, 1, 1, 1, 0.33, 0.33, 0).

Thus, according to the first «perturbed» set of parameters the trajectories of the system ( 1 ) with asymptotically stable equilibrium and permanent coexistence of populations are obtained. In this case, the conditions ( A1  A2 ) , ( A5  A7 ) , A5  ( A1  A2 ) are hold.

Fig 9. Dynamics of populations x1, x2, y under initial conditions: (x1, x2, y) = ( 0,4, 4, 5 ), (a1, a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (16.2, 8, 0.83, 8, 4, 3,5, 1, 1, 1, 0.33, 0.33, 0).

In addition, according to the second «perturbed» set of parameter the trajectories of the system ( 1 ) with permanent coexistence of populations are obtained. In this case, the conditions A5  ( A1  A2 ) are hold.

Computational experiments show, that the specific gravity of speed consumption by the population the conversion rate of the prey's biomass consumed by the predator into its own biomass significantly affect the character of the system ( 1 ) stability.

Optimal control in a three-dimensional population model

For the three-dimensional model ( 1 ), we formulate the optimal control problem. The dynamics of the controlled model is defined by a system of differential equations x1  x1 (a1   11x1   12 x2  b1 y)  u1x1 , x2  x2 (a2   21x1   22 x2  b2 y)  u 2 x2 , y  y(c  d1x1  d 2 x2  y)  u3 y, where ui = ui(t) are control functions, xi  0, y  0, i  1,2. The meaning of the parameters appearing in ( 5 ) is explained in section 2.

Constraints for the model ( 5 ) are set as x1 (0)  x10 , x2 (0)  x20 , x3 (0)  x30 , x1 (T )  x11, x2 (T )  x21, x3 (T )  x31, t  0, T , ( 5 ) ( 6 ) ( 7 ) 0  u1  u11, 0  u2  u21, 0  u3  u31, t  0, T .

In relation to problem ( 5 )–( 7 ) we consider the functional to be minimized as

T 3

J (u)  0 i1 ki ui (t)dt. ( 8 )

The control quality criterion ( 8 ) corresponds to minimizing losses from population size regulation, with ki is value of one population relatively to another.

For the model ( 1 ) the optimal control problem can be formulated as follows: to find the minimum of the functional ( 8 ) under conditions ( 6 ), ( 7 ).

Methods of control theory and artificial intelligence allow us to construct control laws u1, u2, u3. In particular, in some cases, the use of PID controllers or controllers with the use of sliding mode is effective. We suggest using methods based on machine learning and building controllers using artificial intelligence. In particular, it is possible to use machine learning in combination with controllers based on fuzzy logic or artificial neural networks.

In order to search for functions u1, u2, u3 we can consider the construction of a parametric control model. For example, a polynomial control of the form ui (t )  RT , R  (ri1 , ri2 ,..., rin )T , T  (t 0 , t1 ,..., t m ), i  1, 2, 3, used in [22] for a model with migration flows. In this case, the model parameters are coefficients ri1,…, rin the polynomial functions. Global parametric optimization methods, in particular differential evolution, are used to calculate the parameters [22].

Similar to polynomial control, for artificial intelligence-based controllers, parameters are formal characteristics that determine their internal structure. We propose a generalized algorithm for solving the optimal control problem based on machine learning. This algorithm consists of the following steps. ─ To construct a formal parametric control model for the system ( 5 ). ─ To adjust the parameters of the control model using global parametric optimization. ─ To evaluate the control quality criterion. If the required characteristics are achieved, proceed to the next step. Otherwise, go back to step 2. ─ To construct the optimal trajectory.

This algorithm can be used in the framework of solving the optimal control problem formulated in this paper for the model ( 5 ).

5 Discussion A computer study of a nonlinear model of interaction between two competing prey individuals and a predator population made it possible to study the stability of the proposed model under various sets of variables and initial conditions. With the help of developed computer programs, graphs of population dynamics are constructed. With the considered sets of parameters, we obtained the influence estimation of the predator species on the result of prey competition. In a number of identified cases, the presence of trophic chains (predation) has a positive effect on the result of competition and contributes to the coexistence of species. We have formulated a new optimal control problem and proposed a control quality criterion.

The developed algorithm allows to construct a parametric control model, perform correction of the model parameters, evaluate the control quality criterion and obtain optimal trajectories.

The software package developed in Python using the Jupyter development environment has high flexibility and scalability. In the future, it is planned to develop this software package in order to adapt to a wide class of mathematical models with logic controllers.

6 Conclusion

In this paper, we develop an approach to computer research of a nonlinear threedimensional model of population dynamics with trophic chains. This approach is based on the application of stability research methods, numerical methods for solving ordinary differential equations, and modern methods of parametric optimization.

The software package developed in Python using the Jupyter development environment has shown high efficiency for computer research of a multidimensional nonlinear model. The obtained results can be used in computer modeling of ecological and socio-economic systems. A statement of the optimal control problem in the model ( 2 ) is proposed and an algorithm for solving this problem is given. To solve this problem, it is proposed to use numerical optimization methods and intelligent algorithms for symbolic calculations. As a perspective for further research, a computer study of the proposed model with partial control should be noted. 21. Fuhrer, C., Solem, J.E., Verdier, O.: Scientific Computing with Python, 3, Packt

Publishing (2016). 22. Demidova, A.V., Druzhinina, O.V., Masina, O.N., Petrov, A.A.: Computer research of the controlled models with migration flows. Kulyabov, D.S., Samouylov, K.E., Sevastianov, L.A. (eds.) ITTMM-2020, CEUR Workshop Proceedings, 2639, 117–129, CEUR-WS, RWTH Aachen University (2020).

References

1. Bazykin , A.D. : Nonlinear dynamics of interacting populations . Institute of computer research, Moscow-Izhevsk ( 2003 ).

2. Murray , J.D.: Mathematical Biology: I. An Introduction . Springer-Verlag, New York ( 2007 ).

3. Volterra , V.: Mathematical theory of struggle for existence . Nauka , Moscow ( 1976 ).

4. Svirezhev , Yu.M. , Logofet , D.O. : Stability of biological communities . Nauka , Moscow ( 1978 ).

5. Freedman , H.I. , Rai , B. : Uniform Persistence and Global Stability in Models Involving Mutualism Competitor-Competitor-Mutualist Systems . Indian J. Math. 30 , 175 - 186 ( 1988 ).

6. Druzhinina , O.V. , Masina , O.N. , Shcherbakov , A.V. : Structure and Qualitative Analysis of Mathematical Models of Population Dynamics in the Presence of Mutualism . Nonlinear World, 14 ( 6 ), 32 - 42 ( 2016 ).

7. Shestakov , A.A. : Generalized direct method for systems with distributed parameters . URSS , Moscow ( 2007 ).

8. Druzhinina , O.V. , Masina , O.N. : Methods for studying the stability and controllability of fuzzy and stochastic dynamical systems . Dorodnitsyn Computing Center of RAS , Moscow ( 2009 ).

9. Demidova , A.V. , Druzhinina , O.V. , Masina , O.N. : Study of the Stability of the Population Dynamics Model Based on the Construction of Stochastic Self-consistent Models and the Reduction Principle . Bulletin of Peoples' Friendship University of Russia. Series Mathematics. Information Sciences. Physics , 3 , 18 - 29 ( 2015 ).

10. Demidova , А.V., Druzhinina , О.V. , Masina , O.N. , Tarova , E.D. Computer research of nonlinear stochastic models with migration flows . In: Kulyabov, D.S. , Samouylov , K.E. , Sevastianov , L.A . (eds.) ITTMM-2019 , CEUR Workshop Proceedings , 2407 , 26 - 37 , CEUR-WS, RWTH Aachen University ( 2019 ).

11. Lamy , R.: Instant SymPy Starter. Packt Publishing ( 2013 ).

12. Oliphant , T.E.: Guide to NumPy. 2nd edn . CreateSpace Independent Publishing Platform, USA ( 2015 ).

13. Moskalenko , A.I. : Methods of nonlinear mappings in optimal control. Theory and applications to models of natural systems . Nauka, Novosibirsk ( 1983 ).

14. Kuzenkov , O.A. , Kuzenkova , G.V. : Optimal Control of Self-Reproduction Systems . Journal of Computer and Systems Sciences International , 51 , 500 - 511 ( 2012 ).

15. Hutson , V. , Vickers , G.T.: A Criterion for Permanent Coexistence of Species, with an Application to a Two-Prey One-Predator System . Math. Biosci., 63 , 253 - 269 ( 1983 ).

16. Fujii , K. : Complexity-Stability Relationship of Two-Prey One-Predator Species System Model: Local and Global Stability . J. Theor. Biol ., 69 , 613 - 623 ( 1977 ).

17. Hsu , S.B. : Predator-mediated Coexistence and Extinction . Math. Biosci., 54 , 231 - 248 ( 1981 ).

18. Vance , R.R. : Predation and Resource Partitioning in One Predator-Two Prey Mode Communities . Amer. Natur ., 112 , 797 - 813 ( 1978 ).

19. Tarova , E.D.: Computer Modeling and Stability Research of Dynamics Models Interconnected Communities. Continuum. Mathematics. Computer science. Education , 3 ( 15 ), 87 - 94 ( 2019 ).

20. Bashkirceva , I.A. , Karpenko , L.V. , Ryashko , L.B. : Stochastic Sensitivity of Limit Cycles for "Predator-Two Preys" Model . Izvestiya VUZ. Applied Nonlinear Dynamics , 18 ( 6 ), 42 - 64 ( 2010 ).