Numerical Investigation and Factor Analysis of the Spatial-Temporal
Multi-Species Competition Problem
MARIA VASILYEVA1, YOUWEN WANG1, SERGEI STEPANOV2, ALEXEY SADOVSKI1
1Department of Mathematics and Statistics,
Texas A&M University - Corpus Christi,
Corpus Christi, Texas,
UNITED STATES OF AMERICA
2Institute of Mathematics and Information Science,
North-Eastern Federal University,
Yakutsk, Republic of Sakha (Yakutia),
RUSSIA
Abstract: - This work considers the spatial-temporal multi-species competition model. A mathematical model is
described by a coupled system of nonlinear diffusion reaction equations. We use a finite volume approximation
with semi-implicit time approximation for the numerical solution of the model with corresponding boundary
and initial conditions. To understand the effect of the diffusion to solution in one and two-dimensional
formulations, we present numerical results for several cases of the parameters related to the survival scenarios.
We control all non-diffusion parameters, including reproductive growth rate, competition rate, and initial
condition of population density of competing species, and compare the dynamic and equilibrium under regular
diffusion rate and small diffusion rate; we found that competing species with small diffusion rate can reach a
higher equilibrium over the whole geographic domain, but requires more time steps. The random initial
conditions' effect on the time to reach equilibrium is investigated. We control other parameters and examine the
impact of the initial condition of the species population; we found that regardless of the values of initial
conditions in the system, competing species populations will arrive at an equilibrium point. The influence of
diffusion on the survival scenarios is presented. We control other parameters and examine the effect of
diffusion of species; we found that when the ratio of diffusion rates passes some thresholds, the survival status
will change. In real-world problems, values of the parameters are usually unknown yet vary in some range. To
evaluate the impact of parameters on the system stability, we simulate a spatial-temporal model with random
parameters and perform factor analysis for two and three-species competition models. From the perspective of
the numerical experiment, we release control for all parameters and perform factor analysis on simulation
results. We found that the initial population condition has a minimum effect on the final population, which
aligns with the outcome of our controlled numerical experiment on the initial condition. Diffusion is the
dominant factor when diffusion rates are on the same scale as other parameters. This dominant factor aligns
with our controlled numerical experiment on diffusion rate, where the change in diffusion rate leads to different
survival statuses of species. However, when diffusion rates are 1/10 on the scale of other parameters,
reproductive growth rates and competition rates become the dominant factors.
Key-Words: - Spatial-temporal model, Multi-species competition model, System of diffusion-reaction
equations, Numerical investigation, Impact of parameters, Factor analysis
Received: September 19, 2021. Revised: September 6, 2022. Accepted: October 9, 2022. Published: November 4, 2022.
1 Introduction
Understanding the stability of ecosystems is of
fundamental importance to ecology, [1],[2],[3].
Fundamental mathematical models of such systems
are described by a coupled system of ordinary
differential equations (ODEs). The Lotka–Volterra
Competition model (LVC) is a basic model that
describes the dynamics of the species population
competing for some shared resource. The LVC
model has been applied in many areas, including
biological systems, industry, and economics,
[4],[5],[6],[7],[8]. For example, the model can
simulate the marsh ecosystems for the wetlands at
the Nueces River mouth, [9].
The LVC model is based on the logistic
population model
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
731
Volume 21, 2022

󰇛󰇜
where is the size of the population at a given time
and  is the per-capita growth rate.
For the general case of multi-species
competition, we have
󰇛󰇜
 󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜


where󰇛󰇜 is the population of the -th species,
is the -th population growth rate,  is the
interaction coefficient due to competition
(󰇛󰇜compete with 󰇛󰇜), and is the number of
species (equations).
In such systems, a temporal dynamic model can
represent and describe the behaviors of the entire
system. However, real ecosystems interact in
different locations, and spatial structure is essential
and has an enormous impact on the final equilibrium
state, [10],[11]. Systems of the partial differential
equations (PDEs) are used to describe
spatial-temporal systems. A mathematical model is
described by a coupled system of unsteady
nonlinear reaction-diffusion equations. For the
multi-species interaction, we have
󰇛󰇜
 󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜

where the diffusion coefficient. This model
incorporates spatial structure by adding diffusion
terms in equations and considering the system of
equations in the domain , where is the
spatial dimension.
Previous research has put much effort into
modifying the Lotka–Volterra competition model
(LVC). The parameters , , and can be
modeled as functions instead of constants. Diffusion
, more specifically, traveling waves, and their
minimal-speed selection mechanisms in the LV
model can be studied by applying the upper-lower
solution technique on the cooperative system, [12].
The crowding effect of diffusion can also be
modeled within the LVC model, [13]. Advection
rate can be introduced to the LVC model as a
supplement of diffusion to bring much richer
phenomena, [14],[15]. Carrying capacity can be
replaced by function instead of constant and can
vary between species; connection topology can be
modeled as a 'loop,' 'star,' 'chain,' or 'full' connection
when there are more than two species in the system,
[16]. In addition to making parameters endogenous,
some other efforts have also been made to modify
the LVC model to simulate real-world problems.
For example, simulation adding small immigration
into the prey or predator population can stabilize the
LV system, [17]. The random fluctuating
environment can also be modeled within the LV
model to show how switching between
environments can make survival harder, [18].
Stochastic noises can be introduced to the LV
population model and are presented to play an
essential role in the permanence and
characterization of the system, [19].
Multiple methods are used for solving the LVC
model, such as Haar wavelet (HW), Adams-
Bashforth-Moulton (ABM) methods, [20], and finite
element method. After simulation, multiple methods
have been used to analyze the simulated species
interaction data. Real-world data can be used to
evaluate the multi-species interaction model, [21].
For example, the population dynamics model of
individual reefs can be compared with data on coral
reefs in Pilbara, [22]. Numerical models and
machine learning can be combined to identify the
factors influencing Alexandrium catenella blooms,
[23].
In this work, we consider a spatial-temporal
multi- species competition model in one- and two-
dimensional formulations. For the numerical
solution of the model with corresponding boundary
and initial conditions, we construct a finite volume
approximation with a semi-implicit time
approximation. When two or more species compete
for the same limited food source or in some way
influence each other's growth, one or several of the
species usually becomes extinct. To understand the
diffusion effect of the solution, we perform a
numerical investigation for several cases of the
parameters related to the survival scenarios. The
influence of the diffusion and initial conditions on
the survival scenarios is presented for two and three-
species competition models. It turns out that
equilibrium depends only on system parameters
(birth rate, competition, predation, etc.) and does not
depend on the initial conditions. The time it takes
for the system to be in the steady state depends on
how far initial populations of species are located
from the equilibrium point in the phase space of the
model. After understanding the effect of parameters
and initial conditions on the test problems, we
simulate the spatial-temporal model with random
input parameters, i.e., releasing the hold for all
parameters. Simulations with different combinations
of parameters and initial conditions support the
hypothesis that such systems reach equilibrium
sooner or later. At last, we apply factor analysis to
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
732
Volume 21, 2022
evaluate the impact of parameters on the system
stability.
The main novelty of this work is in performing
factor analysis methods for PDE-based data sets.
For the PDE-based mathematical models, the
finding of an analytical solution is limited by special
cases and cannot be performed for the general
spatial-temporal model and consider complex
diffusion and domain structure and effect of the
various boundary conditions, [24],[25]. In this
paper, we construct a discrete system using the
Finite Volume Method with semi-implicit time
approximation. Simulation of the ecological system
with a PDE-based discrete system allows us to (1)
take into the spatial distribution of the species
concerning boundary conditions and (2) use the
limited information about ecosystem parameters by
a given estimated range instead of acquiring specific
value which is usually unknown for the real
systems. The presented approach allows us to
identify which ecosystem parameters are dominant,
only given the range of each parameter. For this
purpose, we use a factor analysis for (a)
understanding the system and (b) predicting the
system. From the goal of understanding the system,
the factor analysis approach allows us to understand
the relative importance of parameters when
information about specific parameter values is
limited. For example, given ranges of parameters in
a food web containing multiple species, we can
identify what parameters (e.g., birth rate,
competition rate, or diffusion rate) or what species
is the" critical factor" that influences the final
solution of the population dynamics. It can be
applied not only to two and three-species competing
systems but also to systems corresponding to more
giant food web (state-of-the-art simulations in terms
of the number of species, [26],[27],[28] simulated
system of four species). Currently, in the ecology
and environmental field, factor analysis is applied
more to real-world data, [29],[30]. Still, we have
shown that it can also be used for PDE-based
simulation and can be applied to more complex
ecosystem simulation when computational power
increase. We can also apply this simulation
approach with the factor analysis method to study
catastrophic events in future works. The extinction
property has been studied analytically, [31],[32] and
numerically, [33],[34], and the presented method
can bring about a novel perspective to study the
topic numerically. From the goal of predicting the
system, the factor analysis approach can be used for
dimension reduction of input data of the predictive
model. We can remove less dominant factors to
speed up the prediction process or replace the input
of the predictive model with the main components.
Some other feature selection methods can also serve
the purpose of dimension reduction [35]. Multiple
predictive models can be applied, afterward, such as
beta regression, [36], explainable prediction, [37],
and neural networks, [38]. In future works, we will
concentrate on the following directions: (1) a real-
world range of parameters for a given application;
(2) apply the result of factor analysis to build
machine learning predictive models, and predict
final population density, as well as time to reach
equilibrium; and (3) use the model to predict the
equilibrium population density and simulate future
catastrophic event (given the assumption that the
system is at equilibrium then catastrophe strokes).
The paper is organized as follows. Section 2 de-
scribes the mathematical model with fine-scale
approximation using the finite volume method and a
semi-implicit scheme for time approximation.
Section 3 presents numerical results for two and
three-species competition models in one and two-
dimensional formulations. In Section 4, we simulate
a spatial-temporal model with random parameters
and perform factor analysis to evaluate the impact of
parameters on the system stability. The paper ends
with a conclusion.
2 Mathematical Models with
Approximation by Space and Time
A mathematical model is the simplified
representation of the complex real-world objects or
systems used to understand the system's complex
interactions to predict possible outcomes of different
changes. In mathematical biology and ecology, the
fundamental mathematical model is the
Lotka-Volterra model, which describes the temporal
dynamic of the species population in various
ecosystem models. However, real ecosystems
interact in different locations, and spatial structure is
essential and dramatically impacts the population.
In this work, we consider the spatial-temporal
multi-species competition problem in one and two-
dimensional domains. Let  be the
computational domain, where for  we have a
one- dimensional case, and for  we obtain
two- dimensional problem. The population of
species in Ω is described by a coupled system of
nonlinear diffusion - reaction equation.
󰇛󰇜
 󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜

WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
733
Volume 21, 2022
 ( 1 )
where where is the number of
species (equations). Here 󰇛󰇜 is the population of
the -th species, is the -th population
reproductive growth rate,  is the interaction
coefficient due to competition (󰇛󰇜 compete with
󰇛󰇜) and is the diffusion coefficient, [39],[40].
The system of equations is considered with
initial conditions
󰇛󰇜
󰇛󰇜   ( 2 )
and boundary conditions
󰇛󰇜   ( 3 )
In this work, we consider the following special
cases:
Two-species competition
󰇛󰇜
 󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜
󰇛󰇜
 󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜
Three-species competition
󰇛󰇜
 󰇛󰇜󰇛󰇜󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜
󰇛󰇜
 󰇛󰇜󰇛󰇜󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜
󰇛󰇜
 󰇛󰇜󰇛󰇜󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜󰇛󰇜
Analytical solution of such spatial-temporal
models is possible only for some simplified cases. A
numerical simulation of the model required the
construction of an accurate approximation for
spatial variation (diffusive operators) and time
approximation. Three common space approximation
techniques exist: The Finite Difference Method, the
Finite Element Method, and the Finite Volume
Method. The main advantage of the Finite Element
Method is the accurate representation of the solution
in the complex computational domains by
constructing unstructured grids. For simplified
geometries, the Finite Volume Method is the
common choice for space approximation because it
provides a conservative approximation.
Let 󰇟󰇠 () be the computational
domain, where we have an interval for the one-
dimensional case and a square domain for the two-
dimensional problem. For the numerical solution of
the model (1) with boundary conditions (2) and
initial conditions (3), we construct a computational
mesh and use a Finite Volume Method for
approximation.
Let be the structured grid for domain

where be the square cell and is the number
of cells, [41]. Here 󰇟󰇛󰇜󰇠 for one-
dimensional case, where . For the two-
dimensional case, we have 󰇟󰇛󰇜󰇠
󰇟󰇛󰇜󰇠 with  and ,
where  is the global cell indexing and
is the number of nodes in x and y
direction. Here we have for the one-
dimensional case and for the two-
dimensional case.
Fig. 1: Illustration of the computational grid for
one and two dimensional cases. Ki is the cell
and Eij is the interface between to cells K
i
and
K
j
To write an approximation by space using the
Finite Volume Method, we integrate equation (1)
over the cell volume and obtain the following semi-
discrete form:
󰇛󰇜


󰇛󰇜
󰇛󰇜󰇛󰇜󰇛󰇜
 󰇛󰇜
( 4 )
In the Finite Volume Method}, we set
󰇛󰇜
󰇛󰇜
where is cell volume and 󰇛󰇜 is the cell
average value of the function 󰇛󰇜 on cell .
For the diffusion operator approximation, we use
a classic two-point flux approximation
󰇛󰇜 󰇛󰇜

󰇡󰇛󰇜󰇛󰇜󰇢
with 
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
734
Volume 21, 2022
where  is the distance between cell center points
and ,  is the length of the interface
between two cells and . Note that, for a
structured uniform grid, we have , 
for the one-dimensional case and  for the
two-dimensional problem.
Therefore, we have
󰇛󰇜
 󰇡󰇛󰇜󰇛󰇜󰇢
󰇛󰇜󰇡󰇛󰇜󰇢󰇛󰇜
 󰇛󰇜
( 5 )
For the time derivative, we use a backward Euler
approximation 󰇛󰇜
 󰇛󰇜󰇛󰇜
where is the time step size and 󰇛󰇜 is the solution
from the previous time layer.
To remove large time step restrictions regarding
the diffusion operator, we approximate the diffusion
part using the solution from the current time level.
However, we use an explicit approximation for the
nonlinear reaction part to linearize the problem.
Finally, we obtain the following discrete problem
using a semi-implicit time approximation scheme,
[42],[43],[44].
󰇛󰇜󰇛󰇜
󰇡󰇛󰇜󰇛󰇜󰇢
󰇛󰇜󰇡󰇛󰇜󰇢󰇛󰇜
 󰇛󰇜
( 6 )
Note that we obtain an uncoupled system of
equations and can solve the equation for each
component separately. Algorithm:
Set initial conditions for each species :
󰇛󰇜
󰇛󰇜  
where is the number of species and is the
number of grid cells.
Solve the linear system of equations on each
time layer till all species solutions converge
to the final equilibrium state
󰇛󰇜
where 󰇛󰇜󰇡󰇛󰇜
󰇛󰇜󰇢 is the vector of
solution with size and
󰇛󰇜󰇡󰇛󰇜󰇛󰇜󰇢
󰇛󰇜
󰇛󰇜󰇡󰇛󰇜󰇢
󰇛󰇜
 󰇛󰇜
for .
Here we have the right-hand side vector that
depends on the solution from the previous time layer.
The matrix is the positive definite and symmetric.
Furthermore, we have a tridiagonal matrix for the
one-dimensional case and a five-diagonal matrix for
the two-dimension. Implementation is performed
using Python programming language using a
standard solver from the NumPy package, [45],[46].
In most results, we perform simulations till both
populations reach equilibrium, 󰇛󰇜󰇛󰇜
with for each .
3 Numerical Results
Before we make factor analysis for randomly
generated parameters, we consider some special
cases with parameters (growth rate "", competition
efficiency " ", and diffusion rate " ") fixed as
constant and examine results such that either one
species survives, two species survive, or three
species survive.
Let 󰇛󰇜
󰇛󰇜
  
  
  
