=Paper= {{Paper |id=Vol-2843/paper5 |storemode=property |title=Computer research of a nonlinear model of population dynamics taking into account trophic chains (paper) |pdfUrl=https://ceur-ws.org/Vol-2843/paper005.pdf |volume=Vol-2843 |authors=Alexey A. Petrov,Alexander V. Shcherbakov,Olga V. Druzhinina,Olga N. Masina }} ==Computer research of a nonlinear model of population dynamics taking into account trophic chains (paper)== https://ceur-ws.org/Vol-2843/paper005.pdf
    Computer research of a nonlinear model of population
       dynamics taking into account trophic chains*

        Alexey A. Petrov1[0000-0002-3011-9278], Alexander V. Shcherbakov1[0000-0002-4752-8422],
          Olga V. Druzhinina2,3[0000-0002-9242-9730] and Olga N. Masina1[0000-0002-0934-7217]
        1
     Bunin Yelets State University, 28, Kommunarov st., Yelets, 399770, Russian Federation
    2
    Federal Research Center “Computer Science and Control” of Russian Academy of Science,
               44, building 2, Vavilov st., Moscow, 119333, Russian Federation,
       3
        V.A. Trapeznikov Institute of Control Sciences of Russian Academy of Sciences,
                  65, Profsoyuznaya st., Moscow, 117997, Russian Federation
                                     ovdruzh@mail.ru



            Abstract. We develop an approach to the computer study of a three-
            dimensional 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.


1           Introduction

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
Attribution 4.0 International (CC BY 4.0).
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 socio-
economic 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: predator-
prey (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 self-
consistent 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        dx
population density, as y the predators population density, as y =      , xi = i time
                                                                     dt         dt
derivatives,        i = 1, 2.     According      to    the   ecological     sense,   we
have y 0  0, x i 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   11 x1   12 x 2  b1 y ),
                              x 2  x 2 (a 2   21 x1   22 x 2  b2 y),                        (1)
                               y  y(c  d1 x1  d 2 x 2  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  b1 d1 22  b2 d 112  b1 d 2  21    11 22   12  21 D       (2)
                                                            ,                        ,

                    ( a1b 2  a 2 b1 ) d 2  (b1 22  b2  12 )c  (a1 22  a 2  12 )
              x1*                                                                          ,
                                                      D
                    (a b  a1b2 )d 1  (b2  11  b1 21 )c  (a 2  11  a1 21 )
              x 2*  2 1                                                                  ,        (3)
                                                     D
                    ( a   a 2  12 )d 1  (a 2  11  a1 21 ) d 2  c
              y *  1 22                                                    .
                                                D
  The equilibrium states of the model (1) are obtained in a general form. These
equilibrium states have the form:
                              a                a                    c
              Е0 (0,0,0), Е1  1 ,0,0  , Е2  0, 2 ,0 , Е3  0,0,     ,
                               11              22                   
                a   b2 c a 2 d 2  c 22            a   b1c       a d  c 22 
           Е4  0, 2          ,               , Е5  1            ,0, 1 1         ,
                 22  b2 d 2  22  b2 d 2         11  b1 d1    11  b1d1             (4)

                        a   a 2  12 a 2  11  a1 21 
                   Е6  1 22
                              
                                       ,
                                                 
                                                                         
                                                         ,0 , Е7 x1* , x 2* , 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: a 2  11  a1 21 ,
   A3: a1 22  a212 ,
   A4: a2 11  a1 21 ,
   A5: D> 0,
    A6: D<0,
               *          *     *     * *                       * *                       * *      * * *
    A7: { 11 x1   22 x 2  y }{ x1 x 2  ( 11  b1 d1 ) x1 y  ( 22  b2 d 2 ) x 2 y }  x1 x 2 y D ,
    A8: { 11 x1*   22 x 2*  y * }{ x1* x *2  ( 11  b1 d1 ) x1* y *  ( 22  b2 d 2 ) x 2* y * }  x1* x 2* 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
E7(1) point ( x1* , x2* , y * )  (0.33, 2.17, 4.67) . Figure 1 shows the phase portrait of the
                                                                                                            (1)
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 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) .




        Fig. 2. Phase portrait of the model (1) relatively to projections (х1, х2) in the
                        neighbourhood of the equilibrium state E 7(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
E1(2) and E7(2) points (2, 0, 0) and (0.57, 1.9, 3.7) respectively. Figure 3 shows the
                                                                                         ( 2)
phase portrait of the model (1) in the neighbourhood of equilibrium states E1                   and
 E7(2) in three-dimensional space. Figure 4 shows the phase portrait of the model (1)
                                                                                   ( 2)
relatively to the projections (х1, х2) in the neighbourhood of equilibrium states E1
        ( 2)
and E7 . It should be noted that two semi-orbits are observed in Fig. 3, 4. One tends
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
          ( 2)                                                   ( 2)
states: E1 under initial conditions (x1, x2, y) = (2, 1, 6) and E7 under initial
conditions (x1, x2, y) = (0.4, 4, 5).




     Fig. 3. Phase portrait of the model (1) in the neighbourhood of equilibrium states and
                                  E 7( 2 ) in three-dimensional space.
          Fig. 4. Phase portrait of the model (1) relatively to projections (х1, х2) in the
                      neighbourhood of the equilibrium state E1( 2 ) и E 7( 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.




 Fig. 5. Dynamics of populations x1, x2, y under initial conditions: (x1, x2, y) = (1.5, 2, 4), (a1,
      a2, c, ε11, ε12, ε21, ε22, b1, b2, d1, d2, γ) = (16, 8, 0.83, 8, 4, 3.5, 1, 1, 1, 0.33, 0.33, 0).
   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.
4       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   11 x1   12 x 2  b1 y )  u1 x1 ,
                         x 2  x 2 (a 2   21 x1   22 x 2  b2 y )  u 2 x 2 ,                          (5)
                         y  y (c  d1 x1  d 2 x 2  y )  u 3 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 , x 2 (0)  x 20 , x3 (0)  x30 , x1 (T )  x11 , x 2 (T )  x 21 , x3 (T )  x31 ,
                                                                                                            (6)
     t  0, T ,

                       0  u1  u11 , 0  u 2  u 21 , 0  u 3  u31 , t  0, T .                         (7)
    In relation to problem (5)–(7) we consider the functional to be minimized as
                                                   T 3
                                        J (u )  
                                                   0     k i u i (t )dt.                                   (8)
                                                     i 1
   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
              u i (t )  RT ,   R  ( ri1 , ri 2 ,..., rin ) T , T  (t 0 , t 1 ,..., 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 three-
dimensional 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.
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).
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).