26 UDC 004.4 Computer research of nonlinear stochastic models with migration flows Anastasia V. Demidova* , Olga V. Druzhinina† , Olga N. Masina‡ , Ekaterina D. Tarova‡ * Department of Applied Probability and Informatics Peoples’ Friendship University of Russia Miklukho-Maklaya str. 6, Moscow, 117198, Russian Federation † Federal Research Center “Computer Science and Control” of Russian Academy of Sciences Vavilov str. 44, building 2, Moscow, 119333, Russian Federation ‡ Bunin Yelets State University Communards str. 28, Yelets, 399770, Russian Federation Email: demidova_av@rudn.university, ovdruzh@mail.ru, olga121@inbox.ru, katerina.tarova@yandex.ru The problems of the design methods extension and computer research of nondeterministic finite-dimensional population models describing migration flows are studied. The significant difficulties arise in the construction of high dimension dynamic models in the course of analytical research. Computer research allows not only to obtain the results of numerical experiments to search for trajectories and estimate the parameters of deterministic models, but also to reveal the effects caused by stochasticization. The model parameters are estimated and local phase portraits are constructed for the initial four-dimensional migration-population model. The transition from the vector ordinary differential equation to the corresponding stochastic differential equation is performed. The structure of the stochastic model is described on the basis of applying the method of constructing self-consistent stochastic models. As a tool for the study of population-migration models, a software package is used to solve numerically the differential equations systems using modified Runge–Kutta methods. The software package allows numerical experiments based on the implementation of algorithms for generating trajectories of multidimensional Wiener processes and multipoint distributions and algorithms for solving stochastic differential equations. The comparative analysis of the computer research results obtained for stochastic models is carried out. The properties of migration-population systems in deterministic and stochastic cases are characterized. The comparison of the results obtained for the three-dimensional and four-dimensional cases is carried out. The effects inherent in models with migration flows are revealed.The obtained results can be applied to the problems of modeling and forecasting the behavior of multidimensional systems describing the migration flows. Key words and phrases: computer modeling, nonlinear models of migration flows, stochasticization of one-step processes, symbolic representation algorithms, symbol computing libraries. Copyright © 2019 for the individual papers by the papers’ authors. Copying permitted for private and academic purposes. This volume is published and copyrighted by its editors. In: K. E. Samouylov, L. A. Sevastianov, D. S. Kulyabov (eds.): Selected Papers of the IX Conference “Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems”, Moscow, Russia, 19-Apr-2019, published at http://ceur-ws.org Demidova A. V. et al. 27 1. Introduction The construction and study of the models with migration are devoted to the work of various authors (see, for example, [1–3]). In [2–6] migration flows in deterministic population is considered. The analysis of distributed multidimensional population models taking into account cross-migration is carried out in papers [7–10]. Migration mechanisms can be described by both linear and non-linear functions, and lead to different effects [2, 11]. In [3, 5, 11, 12] the issues of qualitative behavior and sustainability of population-migration models are studied. In [13–15] a qualitative analysis of a three-dimensional non-deterministic model with migration is performed. To study the model, a combination of known methods of synthesis and analysis of models and a method for constructing stochastic self-consistent models are used [16]. In [12] a methodological support is developed for the analysis and synthesis of multidimensional nonlinear dynamic models describing migration flows taking into account the effects of broadband parametric and additive noise. The stability of stationary states is studied and the effects obtained for stochastic models are interpreted. The model examples show a comparison of the migration-population systems properties in deterministic and stochastic cases and the effects due to stochastic broadband perturbations are revealed. When modeling population-migration systems, various software tools are used that present wide possibilities for building computer models and carrying out computational experiments. However, many software products do not contain libraries for numerical and symbolic calculations and do not have sufficient computational complexity. In this regard, in the study of the population-migration systems models, the application of mathematical packages and general-purpose programming languages [17–19] is relevant. In [20] sufficient conditions are proposed for the uniform stability of the equilibrium states of a four-dimensional non-linear model of population dynamics. In [21] the Maxima mathematical package is used to study the stability of this model. One of the instrumental software tools for studying population-migration models is a software package for the numerical solution of differential equations systems using modified Runge – Kutta methods. The specified complex was developed in [22, 23]. The software package allows numerical experiments based on the implementation of algorithms for generating trajectories of multidimensional Wiener processes and multipoint distributions and algorithms for solving stochastic differential equations. The several types of four-dimensional models with migration flows are studied in this paper. A comparative analysis of the computer research results obtained for three- dimensional and four-dimensional stochastic models with migration flows is carried out. A comparison is made of the qualitative properties of four-dimensional models, taking into account changes in migration rates, as well as intraspecific and interspecific interaction coefficients. The properties of the models in the deterministic and stochastic cases are characterized. As a tool for the study of the models, a software package in the Python language using the NumPy and SciPy libraries is used. 2. Deterministic Migration Models A nonlinear model is considered that is described by a system of ordinary differential equations of the form: 𝑥˙1 = 𝑥1 − 𝑥21 − 𝑝13 𝑥1 𝑥3 − 𝑝14 𝑥1 𝑥4 + 𝛽𝑥2 − 𝛾𝑥1 , 𝑥˙2 = 𝑥2 − 𝑥22 + 𝛾𝑥1 − 𝛽𝑥2 , (1) 𝑥˙3 = 𝑥3 − 𝑝31 𝑥1 𝑥3 − 𝑥23 − 𝑝34 𝑥3 𝑥4 , 𝑥˙4 = 𝑥4 − 𝑝41 𝑥1 𝑥4 − 𝑝43 𝑥3 𝑥4 − 𝑥24 . The following notation is used in (1): 𝑥1 , 𝑥3 and 𝑥4 is the population density of competing species in the area of the forma 𝑥1 (area 1), 𝑥2 is the population density 28 ITTMM—2019 in the area of the form 𝑥2 (area 2), 𝑝13 , 𝑝14 , 𝑝31 , 𝑝34 , 𝑝41 , 𝑝43 > 0 are the coefficients of species competition in the area of 1, 𝛽 > 0 and 𝛾 > 0 are the coefficients of the migration of species between two areas, while area 2 is a refuge and 𝛽 ̸= 𝛾. In the model (1) the first two equations describe the dynamics of one type taking into account migration processes. The first equation sets the dynamics in the first area, the second one sets the dynamics in the second area. The third and fourth equations describe the dynamics of the second and third types, interacting as competitors with the first type in the first area. Previously, we studied the three-dimensional model (two species, the first species migrated to another area, and in the first area it competed with the second type). In the case when 𝑝13 = 𝑝31 , 𝑝14 = 𝑝41 , 𝑝34 = 𝑝43 , the model (1) takes view 𝑥˙1 = 𝑥1 − 𝑥21 − 𝑝13 𝑥1 𝑥3 − 𝑝14 𝑥1 𝑥4 + 𝛽𝑥2 − 𝛾𝑥1 , 𝑥˙2 = 𝑥2 − 𝑥22 + 𝛾𝑥1 − 𝛽𝑥2 , (2) 𝑥˙3 = 𝑥3 − 𝑥23 − 𝑝13 𝑥1 𝑥3 − 𝑝34 𝑥3 𝑥4 , 𝑥˙4 = 𝑥4 − 𝑥24 − 𝑝14 𝑥1 𝑥4 − 𝑝34 𝑥3 𝑥4 . At the same migration rates 𝛽 = 𝛾 = 𝜀 from (2) we get a model of the form 𝑥˙1 = 𝑥1 − 𝑥21 − 𝑝13 𝑥1 𝑥3 − 𝑝14 𝑥1 𝑥4 + 𝜀(𝑥2 − 𝑥1 ), 𝑥˙2 = 𝑥2 − 𝑥22 + 𝜀(𝑥2 − 𝑥1 ), (3) 𝑥˙3 = 𝑥3 − 𝑥23 − 𝑝13 𝑥1 𝑥3 − 𝑝34 𝑥3 𝑥4 , 𝑥˙4 = 𝑥4 − 𝑥24 − 𝑝14 𝑥1 𝑥4 − 𝑝34 𝑥3 𝑥4 . In this paper, deterministic models (1)–(3) and their stochastic generalizations are studied. A three-dimensional model that describes two types of dynamics, the first species migrating to another area, and in the first area competing with the second type, has the form: 𝑦˙1 = 𝑦1 − 𝑦12 − 𝑝13 𝑦1 𝑦3 + 𝛽𝑦2 − 𝛾𝑦1 , 𝑦˙2 = 𝑦2 − 𝑦22 + 𝛾𝑦1 − 𝛽𝑦2 , (4) 𝑦˙3 = 𝑦3 − 𝑦32 − 𝑝13 𝑦1 𝑦3 , where 𝑦1 and 𝑦3 are the density of populations of competing species in the area of the form 𝑦1 (area 1), 𝑦2 is the population density in the area type 𝑦2 (area 2), 𝑝13 > 0 is the coefficient of competition of species in the area 1. The results of the analysis of the (4) model and the generalized stochastic model are presented in [13–15]. The model (4) is a generalization of the model considered in [3], in the case when the migration rates are different. The numerical experiment for the model (4) is carried out with the help of the developed software for the study and numerical solution of systems of ordinary and stochastic differential equations using the Runge–Kutta method [22, 23]. The library is written in Python using the NumPy and SciPy modules. The models (1)–(3) and their stochastic generalizations are studied using this software package. For the numerical experiment of the deterministic model, the following sets of parameters are chosen: the initial values 𝑦1 (0) = 0.5, 𝑦2 (0) = 1.0, 𝑦3 (0) = 0.5 for the three-dimensional model and 𝑥1 (0) = 0.5, 𝑥2 (0) = 1.0, 𝑥3 (0) = 0.5,𝑥4 = 0.7 for a Demidova A. V. et al. 29 four-dimensional model, the values of the parameters 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4, 𝛽 = 1.5, 𝛾 = 0.2. Fig. 1 presents solution trajectories for (2) and (4) models with the specified initial values and parameter values in the time interval [0,40]. The effect of an additional competitor on the migration flow is shown. The introduction of an additional competitor corresponding to the phase variable 𝑥4 , leads to a decrease in both the density of the migrating population in both areas and the density of the competitor corresponding to the variable 𝑥3 . Figure 1. Model solution trajectories (2) and (4) with (𝑦1 (0), 𝑦2 (0), 𝑦3 (0)) = (0.8, 0.1, 1.0), (𝑥1 (0), 𝑥2 (0), 𝑥3 (0), 𝑥4 (0)) = (0.8, 0.1, 1.0, 1.0), 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4, 𝛽 = 1.5, 𝛾 = 0.2. With the values of the parameters 𝑝13 = 0.9, 𝑝14 = 0.4, 𝑝34 = 0.8, 𝑝34 = 0.6, 𝑝41 = 0.5, 𝑝43 = 0.7, 𝛽 = 0.25, 𝛾 = 0.3 in the model (1) there are the following equilibria: 𝑃1 (0, 0, 0, 0), 𝑃2 (0, 0, 0, 1), 𝑃3 (0, 0, 1, 0), 𝑃4 (0, 0, 0.69, 0.517), 𝑃5 (0.63, 0.949, 0.4957, 0), 𝑃6 (0.75, 0.98, 0.043, 0.595). Here 𝑃1 , 𝑃 5 are unstable nodes, 𝑃2 − 𝑃4 are saddles, 𝑃6 is a stable node. Consider the models (1) and (2) with the same set of migration velocity values for these models, but with different sets of intraspecific and interspecific interaction coefficients. Solution trajectories for the model (1) with the values of the parameters 𝑝13 = 0.9, 𝑝14 = 0.4, 𝑝31 = 0.8, 𝑝34 = 0.6, 𝑝41 = 0.5, 𝑝43 = 0.7 and the trajectories of solutions for the model (2) with 𝑝13 = 1, 𝑝14 = 0.6, 𝑝34 = 0.7 in the time interval [0, 25] is presented in fig. 2. For models (1) and (2) the initial values are (𝑥1 (0), 𝑥2 (0), 𝑥3 (0), 𝑥4 (0)) = (0.3, 0.7, 1.1, 0.8), values of migration speeds 𝛽 = 0.25, 𝛾 = 0.3. In fig. 2 a record of the form 𝑥𝑖 (1), 𝑖 = 1, 2, 3, 4, indicates the trajectory corresponding to the phase variable 𝑥𝑖 for the model (1), a record of the form 𝑥𝑖 (2) denotes the trajectory corresponding to the phase variable 𝑥𝑖 for the model (2). As a result of a comparative analysis of the solution trajectories of (1) and (2) presented in fig. 2: 1) in the model (2) the population of 𝑥4 is dying out, the population density of 𝑥1 is increasing, and the population density of 𝑥3 is decreasing; 2) for the models (1) and (2) favorable conditions are created in the second area, and the population density of 𝑥2 increases; 3) in the model (2) competition intensified, which led to the extinction of the 𝑥4 population and the improvement of living conditions for the 𝑥1 and 𝑥3 populations. For the values of the parameters 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4, 𝜀 = 0.3 in model (3) there are the following equilibrium positions: 𝑃1 (0, 0, 0, 0), 𝑃2 (0, 0, 0, 1), 𝑃3 (0, 0, 1, 0), 𝑃4 (0, 0, 0.42, 0.42), 𝑃5 (1, 1, 0, 0), 𝑃6 (0.76, 0.94, 0, 0.62). Here 𝑃2 − 𝑃5 are saddles, 𝑃1 is an unstable node, 𝑃6 is a stable node. 30 ITTMM—2019 Figure 2. Model solution trajectories (1) and (2) with (𝑥1 (0), 𝑥2 (0), 𝑥3 (0), 𝑥4 (0)) = (0.3, 0.7, 1.1, 0.8), 𝛽 = 0.25, 𝛾 = 0.3. Consider the models (2) and (3) taking into account the coincidence of the intraspecific and interspecific interaction for these models. Note that in the model (2) in the general case 𝛽 ̸= 𝛾, and in the model (3) 𝛽 = 𝛾. Solution trajectories for models (2) and (3) with the values of the parameters 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4 and the initial values (𝑥1 (0), 𝑥2 (0), 𝑥3 (0), 𝑥4 (0)) = (0.5, 1, 0.8, 1) in the time interval [0, 25] are presented in fig. 3. For the model (2) the values of migration speeds 𝛽 = 1.5, 𝛾 = 0.2, for the model (3) the value of migration speed 𝜀 = 0.3. In fig. 3 the record of the form 𝑥𝑖 (3) denotes the trajectory corresponding to the phase variable 𝑥𝑖 for the model (3). Figure 3. Model solution trajectories (2) and (3) (𝑥1 (0), 𝑥2 (0), 𝑥3 (0), 𝑥4 (0)) = (0.5, 1, 0.8, 1), 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4, 𝛽 = 1.5, 𝛾 = 0.2, 𝜀 = 0.3. As a result of a comparative analysis of the solution trajectories of the (2) and (3) models shown in fig. 3, the following conclusions are obtained: 1) in the model (3)the population density 𝑥4 increased insignificantly, and the pop- ulation density 𝑥1 decreased, while the population 𝑥3 in both models is rapidly dying out; 2) as the migration rate increases, the model (3) in the second area creates more favorable conditions than in the model (2), so the population density 𝑥2 increases significantly. For a series of parameter sets, stationary states of the model (3) are found. It is established that for the values of the parameters 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4, 𝛽 = 1.5, 𝛾 = 0.2 in the model (2) the following stationary states exist: 𝑃1 (0, 0, 0, 0), 𝑃2 (0, 0, 0, 1), Demidova A. V. et al. 31 𝑃3 (0, 0, 1, 0), 𝑃4 (1.18, 0.29, 0, 0), 𝑃5 (0, 0, 0.42, 0.42), 𝑃6 (0.93, 0.25, 0, 0.53). In addition, the system is linearized in the vicinity of stationary states and the stability on the first approach is investigated. For the model (2) the projections of the phase portraits (𝑥1 , 𝑥2 ) are shown in fig. 4. The projections of the phase portraits (𝑥1 , 𝑥2 ) in the vicinity of equilibrium positions are shown in fig. 5. The projections of the phase portraits (𝑥1 , 𝑥3 ) are shown in fig. 6. a) 𝑃1 b) 𝑃2 c) 𝑃3 d) 𝑃4 d) 𝑃5 d) 𝑃6 Figure 4. Phase portraits for the model (2). 32 ITTMM—2019 a) 𝑃1 b) 𝑃2 c) 𝑃3 d) 𝑃4 d) 𝑃5 d) 𝑃6 Figure 5. Phase portraits of the projections (𝑥1 , 𝑥2 ) for the model (1) in the vicinity of the equilibrium positions 𝑃1 –𝑃6 . 3. Stochastic Migration Models Questions of constructing and studying a stochastic model for a three-dimensional system with regard to migration are considered in the works [13–15]. In this paper, stochastization of (1)–(3) models is carried out using the method of constructing self- consistent stochastic models [16]. The basis of this method is combinatorial methodology. This method involves recording the system under study as an interaction scheme, i.e. a symbolic record of all possible interactions between elements of the system. For this, the Demidova A. V. et al. 33 a) 𝑃1 b) 𝑃2 c) 𝑃3 d) 𝑃4 d) 𝑃5 d) 𝑃6 Figure 6. Phase portraits for the model (3). 34 ITTMM—2019 system state operators and the system state change operator are used. Then we can get the drift and diffusion coefficients for the Fokker–Planck equation, which allows to write the equation itself and its equivalent stochastic differential equation in the Langevin form. However, an increase in the dimension and complexity of the system under study leads to the complexity of the analytical conclusion of the necessary coefficients of the Fokker– Planck equation. As a solution to this problem, a software implementation is developed for obtaining the coefficients of the Fokker–Planck equation from interaction schemes using a symbolic computation system. This software implementation is a modification of the stochastization method for one-step processes in the computer algebra system described in [24]. This implementation is introduced as a module into a software package developed earlier for the numerical study of deterministic and stochastic models [22]. The following is the algorithm for obtaining the symbolic notation of a stochastic differential equation (Algorithm 1). Algorithm 1. Getting the system of differential equations from the interaction scheme Initial parameters: interaction scheme. Result: system of differential equations in the form of Langevin. Start: 1. Getting the system state operators from the interaction scheme. 2. Getting the change of the system state. 3. Getting the transition intensities. 4. Record the coefficients of the Fokker–Planck equation. 5. Record the system of differential equations. end To implement the described algorithm, SymPy [17], computer computing system is used, which is a powerful symbolic computation library for the Python language. In addition, the output data obtained using the SymPy library can be transferred for numerical calculations using the NumPy [19] library and SciPy [18]. To obtain a stochastic model corresponding to the models (1)–(3), it is necessary to write the interaction scheme, which is as follows: 𝑎 𝑖 𝑋𝑖 −→ 2𝑋𝑖 , 𝑖 = 1, 4; 𝑏𝑖 𝑋𝑖 + 𝑋𝑖 −→ 𝑋𝑖 , 𝑖 = 1, 4; 𝑐𝑖𝑗 (5) 𝑋𝑖 + 𝑋𝑗 −−→ 𝑋𝑗 , 𝑖, 𝑗 = 1, 3, 4, 𝑖 ̸= 𝑗; 𝑑 𝑖 𝑋𝑖 −→ 𝑋𝑗 , 𝑖, 𝑗 = 1, 2, 𝑖 ̸= 𝑗. The values of the coefficients for each of the models are given in the table 1. The first line corresponds to the natural reproduction of the species in the absence of other factors, the second line symbolizes intraspecific competition, and the third one is interspecific competition. The last line corresponds to the description of the species migration process from one area to another. Then, expressions for the coefficients of the Fokker–Planck equation are obtained for this interaction scheme using the developed software (fig. 7). Further, the coefficients are transferred to the appropriate module of the software complex for the numerical solution of the constructed stochastic differential equation. For the numerical experiment of the obtained stochastic model, the parameters are chosen, the same as for the numerical analysis of the deterministic model (3). Fig. 8 presents the trajectories of average values for 100 realizations for the specified initial values and parameter values in the time interval [0, 40]. A numerical experiment showed the closeness of trajectories for stochastic and deterministic cases. Similarly to the deterministic case, the mean values trajectories of different implementations of the stochastic differential equation solutions go to the Demidova A. V. et al. 35 Table 1 Coefficients for the interaction scheme (5) Model (1) Model (2) Model (3) 𝑎𝑖 𝑎𝑖 = 1, 𝑖 = 1, 4 𝑎𝑖 = 1, 𝑖 = 1, 4 𝑎𝑖 = 1, 𝑖 = 1, 4 𝑏𝑖 𝑏𝑖 = 1, 𝑖 = 1, 4 𝑏𝑖 = 1, 𝑖 = 1, 4 𝑏𝑖 = 1, 𝑖 = 1, 4 𝑐𝑖𝑗 𝑐13 = 𝑝13 , 𝑐14 = 𝑐13 = 𝑐31 = 𝑝13 , 𝑐13 = 𝑐31 = 𝑝13 , 𝑝14 , 𝑐31 = 𝑝31 , 𝑐14 = 𝑐41 = 𝑝14 , 𝑐14 = 𝑐41 = 𝑝14 , 𝑐34 = 𝑝34 , 𝑐41 = 𝑐34 = 𝑐43 = 𝑝34 𝑐34 = 𝑐43 = 𝑝34 𝑝41 , 𝑐43 = 𝑝43 , 𝑑𝑖 𝑑1 = 𝛾, 𝑑2 = 𝛽 𝑑1 = 𝛾, 𝑑2 = 𝛽 𝑑1 = 𝑑2 = 𝜀 Figure 7. The coefficients of the equation of Fokker–Planck Figure 8. Model Solution Trajectories (5) with (𝑥1 (0), 𝑥2 (0), 𝑥3 (0), 𝑥4 (0)) = (0.8, 0.1, 1.0, 1.0), 𝑝13 = 1.2, 𝑝14 = 0.5, 𝑝34 = 1.4, 𝛽 = 1.5, 𝛾 = 0.2. stationary mode, retain the character and values of the stationary state are weakly distinguished. A series of numerical experiments is carried out using different sets of model parame- ters. For each set, we analyzed the stationary states and constructed the trajectories of solutions. Taking into account the characteristic cases, the peculiarities of the influence 36 ITTMM—2019 of an additional competitor on the dynamics of the migrating population are revealed and the effects of stochastization of one-step processes are described. 4. Conclusions A computer study of nonlinear models with migration flows made it possible obtain the results of numerical experiments on trajectory search and parameter estimation in the case of high dimensionality of models, as well as to reveal effects due to stochasticization. The estimation of the influence of an additional competitor type on the dynamics, stability and level of migration in the multidimensional model is obtained. A comparative analysis of computer simulation results showed that with an increase in the rate of migration in the second area, more favorable conditions are created, which leads to a significant increase in the number of the second population. At the same time, competition in the first area is increasing, which leads to the extinction of one of the population. The comparison of the results obtained for the three-dimensional and four-dimensional cases. The software package developed in the Python language using the NumPy and SciPy libraries is demonstrated sufficient efficiency for computer studies of multidimensional nonlinear migration models. Numerical experiments using problem-oriented software show the closeness of the trajectories types for stochastic and deterministic cases. The results obtained can be applied to the problems of computer modeling of ecological, demographic and socio-economic systems. Acknowledgments The publication has been prepared with the support of the “RUDN University Program 5-100” (Demidova A.V., numerical analysis). References 1. N. Tkachenko, J. D. Weissmann, W. P. Petersen, G. Lake, C. P. E. Zollikofer, S. Callegari, Individual-based modelling of population growth and diffusion in discrete time 2 (4) (2017) e0176101. doi:10.1371/journal.pone.0176101. 2. H. C. Tuckwell, A study of some diffusion models of population growth, Theoretical Population Biology 5 (3) (1974) 345–357. doi:10.1016/0040-5809(74)90057-4. 3. X.-a. Zhang, L. Chen, The linear and nonlinear diffusion of the competitive Lotka–Volterra model, Nonlinear Analysis 66 (2007) 2767–2776. doi:10.1016/j.na. 2006.04.006. 4. Y. Svirezhev, Nonlinear waves, dissipative structures and disasters in ecology, Moscow: Science, 1987. 5. Z. Lu, Y. Takeuchi, Global asymptotic behavior in single-species discrete diffu- sion systems, Journal of Mathematical Biology 32 (1) (1993) 67–77. doi:10.1007/ BF00160375. 6. J. Cui, L. Chen, The effect of diffusion on the time varying logistic population growth, Computers & Mathematics with Applications 36 (1998) 1–9. doi:10.1016/ S0898-1221(98)00124-2. 7. L. Chen, A. Jüngel, Analysis of a multi-dimensional parabolic population model with strong cross-diffusion, SIAM J. Math. Anal. 36 (2004) 301–322. doi:10.1137/ S0036141003427798. 8. Z. Wen, S. Fu, Global solutions to a class of multi-species reaction-diffusion systems with cross-diffusions arising in population dynamics, J. Comput. Appl. Math. 230 (2) (2009) 34–43. doi:10.1016/j.cam.2008.10.064. 9. N. Zamponi, A. Jüngel, Analysis of degenerate cross-diffusion population models with volume filling, Annales de l’Institut Henri Poincaré C, Analyse non linéaire 34 (1) (2017) 1–29. doi:10.1016/j.anihpc.2015.08.003. Demidova A. V. et al. 37 10. X. Chen, E. S. Daus, A. Jüngel, Global existence analysis of cross-diffusion popu- lation systems for multiple species, Archive for Rational Mechanics and Analysis 227 (2) (2018) 715–747. doi:10.1007/s00205-017-1172-6. 11. L. J. S. Allen, Persistence and extinction in single-species reaction-diffusion models, Bulletin of Mathematical Biology 45 (2) (1983) 209–227. doi:10.1007/BF02462357. 12. I. N. Sinitsyn, O. V. Druzhinina, O. N. Masina, Analytical modeling and stability analysis of nonlinear broadband migration flow, Nonlinear world 16 (3) (2018) 3–16. 13. A. Demidova, O. Druzhinina, M. Jacimovic, O. Masina, Construction and analysis of nondeterministic models of population dynamics, Vol. 678, 2016, pp. 498–510. doi:10.1007/978-3-319-51917-3_43. 14. A. V. Demidova, O. V. Druzhinina, O. N. Masina, Design and stability analysis of nondeterministic multidimensional populations dynamics models, Vol. 1995, 2017, pp. 14–21. 15. A. V. Demidova, O. V. Druzhinina, M. Jacimovic, O. N. Masina, N. Mijajlovic, Synthesis and analysis of multidimensional mathematical models of population dynamics, in: 2018 10th International Congress on Ultra Modern Telecommu- nications and Control Systems and Workshops (ICUMT), 2018, pp. 361–366. doi:10.1109/ICUMT.2018.8631252. 16. A. V. Demidova, M. N. Gevorkyan, A. D. Egorov, D. S. Kulyabov, A. V. Ko- rolkova, S. L.A., Influence of stochastization on one–step models, RUDN Journal of Mathematics, Information Sciences and Physics (1) (2014) 71–85. 17. R. Lamy, Instant SymPy Starter, Packt Publishing, 2013. 18. T. E. Oliphant, Python for scientific computing, Computing in Science and Engg. 9 (3) (2007) 10–20. doi:10.1109/MCSE.2007.58. 19. T. E. Oliphant, Guide to NumPy, 2nd Edition, CreateSpace Independent Publishing Platform, USA, 2015. 20. E. D. Tarova, Application the divergent method to the stability analysis of the dynamic population nonlinear model, Nonlinear world 15 (3) (2017) 59–64. 21. E. D. Tarova, Qualitative research of the population dynamics four-dimensional model with the help of lyapunov forst method, Nonlinear world 16 (3) (2018) 17–24. 22. E. G. Eferina, A. V. Korolkova, M. N. Gevorkyan, D. S. Kulyabov, L. A. Sevastyanov, One–step stochastic processes simulation software package, RUDN Journal of Math- ematics, Information Sciences and Physics (3) (2014) 46–59. 23. M. N. Gevorkyan, T. R. Velieva, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastyanov, Stochastic runge–kutta software package for stochastic differential equations, in: W. Zamojski, J. Mazurkiewicz, J. Sugier, T. Walkowiak, J. Kacprzyk (Eds.), De- pendability Engineering and Complex Systems, Springer International Publishing, Cham, 2016, pp. 169–179. doi:10.1007/978-3-319-39639-2_15. 24. M. N. Gevorkyan, A. V. Demidova, T. R. Velieva, A. V. Korol’kova, D. S. Kulyabov, L. A. Sevast’yanov, Implementing a method for stochastization of one-step processes in a computer algebra system, Programming and Computer Software 44 (2) (2018) 86–93. doi:10.1134/S0361768818020044.