We present numerical results for multi-species
competition in domain . We consider the following
domain and boundary conditions configurations:
: 󰇟󰇠 with zero (fixed) boundary
conditions on .
󰇛󰇜: 󰇟󰇠 with zero (fixed)
boundary conditions on .
󰇛󰇜: 󰇟󰇠 with zero (fixed)
boundary conditions on the left and right
boundaries and zero flux (free) boundary
conditions on the top and bottom boundaries
In simulations, we use a grid with  nodes for
the one-dimensional problem and a  grid
for the two-dimensional case. We simulate with
 and set initial conditions
󰇛󰇜
󰇛󰇜 for
two-species and
󰇛󰇜
󰇛󰇜
󰇛󰇜 for three-
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
735
Volume 21, 2022
species competition models. As for diffusion rate,
we consider cases with regular diffusion 
(between 0.01 and 0.1 and has the same scale as
other parameters) and small diffusion .
To represent the result and compare the final
equilibrium state, we calculate the average solution
over computational domain for each species
󰇛󰇜󰇛󰇜
󰇛󰇜󰇛󰇜

where is the volume of domain . In most
results, we perform simulations till both
populations reach equilibrium, 󰇛󰇜󰇛󰇜
with  for each .
First, we present results for the two-species
interaction problem, where we simulate two sets of
parameters related to the species survival: Case 1
(one species survive) and Case 2 (both species
survive). We compare the dynamics of the average
solution and final state for three configurations: ,
󰇛󰇜, and 󰇛󰇜. Furthermore, the investigation
of the diffusion parameters scale to the final state is
presented. Next, we consider the three-species
competition model with three cases of the test
parameters: Case 1 (one species survive), Case 2
(two species survive), and Case 3 (all species
survive). Similar to the previous problems, we
investigate the final state and dynamic for ,
󰇛󰇜, and 󰇛󰇜 for small and regular diffusion.
After that, we present results for random
diffusion to investigate the effect of the species
diffusion coefficient on the final equilibrium state,
survival group, and time to reach an equilibrium
state. Next, we consider the influence of the initial
conditions on the equilibrium state and the time to
reach it.
3.1 Numerical Results for TwoSpecies
Competition Model
We consider a two-species competition model and
simulate two cases of the parameters:
Case 1 (one species survive)
󰇛󰇜
󰇛󰇜
󰇡  
  󰇢
Case 2 (both species survive)
󰇛󰇜
󰇛󰇜
󰇡  
  󰇢
As for diffusion rate, we set regular diffusion 
and ten times smaller diffusion .
First, we present the population density at the
final equilibrium when all the parameters are fixed.
In Fig. 2, we plot the solution for one and two-
dimensional formulations at the final time. Note that
for 2d, we plot the solution over the middle line
(). We observed that the effect of boundary
constraint (set at zero) is more severe in regular
diffusion groups, causing a lower final population
density than in small diffusion groups. With the
same survival status (one species survives, or both
species survive) and boundary conditions ( ,
󰇛󰇜, or 󰇛󰇜), the regular diffusion group
always arrives at a lower population density
equilibrium, compared to the small diffusion group.
In other words, when the diffusion is smaller, the
species can reach a higher population density
equilibrium.
Second, we present the dynamic of average
population density over time steps when all the
parameters are fixed. In Fig. 3, we present the
dynamic of the solution average over time for two
species systems. We observed that the time to reach
equilibrium is different among different diffusion
conditions (regular diffusion , or small diffusion
), survival status (one species survives, or both
species survive), and boundary conditions ( ,
󰇛󰇜, or 󰇛󰇜). In general,  and 󰇛󰇜 give
similar solutions, while 󰇛󰇜 gives a different
solution. In general, when the diffusion is small, it
usually takes more time to reach equilibrium.
Third, we present the population density over
󰇛󰇜 (imitation of the pond) and 󰇛󰇜 (imitation
of the river) at final equilibrium when all the
parameters are fixed. In Fig. 4 and Fig. 5, we
present the solution for 2d problems at the final time
in whole domain for regular diffusion 
and small diffusion , respectively. Like
Fig. 2, the effect of boundary constraint (set at zero)
is more severe in regular diffusion groups. If the
diffusion is regular (Fig. 4), the final population is
only dense in the middle, and the area of zero
population density is wider at the boundaries. If the
diffusion is small (Fig. 5), the final population is
more spread across the whole spatial domain.
To sum up, in a two-species competing model,
when we compare solutions of final population
density from two fixed groups of parameters
(diffusion rates/ diffusion rates*0.1 / growth rates/
competition rates) that leads to two survival status
(both species survive/ only one species survive), we
observed that species with diffusion rates that are 10
times smaller in scale than other parameters could
reach a higher equilibrium population density in the
whole spatial domain, and it usually takes the
species with small diffusion rates more time steps to
reach equilibrium.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
736
Volume 21, 2022
3.2 Numerical Results for Threespecies
Competition Model For the three-species competition model, we
consider three cases:

󰇛󰇜
󰇛󰇜
Fig. 2: Solution at the final time for regular diffusion ε  (solid line) and small diffusion ε  (dashed
line).

󰇛󰇜
󰇛󰇜
Fig. 3: Dynamic of the average solution over the domain for regular diffusion  (solid line) and small
diffusion  (dashed line).
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
737
Volume 21, 2022
󰇛󰇜
󰇛󰇜
Fig. 4: Solution at the final time for regular diffusion . First picture: first species. Second picture: second
species
󰇛󰇜
󰇛󰇜
Fig. 5: Solution at the final time for small diffusion . First picture: first species. Second picture:
second species
Case 1 (one species survive)
󰇛󰇜
󰇛󰇜
  
  
  
Case 2 (two species survive)
󰇛󰇜
󰇛󰇜
  
  
  
Case 3 (all three species survive)
󰇛󰇜
󰇛󰇜
  
  
  
First, we present the population density at the
final equilibrium when all the parameters are fixed
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
738
Volume 21, 2022
under the three-species competing model. In Fig. 6,
we plot the solution for one and two-dimensional
formulations at the final time. Similar to the
previous results, we plot the solution over the
middle line () for the 2d formulation. We
observe that the effect of boundary conditions is
generally more severe in regular diffusion groups
than in small diffusion groups. Like a two-species
competing system, in a three-species competing
system, when the diffusion rate is lower, species can
generally reach a higher final population density
both at the center of the spatial domain and near the
boundary. This shows that the boundary constraint
is lower for species with lower diffusion. More
investigation of the diffusion to the final state will
be presented in the next part of the results.
Then, we present the dynamic of average
population density over time steps when all the
parameters are fixed under a three-species
competing model. In Fig. 7, we present the dynamic
of the solution average over time for three species
system. Like a two-species competing system, in a
three-species competing system, generally, when the
diffusion is smaller, it usually takes more time to
reach equilibrium.
To sum up, in a three-species competing model,
when we compare solutions of final population
density from two fixed groups of parameters
(diffusion rates/ diffusion rates*0.1 / growth rates/
competition rates) that leads to two survival status
(both species survive/ only one species survive), we
observed that, similar to two-species competing
model, species with diffusion rates that are 10 times
smaller in scale than other parameters can reach a
higher equilibrium population density in whole
spatial domain with very few exceptions, and it
usually takes the species with small diffusion rates
more time steps to reach equilibrium.
3.3 Effect of the Diffusion
Next, we consider the influence of diffusion on the
equilibrium state. We set all other parameters
(growth rate "" and competition efficiency "") the
same and examine the diffusion rate. We perform
1000 simulations for each case with random
diffusion coefficients

We simulate with fixed initial conditions
󰇛󰇜
till both populations reach equilibrium, 󰇛󰇜
󰇛󰇜 with  for each .
First, we consider two-species competition
models in one and two-dimensional formulations.
Groups are represented by which species survive. In
the two-species competition model, we have the
following:
10 - first species survive,
01 - second species survived,
11 - both species survived and
00 - no one survived.
In Fig. 8 we plot scatter plots of the diffusion
rate of two species and colored survival status at the
final time for one and two-dimensional problems, as
well as corresponding time steps to reach
equilibrium. We observed that diffusion rate has a
huge impact on the final survival status of species
and the time step needed to reach equilibrium.
Different combination of diffusion rates of two
species leads to different final survival status. Under
 or 󰇛󰇜 spatial boundary condition, when the
diffusion rates of both species are above 0.07, no
species will survive. When one of the species'
diffusion rates is smaller than 0.07, the species with
a larger diffusion rate will survive. However, if two
species have approximately the same diffusion rates
and are both less than 0.07, both species will
survive. In comparison, under 󰇛󰇜 spatial
boundary condition, the threshold is reduced from
0.07 to around 0.04.
Different combination of diffusion rates of two
species also leads to different time steps needed to
reach equilibrium. When the combination of both
species' diffusion rates is on the border of varying
survival groups, it takes more time steps to reach
equilibrium. In comparison, it takes fewer time steps
to reach equilibrium when the combination of
diffusion rates of both species lies inside the
survival group.
Next, we consider three-species competition
models, where we have eight groups:
100 - first species survive,
010 - second species survived,
001 - third species survived,
011 - second and third survived,
101 - first and third survived,
110 - first and second survived,
111 - all species survived and
000 - no one survived.
In Fig. 9 we plot scatter plot of the diffusion rate
of three species and colored survival status at the
final time for one and two-dimensional spatial
boundary conditions, as well as corresponding time
steps to reach equilibrium. When we expand from
two species to three species in the system, we
observed that  or 󰇛󰇜 spatial boundary
condition still gives similar results, while 󰇛󰇜
gives a different result. In general, when the
diffusion rates of all three species are small, three
species can all survive. When all species have high
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
739
Volume 21, 2022
diffusion rates, no species can survive. When not all
three species have low diffusion rates, the species
with higher diffusion rates can survive. Similarly,
the different combination of diffusion rates of the
three species also leads to different time steps
needed to reach equilibrium and maximum located
on the border of varying survival groups.
3.4 Effect of the Initial Conditions
Previously we set the initial condition for both
species to be 0.5 and then got a general sense of the
final equilibrium difference. Next, we remove the
control for the initial condition ( 
󰇛󰇜
) and examine how changes in the initial
population can impact the final equilibrium
population density and time to reach it. With the
birth rate, competition rate, and diffusion rate fixed,
we change the initial population density for both
species.
In Fig. 10 we show a scatter plot of the initial
population density of two species and use a blue line
to represent a dynamic from initial condition points
(green color) to final equilibrium points of each
simulation for one and two-dimensional spatial
boundary conditions, as well as corresponding time
steps to reach equilibrium. We observed that
disregarding the value of the initial population of
both species, the final equilibrium population will
rest at a point of focus. The farther the combination
point of initial population density is from the final
equilibrium combination point, the more time steps
it would take to reach the equilibrium. We also
observed that the farther the initial population
density point is from the final equilibrium
population density point, the more time steps it
needs to reach the equilibrium.
In Fig. 11, we plot scatter plots of the initial
population density of three species and use a blue
line to connect initial condition points to the final
equilibrium points of each simulation and
corresponding time steps to reach equilibrium.
When we expand from two species to three species
in the system, we observed that  or 󰇛󰇜 spatial
boundary condition still gives a similar result. In
contrast, 󰇛󰇜 gives a different result. In general,
regardless of the values of the initial condition of
population density of three species, the final
equilibrium of population density for all three
species in the system will rest at a final focus point.
The farther a combination of the initial condition of
population density is from the final equilibrium
point, the more time steps it needs to take to reach
that equilibrium point.
4 Factor Analysis
Factor analysis method has long been applied to the
analysis of Population Dynamics, [47],[48] as well
as Ecosystem topics, [49],[50]. Finally, we present
factor analysis for the given model.
Previous literature about factor analysis suggests
that as sample size increases, the standard error in
factor loadings across repeated samples will
decrease, [51]. Therefore, we perform 10,000
simulations for each case with random parameters to
satisfy the sample size requirement.
We consider one and two-dimensional test
problems and simulate the system with random
coefficients , ,  with and initial
conditions
󰇛󰇜const:
Two-species:
󰇣
󰇛󰇜
󰇛󰇜󰇤
Three-species:

󰇣
󰇛󰇜
󰇛󰇜
󰇛󰇜󰇤.
We perform 10,000 simulations for each case
with random parameters

and 
󰇛󰇜. Two cases of the diffusion
coefficient are considered: regular diffusion (
) and small diffusion ().
In the investigation, we use the influence of the
following parameters to the system solution: ,
and for and  for two-species
competition and  for three species
competition.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
740
Volume 21, 2022

󰇛󰇜
󰇛󰇜
Case 1 (one species survive)
Case 2 (two species survive)
Case 3 (three species survive)
Fig. 6: Solution at the final time for regular diffusion ε (solid line) and small diffusion ε  (dashed
line). Red color: first species. Blue color: second species. Green color: third species

󰇛󰇜
󰇛󰇜
Case 1 (one species survive)
Case 2 (two species survive)
Case 3 (three species survive)
Fig. 7: Solution at the final time for regular diffusion ε (solid line) and small diffusion ε  (dashed
line). Red color: first species. Blue color: second species. Green color: third species
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
741
Volume 21, 2022

󰇛󰇜
󰇛󰇜
Case 1 (one species survive)
Case 2 (two species survive)
Fig. 8: Numerical results for random diffusion. Two-species model. First picture: 00 (grey), 01 (blue), 10 (red)
and 11 (green) groups. Second picture: number of time steps to reach the equilibrium

󰇛󰇜
󰇛󰇜
Fig. 9: Numerical results for random diffusion. Three-species model. First picture: 000 (light grey), 001 (grey),
010 (cyan), 011 (blue), 100 (salmon), 101 (red), 110 (lime) and 111 (green) groups. Second picture: number of
time steps to reach the equilibrium
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
742
Volume 21, 2022

󰇛󰇜
󰇛󰇜
Fig. 10: Numerical results for random initial conditions. Two-species model. First picture: the dynamic of the
solution average. Second picture: number of time steps to reach the equilibrium

󰇛󰇜
󰇛󰇜
Fig. 11: Numerical results for random initial conditions. Three-species model. First picture: the dynamic of the
solution average. Second picture: number of time steps to reach the equilibrium
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
743
Volume 21, 2022
4.1 Factor Analysis for TwoSpecies
Competition Model
To sum up the factor analysis for the final
population and time steps in a two-species
competing system, we observed that diffusion is
more dominant only when diffusion rates are on the
same scale with other parameters when diffusion
rates are ten times small in scale than other
parameters, growth rates/ competition rate/ initial
population density becomes the dominant factor. A
possible explanation for the change of the dominant
factor when the diffusion rate changes is the
existence of the boundary conditions (set at zero).
The boundary effect is more severe when the
diffusion rate is larger, making it a more critical
factor.
4.1.1 Factor Analysis for Final Population
Density in TwoSpecies Competition Model
In Fig. 12, we present correlation matrices as well as
loading matrix for each factor analysis for all three
cases , 󰇛󰇜 and 󰇛󰇜 with regular and small
diffusion, we only include all parameters and final
population density. From the correlation matrix, we
observed that under regular diffusion, final
population density is strongly correlated with both
diffusion rates and growth rates, while under small
diffusion, final population density is mostly strongly
correlated with growth rates.
In
Table 1, we present a summary for factor
analysis of final population density in a two-species
competition model. We observed that under regular
diffusion, diffusion rates are the most dominant
factors, while under small diffusion, growth rates
and competition rates become the most dominant
factors.
4.1.2 Factor Analysis for Time Steps to Reach
Equilibrium in TwoSpecies Competition Model
In Fig. 13, we present the correlation matrix and
loading matrix when we only include all parameters
and the number of time steps to reach equilibrium.
We observed that under both regular and small
diffusion, the time steps to reach equilibrium do not
have much correlation with all parameters. Under
regular diffusion, we only observed weak
correlations 󰇛󰇜 between time steps and
diffusion rates. Under small diffusion, we only
observed weak correlations 󰇛󰇜 between
time steps and growth rates.
In
Table 2, we present a summary for factor
analysis of the time steps until equilibrium in the
two-species competition model. We observed that
diffusion rates are essential factors only when
diffusion has the same scale as other parameters.
Under small diffusion, in which diffusion is ten
times smaller in scale than other parameters, the
dominant factors are growth rates or initial
population density.
4.2 Factor Analysis for ThreeSpecies
Competition Model
To sum up the factor analysis for the final
population and time steps in a three-species
competing system, we observed that these two-
factor analyses differ. For factor analysis for the
final population, we observed that diffusion is more
dominant only when diffusion rates are on the same
scale as other parameters. Growth rates are the
dominant factor when diffusion rates are ten times
smaller in scale than other parameters. For factor
analysis for time steps to reach equilibrium, we
observed that the importance of diffusion is similar,
disregarding the scale of diffusion rates. The
dominant factors are always competition rates/
growth rates/ initial population density.
4.2.1 Factor Analysis for Final Population
Density in ThreeSpecies Competition Model
In Fig. 14, we present three-species competition
system correlation matrices and a loading matrix for
each factor analysis for all three cases , 󰇛󰇜
and 󰇛󰇜 with regular and small diffusion. From
the correlation matrix, we observed that under
regular diffusion, final population density is
strongly correlated with both diffusion rates and
growth rates, while under small diffusion, final
population density is mostly strongly correlated with
growth rates and weakly correlated with diffusion
rates and competition rates.
In Table 3, we present a summary for factor
analysis of the final population in a three-species
competition model. When the diffusion rate is
regular, the dominant factors that cause variation are
diffusion rates; in contrast, when the diffusion rates
are low, the dominant factors become the growth
rates of the three species.
4.2.2 Factor Analysis for Time Steps to Reach
Equilibrium in ThreeSpecies Competition
Model
In Fig. 15, we present three-species competition
system correlation matrices and a loading matrix for
each factor analysis for all three cases , 󰇛󰇜
and 󰇛󰇜 with regular and small diffusion, we
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
744
Volume 21, 2022
only include the number of time steps to reach
equilibrium and all other parameters. We observed
that under both regular and small diffusion, the time
steps to reach equilibrium correlate very little with
all parameters. Under regular diffusion, we only
observed weak correlations 󰇛󰇜 between
time steps and diffusion rates. Under small
diffusion, we only observed weak
correlations 󰇛󰇜 between time steps and
growth rates.

󰇛󰇜
󰇛󰇜
Regular diffusion
Small diffusion
Fig. 12: Two-species competition model. The correlation matrix (first row) and loading (second row) only
include parameters and final population density
In
Table 4, we present a summary for factor analysis of the time steps to reach equilibrium in the three-species
competition model. We observed that the three-species competing system is different from the two-species
competing system in terms of dominant factors influencing time steps. In a two-species competing system,
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
745
Volume 21, 2022
when diffusion is regular, it is the dominant factor, and diffusion rates are unimportant when diffusion is small.
However, In a three-species competing system, the importance of diffusion is similar, disregarding the scale of
diffusion rates. The dominant factors are always competition rates/ growth rates/ initial population density.

󰇛󰇜
󰇛󰇜
Regular diffusion
Small diffusion
Fig. 13 Two-species competition model. Correlation matrix (first row) and loading (second row) only including
parameters and time steps to reach equilibrium
Table 1. Factor Analysis of the final population in two species system
Dimension
1D
2D(a)
2D(b)
1D
2D(a)
2D(b)
Diffusion
Regular
Regular
Regular
Small
Small
Small
Cum. Var.
66.83%
62.50%
66.43%
67.98%
68.33%
68.01%
Factor 1
Fin. Pop. 2 Diffusion 1
Fin. Pop. 2 Diffusion 1
Diffusion 2
Grow 1
Compete 2
Grow 1
Factor 2
Fin. Pop. 2 Diffusion 2
Fin. Pop. 2 Diffusion 2
Diffusion 1
Grow 2
Compete 1
Grow 2
Factor 3
Grow 2
Compete 2
Grow 2
Compete 1
Diffusion 2
Compete 1
Factor 4
Grow 1
Grow 1
Grow 1
Fin. Pop. Compete 2
Fin. Pop. Diffusion 1
Fin. Pop. 2 Compete 2
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
746
Volume 21, 2022
Factor 5
Compete 2
Fin. Pop. 2 Grow 2
Fin. Pop. Compete 1
Diffusion 1
Diffusion 1
Diffusion 1
Table 2. Factor Analysis of time steps in two species system
Dimension
1D
2D(a)
2D(b)
1D
2D(a)
2D(b)
Diffusion
Regular
Regular
Regular
Small
Small
Small
Cum. Var.
57.54%
61.85%
57.91%
60.21%
59.99%
60.77%
Factor 1
Diffusion 1
Diffusion 1
Time step n Diffusion 2
Init. Pop. 1
Grow 1
Grow 1
Factor 2
Diffusion 2
Time step n Diffusion 2
Compete 1
Grow 1
Time step n Grow 2
Init. Pop. 1
Factor 3
Compete 2
Grow 2
Time step n Grow 1
Time step n Grow 1
Compete 2
Compete 2
Factor 4
Init. Pop. 2
Compete 2
Init. Pop. 2
Grow 2
Init. Pop. 1
Time step n Grow 2
Factor 5
Time step n Diffusion 2
Time step n Grow 1
Diffusion 1
Compete 2
Time step n Grow 2
Time step n Grow 2
Table 3. Factor Analysis of the final population in three species system
Dimension
1D
2D(a)
2D(b)
1D
2D(a)
2D(b)
Diffusion
Regular
Regular
Regular
Small
Small
Small
Cum. Var.
58.45%
60.40%
58.59%
59.10%
58.92%
59.10%
Factor 1
Init. Pop. 3
Diffusion 1
Compete 3
Grow 2
Grow 2
Init. Pop. 3
Factor 2
Grow 2
Compete 3
Init. Pop. 3
Compete 3
Init. Pop. 3
Compete 3
Factor 3
Init. Pop. 1
Init. Pop. 2
Grow 2
Diffusion 3
Compete 3
Init. Pop. 1
Factor 4
Compete 3
Grow 3
Init. Pop. 1
Init. Pop. 1
Init. Pop. 1
Time step n Grow 3
Factor 5
Diffusion 1
Diffusion 3
Diffusion 1
Time step n Grow 3
Time step n Grow 3
Grow 2
Factor 6
Time step n Diffusion 3
Diffusion 2
Time step n Diffusion 3
Init. Pop. 3
Diffusion 3
Diffusion 3
Factor 7
Compete 2
Compete 2
Compete 2
Time step n Compete 2
Compete 1
Compete 1
Factor 8
Diffusion 2
Compete 2
Diffusion 2
Compete 1
Compete 2
Time step n Compete 2
Factor 9
Compete 1
Grow 2
Compete 1
Diffusion 2
Compete 1
Compete 1
Table 4. Factor Analysis of time steps in three species system
Dimension
1D
2D(a)
2D(b)
1D
2D(a)
2D(b)
Diffusion
Regular
Regular
Regular
Small
Small
Small
Cum. Var.
58.45%
60.40%
58.59%
59.10%
58.92%
59.10%
Factor 1
Init. Pop. 3
Diffusion 1
Compete 3
Grow 2
Grow 2
Grow 2
Factor 2
Grow 2
Compete 3
Init. Pop. 3
Init. Pop. 3
Init. Pop. 3
Init. Pop. 3
Factor 3
Init. Pop. 1
Diffusion 1
Grow 2
Compete 3
Compete 3
Compete 3
Factor 4
Compete 3
Init. Pop. 2
Init. Pop. 1
Diffusion 3
Init. Pop. 1
Init. Pop. 1
Factor 5
Diffusion 1
Grow 3
Diffusion 1
Init. Pop. 1
Time step n Grow 3
Time step n Grow 3
Factor 6
Time step n Diffusion 3
Diffusion 2
Time step n Diffusion 3
Time step n Grow 3
Diffusion 3
Diffusion 3
Factor 7
Compete 2
Compete 2
Compete 2
Time step n Compete 2
Compete 1
Time step n Compete 2
Factor 8
Diffusion 2
Compete 2
Diffusion 2
Compete 1
Compete 2
Compete 1
Factor 9
Compete 1
Grow 2
Compete 1
Diffusion 2
Compete 1
Compete 1
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
747
Volume 21, 2022

󰇛󰇜
󰇛󰇜
Regular diffusion
Small diffusion
Fig. 14: Three-species competition model. Correlation matrix (first row) and loading (second row) only
including all parameters and final population density
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
748
Volume 21, 2022

󰇛󰇜
󰇛󰇜
Regular diffusion
Small diffusion
Fig. 15: Three-species competition model. Correlation matrix (first row) and loading (second row) only
including all parameters and time steps to reach equilibrium
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
749
Volume 21, 2022
4.3 Statistical Summary for Survival Status
in All Simulations
For a two-species competing system, we released
control for all variables, including diffusion rates,
growth and competition efficiency, and initial
population density, and ran 10,000 simulations for
all , 󰇛󰇜, and 󰇛󰇜 spatial boundary
conditions, and summarized the proportion of
different survival status groups, as is in
Table 5. We observed that  and 󰇛󰇜 spatial
boundary conditions give the similar proportion of
survival groups, while 󰇛󰇜 gives more different
proportions.
When there is regular diffusion, for  and
󰇛󰇜 spatial boundary conditions, there is around
 of simulations that result in either 00(no
species survives), 01 (species two survives), or 10
(species one survives). n comparison, only around
 of simulations result in both species surviving
together. or 󰇛󰇜 spatial boundary conditions,
 of simulations result in 00 (both species go
extinct), while around  and  result in
01 (species two survives) or 10 (species one
survives), and very rare, about  of the
simulation result in 11 (both species survive). To
conclude, when there is regular diffusion, it's
generally harder to arrive at a final solution where
both species survive; it's even harder, almost
impossible, for both species to survive when we are
approximating a pond (󰇛󰇜 scenario).
When there is small diffusion, all three spatial
boundary conditions (, 󰇛󰇜, and 󰇛󰇜) give a
similar proportion of survival status. t is most
likely that both species survive together (around
 of simulations). However, it is almost
impossible that both species will go extinct.
Furthermore, it is equally likely that only one
species survives (the proportions of simulation are
both around ).
For the three-species competing system, we
released control for all variables, including diffusion
rates, growth and competition efficiency, and initial
population density, ran 10,000 simulations for all
, 󰇛󰇜, and 󰇛󰇜 spatial boundary conditions,
and summarized the proportion of different survival
status groups, as is in Σφάλμα! Το αρχείο
προέλευσης της αναφοράς δεν βρέθηκε.. we
observed that  and 󰇛󰇜 spatial boundary
conditions give the similar proportion of survival
groups, while 󰇛󰇜 shows more different
proportion. When there is regular diffusion, the
probability of one species surviving (001, 010, and
100) is similar to one another (around  for 
or 󰇛󰇜 scenario, and approximately  for
󰇛󰇜 scenario), while the probability of two
species survive (011, 110 and 101) is similar to one
another(around  for  or 󰇛󰇜 scenario, and
about  for 󰇛󰇜 scenario). In general, when
there are three species, it is almost impossible for all
three species to survive together, and no species go
extinct; what's worse, for 󰇛󰇜 scenario, which
approximates the lake, around  of the
simulation result in 000 (no species survive). When
there is small diffusion, all three spatial boundary
conditions (, 󰇛󰇜, and 󰇛󰇜) give a similar
proportion of survival status. It's equally likely that
two or all species survive (the proportion for
survival group 011/110/101/111 are all around
). It's equally likely that only one species
survives (the proportion for survival group
001/010/100 are all around ), and it's almost
impossible that all three species will go extinct.
5 Discussion
Numerical results with fixed values of
parameters that cause different survival
status
Species with small diffusion rates (10 times smaller
in scale than other parameters) can reach a higher
equilibrium population density in the whole spatial
domain with very few exceptions. It usually takes
the species with small diffusion rates more time
steps to reach equilibrium.
Effect of diffusion examined by
numerically releasing the control for
diffusion rates
Different combinations of diffusion rates align with
varying statuses of survival. Groups of the
combination of diffusion rates in scatter plot overlap
with groups of varying survival statuses. This
overlap can be interpreted as diffusion rates strongly
affecting the survival status. In a two-species
competing system, under  or 󰇛󰇜 spatial
boundary condition, when the diffusions of both
species are above 0.07, no species will survive;
under 󰇛󰇜 spatial boundary condition, the
threshold is reduced from 0.07 to around 0.04. In
three-species competing species, we observed a
similar effect of diffusion on survival status as in a
two-species system. Furthermore, species with a
combination of diffusion rates on the border of
different survival groups take more time steps to
reach equilibrium than those whose combination is
inside the groups.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
750
Volume 21, 2022
Effect of the initial condition examined by
numerically releasing the control for
initial population density
For both two-species competing model and the
three-species competing model, regardless of the
values of the initial condition of population density
of species, the final equilibrium of population
density for all three species in the system will rest at
a final focus point. The farther a combination of the
initial condition of population density is from the
final equilibrium point, the more time steps it takes
to reach that equilibrium point.
Table 5. Survival Status for two species system
00
01
10
11
Regular diffusion
1d
2796
27.96%
3059
30.59%
3064
30.64%
1081
10.81%
2d(a)
6888
68.88%
1487
14.87%
1489
14.89%
136
1.36%
2d(b)
3125
31.25%
2904
29.04%
3043
30.43%
928
9.28%
Small diffusion
1d
0
0.00%
2251
22.51%
2279
22.79%
5470
54.70%
2d(a)
26
0.26%
2522
25.22%
2533
25.33%
4919
49.19%
2d(b)
0
0.00%
2184
21.84%
2192
21.92%
5624
56.24%
Table 6. Survival Status for three species system
Group
000
001
010
100
011
110
101
111
Regular diffusion
1d
1442
14.42%
2026
20.26%
2031
20.31%
1987
19.87%
806
8.06%
822
8.22%
730
7.30%
156
1.56%
2d(a)
5800
58.00%
1262
12.62%
1319
13.19%
1234
12.34%
131
1.31%
130
1.30%
118
1.18%
6
0.06%
2d(b)
1754
17.54%
2012
20.12%
2034
20.34%
1985
19.85%
712
7.12%
728
7.28%
651
6.51%
124
1.24%
Small diffusion
1d
0
0.00%
639
6.39%
648
6.48%
556
5.56%
2031
20.31%
2023
20.23%
1943
19.43%
2160
21.60%
2d(a)
3
0.03%
858
8.58%
835
8.35%
773
7.73%
1957
19.57%
2001
20.01%
1863
18.63%
1710
17.10%
2d(b)
0
0.00%
651
6.51%
655
6.55%
566
5.66%
2030
20.30%
2023
20.23%
1941
19.41%
2134
21.34%
Factor Analysis after releasing control for
all parameters
In a two-species competing model, diffusion rates
are the more dominant factor for both final
population density and time to reach equilibrium
only when diffusion rates are on the same scale as
other parameters. When diffusion rates are ten times
small in scale than other parameters, growth rates/
competition rate or initial population density
becomes the dominant factor. In a three-species
competing system, for factor analysis for the final
population, we observed that diffusion is more
dominant only when diffusion rates are in the same
scale with other parameters, when diffusion rates are
ten times small in scale than other parameters,
growth rates are the dominant factor; As of factor
analysis for time steps to reach equilibrium, the
importance of diffusion is similar disregard the scale
of diffusion rates. The dominant factors are always
competition rates/ growth rates/ initial population
density.
6 Conclusion
A mathematical model of the spatial-temporal multi-
species competition is considered. A discrete system
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
751
Volume 21, 2022
is constructed using a finite volume approximation
with a semi-implicit time approximation. The
numerical results for two- and three-species models
are presented for several exceptional cases of the
parameters related to the survival status. We
considered the one and two-dimensional model
problems with two cases of the boundary conditions
for the two-dimensional case. First, the effect of
diffusion is investigated numerically. In special
cases where parameters are fixed, we observed that
the impact of boundary constraint is more severe in
regular diffusion groups than in small diffusion
groups, causing a lower population density both in
the middle and near the boundary of the domain. We
also observed that the dynamic of  is similar to
󰇛󰇜, while 󰇛󰇜 gives different dynamic. This
suggests that we can use  to approximate 󰇛󰇜
to save computation time. Furthermore, from
general cases where we release holds of diffusion
while keeping other parameters fixed, we observed
that different combinations of species' diffusion
rates in the system lead to different final survival
statuses of species. When the combination of both
species' diffusion rates is on the border of varying
survival groups, it takes more time for these groups
of species to reach equilibrium. Second, the effect of
the initial condition is investigated numerically. We
release holds of initial conditions while keeping
other parameters fixed. We observed similar
patterns for both two-species competing systems
and three-species competing systems. Take a two-
species competing system as an example; we
observed that disregarding the value of the initial
population of both species, the final equilibrium
population will rest at a point of focus. The farther
the combination point of initial population density is
from the final equilibrium combination point, the
more time steps it would take to reach the
equilibrium.
Finally, the impact of parameters on the system
stability is considered by simulating the spatial-
temporal model with random input parameters.
Factor analysis and statistical summary of the
survival status of species in the system were
performed. We observed that diffusion rates are the
dominant factor. When diffusion rates are regular
and on the same scale as other parameters. In
contrast, when diffusion rates are small, which are
ten times smaller in the scale of other parameters,
growth and competition rates become the dominant
factors. In a statistical summary of species survival
status, we observed a similar pattern for both two-
species competing systems and three-species
competing systems. For both systems, in each
spatial boundary condition ( , 󰇛󰇜, and
󰇛󰇜), when the number of survived species is the
same in the system (001, 010, 100 survival group in
which only one species survive), the proportion of
simulation is similar. However, the proportion of
simulation for different survival groups varies when
the diffusion rate is in different scale, this is the case
for both competing systems. Take a two-species
competing system as an example. When there is
small diffusion, all three spatial boundary conditions
(, 󰇛󰇜, and 󰇛󰇜) give a similar proportion
of survival status. It is most likely that both species
survive together ( -  of simulation), it is
almost impossible that both species will go extinct
( - ), and it is equally likely that only one
species survives (both 01 and 10 survival groups
take  -  of simulation). In contrast, when
there is regular diffusion in two species competing
systems, we observed that  and 󰇛󰇜 spatial
boundary conditions give a different proportion of
simulation for survival groups from 󰇛󰇜. Take a
two-species competing system as an example; when
there is regular diffusion, it is generally harder to
arrive at a final solution where both species survive
(around  of simulation); it is even harder,
almost impossible ( of simulation), for both
species to survive when we are approximating a
pond (󰇛󰇜 scenario).
In future works, we will concentrate on the
following directions: (1) a real-world range of
parameters for a given application; (2) apply the
result of factor analysis to build machine learning
predictive models, and predict final population
density, as well as time to reach equilibrium; and (3)
use the model to predict the equilibrium population
density and simulate future catastrophic event
(given the assumption that the system is at
equilibrium then catastrophe strokes).
References:
[1] A. Okubo and S. A. Levin, Diffusion and
Ecological Problems: Modern Perspectives,
vol. 14. Springer, 2001. Accessed: Jul. 30,
2022. [Online]. Available:
https://link.springer.com/book/10.1007/978-1-
4757-4978-6
[2] J. D. Murray, Mathematical biology: I. An
introduction. Springer, 2002.
[3] G. I. Marchuk, Mathematical models in
environmental problems. Elsevier, 2011.
[4] A. Marasco, A. Picucci, and A. Romano,
“Market share dynamics using Lotka–Volterra
models,” Technological forecasting and
social change, vol. 105, pp. 49–62, 2016.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
752
Volume 21, 2022
[5] W. Zhang and J. S. L. Lam, “Maritime cluster
evolution based on symbiosis theory and
Lotka–Volterra model,” Maritime Policy &
Management, vol. 40, no. 2, pp. 161–176,
2013.
[6] W. Windarto and E. Eridani, “On
modification and application of Lotka-
Volterra competition model,” in Aip
conference proceedings, 2020, vol. 2268, p.
050007.
[7] S.-Y. Wang, W.-M. Chen, and X.-L. Wu,
“Competition analysis on industry populations
based on a three-dimensional lotka–volterra
model,” Discrete Dynamics in Nature and
Society, vol. 2021, 2021.
[8] M. A. Khan, M. Azizah, S. Ullah, and others,
“A fractional model for the dynamics of
competition between commercial and rural
banks in Indonesia,” Chaos, Solitons &
Fractals, vol. 122, pp. 32–46, 2019.
[9] P. A. Montagna, A. L. Sadovski, S. A. King,
K. K. Nelson, T. A. Palmer, and K. H.
Dunton, “Modeling the effect of water level
on the Nueces Delta marsh community,”
Wetlands Ecol Manage, vol. 25, no. 6, pp.
731–742, Dec. 2017, doi: 10.1007/s11273-
017-9547-x.
[10] Q. Chen, R. Han, F. Ye, and W. Li, “Spatio-
temporal ecological models,” Ecological
Informatics, vol. 6, no. 1, pp. 37–43, Jan.
2011, doi: 10.1016/j.ecoinf.2010.07.006.
[11] Y. R. Zelnik, J.-F. Arnoldi, and M. Loreau,
“The Impact of Spatial and Temporal
Dimensions of Disturbances on Ecosystem
Stability,” Frontiers in Ecology and
Evolution, vol. 6, 2018, Accessed: Jul. 30,
2022. [Online]. Available:
https://www.frontiersin.org/articles/10.3389/f
evo.2018.00224
[12] A. Alhasanat and C. Ou, “Minimal-speed
selection of traveling waves to the Lotka
Volterra competition model,” Journal of
Differential Equations, vol. 266, no. 11, pp.
7357–7378, May 2019, doi:
10.1016/j.jde.2018.12.003.
[13] M. K. A. Gavina et al., “Multi-species
coexistence in Lotka-Volterra competitive
systems with crowding effects,” Sci Rep, vol.
8, no. 1, p. 1198, Dec. 2018, doi:
10.1038/s41598-017-19044-9.
[14] P. Zhou, “On a Lotka-Volterra competition
system: diffusion vs advection,” Calc. Var.,
vol. 55, no. 6, p. 137, Oct. 2016, doi:
10.1007/s00526-016-1082-8.
[15] X.-Q. Zhao and P. Zhou, “On a Lotka–
Volterra competition model: the effects of
advection and spatial variation,” Calc. Var.,
vol. 55, no. 4, p. 73, Jun. 2016, doi:
10.1007/s00526-016-1021-8.
[16] V. Dakos, “Identifying best-indicator species
for abrupt transitions in multispecies
communities,” Ecological Indicators, vol. 94,
pp. 494–502, Nov. 2018, doi:
10.1016/j.ecolind.2017.10.024.
[17] T. Tahara et al., “Asymptotic stability of a
modified Lotka-Volterra model with small
immigrations,” Sci Rep, vol. 8, no. 1, Art. no.
1, May 2018, doi: 10.1038/s41598-018-
25436-2.
[18] M. Benaïm and C. Lobry, “Lotka–Volterra
with randomly fluctuating environments or
‘how switching between beneficial
environments can make survival harder,’” The
Annals of Applied Probability, vol. 26, no. 6,
pp. 3754–3785, Dec. 2016, doi: 10.1214/16-
AAP1192.
[19] M. Liu and M. Fan, “Permanence of
Stochastic Lotka–Volterra Systems,” J
Nonlinear Sci, vol. 27, no. 2, pp. 425–452,
Apr. 2017, doi: 10.1007/s00332-016-9337-2.
[20] S. Kumar, R. Kumar, R. P. Agarwal, and B.
Samet, “A study of fractional Lotka-Volterra
population model using Haar wavelet and
Adams-Bashforth-Moulton methods,”
Mathematical Methods in the Applied
Sciences, vol. 43, no. 8, pp. 5564–5578, 2020,
doi: 10.1002/mma.6297.
[21] K. Devarajan, T. L. Morelli, and S. Tenan,
“Multispecies occupancy models: review,
roadmap, and recommendations,” Ecography,
vol. 43, no. 11, pp. 1612–1624, Nov. 2020,
doi: 10.1111/ecog.04957.
[22] F. Boschetti et al., “Setting priorities for
conservation at the interface between ocean
circulation, connectivity, and population
dynamics,” Ecol Appl, vol. 30, no. 1, Jan.
2020, doi: 10.1002/eap.2011.
[23] S.-S. Baek, Y. S. Kwon, J. Pyo, J. Choi, Y. O.
Kim, and K. H. Cho, “Identification of
influencing factors of A. catenella bloom
using machine learning and numerical
simulation,” Harmful Algae, vol. 103, p.
102007, Mar. 2021, doi:
10.1016/j.hal.2021.102007.
[24] B. C. T. Cabella, A. S. Martinez, and F.
Ribeiro, “Full analytical solution and
complete phase diagram analysis of the
Verhulst-like two-species population
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
753
Volume 21, 2022
dynamics model,” arXiv preprint
arXiv:1010.3361, 2010.
[25] K. Murty and D. Rao, “Approximate
analytical solutions of general Lotka-Volterra
equations,” Journal of mathematical analysis
and applications, vol. 122, no. 2, pp. 582–
588, 1987.
[26] R. BHARDWAJ and S. DAS,
“SYNCHRONIZATION OF CHAOTIC
FOOD CHAIN WITH COMPETITIVE
SPECIES,” Bull. Cal. Math. Soc, vol. 111, no.
1, pp. 53–64, 2019.
[27] J. O. Ojonubah and M. H. Mohd, “Impacts of
asymmetric biotic interactions and
environmental factors on the presence-
absence of multispecies.,” Pertanika Journal
of Science & Technology, vol. 28, no. 1, 2020.
[28] Z. R. Miller, P. Lechón-Alonso, and S.
Allesina, “No robust multispecies coexistence
in a canonical model of plant–soil feedbacks,”
Ecology Letters, 2022.
[29] P. S. Jamwal, M. Di Febbraro, M. L.
Carranza, M. Savage, and A. Loy, “Global
change on the roof of the world: Vulnerability
of Himalayan otter species to land use and
climate alterations,” Diversity and
Distributions, vol. 28, no. 8, pp. 1635–1649,
2022.
[30] R. Engler, A. Guisan, and L. Rechsteiner, “An
improved approach for predicting the
distribution of rare and endangered species
from occurrence and pseudo-absence data,”
Journal of applied ecology, vol. 41, no. 2, pp.
263–274, 2004.
[31] L. Wang, “Study on asymptotic behavior of
stochastic Lotka–Volterra system in a polluted
environment,” Advances in Difference
Equations, vol. 2021, no. 1, pp. 1–18, 2021.
[32] Z. Luo and X. Fan, “Optimal control for an
age-dependent competitive species model in a
polluted environment,” Applied Mathematics
and Computation, vol. 228, pp. 91–101, 2014.
[33] F. Vadillo, “Comparing stochastic Lotka
Volterra predator-prey models,” Applied
Mathematics and Computation, vol. 360, pp.
181–189, 2019.
[34] M. Liu, K. Wang, and Q. Wu, “Survival
analysis of stochastic competitive models in a
polluted environment and stochastic
competitive exclusion principle,” Bulletin of
mathematical biology, vol. 73, no. 9, pp.
1969–2012, 2011.
[35] R. Lou, Z. Lv, S. Dang, T. Su, and X. Li,
“Application of machine learning in ocean
data,” Multimedia Systems, Feb. 2021, doi:
10.1007/s00530-020-00733-x.
[36] C. A. Johnson et al., “Science to inform
policy: linking population dynamics to habitat
for a threatened species in Canada,” Journal
of Applied Ecology, vol. 57, no. 7, pp. 1314–
1327, 2020.
[37] H.-C. Thorsen-Meyer et al., “Dynamic and
explainable machine learning prediction of
mortality in patients in the intensive care unit:
a retrospective study of high-frequency data in
electronic patient records,” The Lancet Digital
Health, vol. 2, no. 4, pp. e179–e191, 2020.
[38] M. R. Keshtkaran et al., “A large-scale neural
network training framework for generalized
estimation of single-trial population
dynamics,” BioRxiv, pp. 2021–01, 2022.
[39] Z. P. Chairez, “Spatial-temporal models of
multi-species interaction to study impacts of
catastrophic events,” Texas A&M University-
Corpus Christi, 2020.
[40] M. Vasilyeva, A. Sadovski, and D.
Palaniappan, “Multiscale solver for multi-
component reaction-diffusion systems in
heterogeneous media,” arXiv preprint
arXiv:2209.04495, 2022.
[41] A. A. Samarskii, The theory of difference
schemes. CRC Press, 2001.
[42] A. A. Samarskii and P. N. Vabishchevich,
Computational heat transfer. 1995.
[43] P. Vabishchevich, “Additive schemes
(splitting schemes) for some systems of
evolutionary equations,” Mathematics of
Computation, vol. 83, no. 290, pp. 2787–
2797, 2014.
[44] N. Afanasyeva, P. N. Vabishchevich, and M.
Vasilyeva, “Unconditionally stable schemes
for non-stationary convection-diffusion
equations,” in International conference on
numerical analysis and its applications, 2012,
pp. 151–157.
[45] G. Strang, Linear algebra and its
applications. Belmont, CA: Thomson,
Brooks/Cole, 2006.
[46] C. R. Harris et al., “Array programming with
NumPy,” Nature, vol. 585, no. 7825, pp. 357–
362, 2020.
[47] G. C. Varley and G. R. Gradwell, “Recent
Advances in Insect Population Dynamics,”
Annual Review of Entomology, vol. 15, no. 1,
pp. 1–24, 1970, doi:
10.1146/annurev.en.15.010170.000245.
[48] S. D. Albon, T. N. Coulson, D. Brown, F. E.
Guinness, J. M. Pemberton, and T. H.
Clutton-Brock, “Temporal changes in key
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
754
Volume 21, 2022
factors and key age groups influencing the
population dynamics of female red deer,”
Journal of Animal Ecology, vol. 69, no. 6, pp.
1099–1110, 2000, doi: 10.1111/j.1365-
2656.2000.00485.x.
[49] N. Bernier and F. Gillet, “Structural
relationships among vegetation, soil fauna and
humus form in a subalpine forest ecosystem: a
Hierarchical Multiple Factor Analysis
(HMFA),” Pedobiologia, vol. 55, no. 6, pp.
321–334, Nov. 2012, doi:
10.1016/j.pedobi.2012.06.004.
[50] P. Petitgas et al., “Ecosystem spatial structure
revealed by integrated survey data,” Progress
in Oceanography, vol. 166, pp. 189–198, Sep.
2018, doi: 10.1016/j.pocean.2017.09.012.
[51] R. C. MacCallum, K. F. Widaman, S. Zhang,
and S. Hong, “Sample size in factor analysis,”
Psychological Methods, vol. 4, no. 1, pp. 84–
99, 1999, doi: 10.1037/1082-989X.4.1.84.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.85
Maria Vasilyeva, Youwen Wang,
Sergei Stepanov, Alexey Sadovski
E-ISSN: 2224-2880
755
Volume 21, 2022