Numerical investigation and factor analysis of the spatial-temporal multi-species competition problem

Understanding the stability of ecosystems is of fundamental importance to ecology [21, 20, 18]. 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 [17, 31, 29, 28, 13]. For example, the model can be used to simulate the marsh ecosystems for the wetlands at the Nueces River mouth [19].

The LVC model is based on the logistic population model

where $u$ is the size of the population at a given time $t$and $r>0$ is per-capita growth rate.

For the general case of multispecies competition, we have

where $u^{(k)}$ is the population of the $k$-th species, $r_{k}$ is the $k$-th population growth rate, $\alpha_{kl}$ is the interaction coefficient due to competition ($u^{(l)}$ compete with $u^{(k)}$), and $M$ 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 important and has a large impact on the final equilibrium state [9, 30]. 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 multispecies interaction, we have

where $\varepsilon_{k}$ is the diffusion coefficient. This model incorporates spatial structure by adding diffusion terms in equations and considering the system of equations in domain $\Omega\subset R^{d}$, where $d$ is the spatial dimension.

Previous research has put much effort into modifying the Lotka–Volterra competition model (LVC). The parameters $r_{k}$, $\alpha_{kl}$, and $\varepsilon_{k}$ can be modeled as functions instead of constant. Diffusion $\varepsilon_{k}$, 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 [3]. The crowding effect of diffusion $\varepsilon_{k}$ can also be modeled within LVC model [12]. Advection rate can be introduced to the LVC model as a supplement of diffusion to bring much richer phenomena [33, 32]. Carrying capacity can be replaced by function instead of constant and can vary between species; connection topology can be modeled as ‘loop’, ‘star’, ‘chain’, or ‘full’ connection when there are more than two species in the system [10]. 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 [25]. The random fluctuating environment can also be modeled within the LV model to show how switching between environments can make survival harder [5]. 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 [15].

Multiple methods are used for solving the LVC model, such as Haar wavelet (HW), Adams-Bashforth-Moulton (ABM) methods [14], 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 [11]. For example, the population dynamics model of individual reefs can be compared with data of coral reefs in Pilbara [7]. Numerical models and machine learning can be combined to identify the factors that influence Alexandrium catenella blooms [4].

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. After understanding the effect of parameters and initial conditions on the test problems, we evaluate the impact of parameters on the system stability by simulating the spatial-temporal model with random input parameters and performing factor analysis. Simulations with different combinations of parameters and initial conditions support the hypothesis that such systems reach equilibrium sooner or later. 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.

The paper is organized as follows. Section 2 describes the mathematical model with fine-scale approximation using the finite volume method and a semi-implicit scheme for time approximation. In Section 3, we present 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.

The mathematical model is described by a coupled system of nonlinear diffusion - reaction equations in domain $\Omega\subset R^{d}$ ($d=1,2$)

where $k=1,..,M$, where $M$ is the number of species (equations). Here $u^{(k)}$ is the population of the $k$-th species, $r_{k}$ is the $k$-th population reproductive growth rate, $\alpha_{kl}$ is the interaction coefficient due to competition ($u^{(l)}$ compete with $u^{(k)}$) and $\varepsilon_{k}$ is the diffusion coefficient c[8].

The system of equation is considered with initial conditions

and boundary conditions

In this work, we consider following special cases:

• •

  Two-species competition

  $$\begin{split}\frac{\partial u^{(1)}}{\partial t}-\nabla\cdot(\varepsilon_{1}\nabla u^{(1)})&=r_{1}u^{(1)}(1-u^{(1)})-\alpha_{12}u^{(1)}u^{(2)},\\ \frac{\partial u^{(2)}}{\partial t}-\nabla\cdot(\varepsilon_{2}\nabla u^{(2)})&=r_{2}u^{(2)}(1-u^{(2)})-\alpha_{21}u^{(1)}u^{(2)},\end{split}$$

• •

  Three-species competition

  $$\begin{split}\frac{\partial u^{(1)}}{\partial t}-\nabla\cdot(\varepsilon_{1}\nabla u^{(1)})&=r_{1}u^{(1)}(1-u^{(1)})-\alpha_{12}u^{(1)}u^{(2)}-\alpha_{13}u^{(1)}u^{(3)},\\ \frac{\partial u^{(2)}}{\partial t}-\nabla\cdot(\varepsilon_{2}\nabla u^{(2)})&=r_{2}u^{(2)}(1-u^{(2)})-\alpha_{21}u^{(1)}u^{(2)}-\alpha_{23}u^{(3)}u^{(2)},\\ \frac{\partial u^{(3)}}{\partial t}-\nabla\cdot(\varepsilon_{3}\nabla u^{(3)})&=r_{3}u^{(3)}(1-u^{(3)})-\alpha_{31}u^{(1)}u^{(3)}-\alpha_{32}u^{(2)}u^{(3)},\end{split}$$

For numerical solution of the problem (1) with initial and boundary conditions (2)-(3), we construct the structured grid for domain $\Omega=[0,L]^{d}$ ($d=1,2$) with $h=L/(N-1)$, where $N$ is the number of nodes in each direction [24]. We set

where $|K_{i}|$ is cell volume and $u_{i}^{(k)}$ is the cell average value of the function $u^{(k)}$ on cell $K_{i}$.

For the finite volume method, we have

with

where $d_{ij}$ is the distance between to cell center points $x_{i}$ and $x_{j}$, $|e_{ij}|$ is the length of the interface between two cells $K_{i}$ and $K_{j}$. Note that, for structured uniform grid we have $d_{ij}=h$, $|e_{ij}|=1$ for one-dimensional case and $|e_{ij}|=h$ for two-dimensional problem.

For the time approximation, we use a semi-implicit scheme [23, 26, 1]

where $\tau$ is the time step size and $\check{u}^{k}_{i}$ is the solution from previous time layer. Note that, we obtain uncoupled system of equations and can solve equation for each component separately.

Before we make factor analysis for randomly generated parameters, we consider some special cases with parameters (growth rate ”r”, competition efficiency ”$\alpha$”, and diffusion rate ”$\varepsilon$” ) 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 multispecies competition in domain $\Omega$. We consider the following domain and boundary conditions configurations:

• •

  1D: $\Omega=[0,1]$ with zero (fixed) boundary conditions on $\partial\Omega$.

• •

  2D(a): $\Omega=[0,1]^{2}$ with zero (fixed) boundary conditions on $\partial\Omega$.

• •

  2D(b): $\Omega=[0,1]^{2}$ with zero (fixed) boundary conditions on left and right boundaries, and zero flux (free) boundary conditions on top and bottom boundaries

In simulations, we use grid with $100$ nodes for one-dimensional problem and $25\times 25$ grid for two-dimensional case. We simulate with $\tau=1$ and set initial conditions $u_{0}^{(1)}=u_{0}^{(2)}=0.5$ for two-species and $u_{0}^{(1)}=u_{0}^{(2)}=u_{0}^{(3)}=0.5$ for three-species competition models. As of diffusion rate, we consider cases with regular diffusion $\varepsilon=D$ (which is between 0.01 and 0.1, and has same scale as other parameters) and small diffusion $\varepsilon=D/10$.

To represent result and compare final equilibrium state, we calculate average solution over computational domain $\Omega$ for each species

where $|\Omega|$ is the volume of domain $\Omega$. In most results, we perform simulations till both population reach equilibrium, $|\bar{u}^{(k)}-\check{\bar{u}}^{(k)}|<\epsilon$ with $\epsilon=10^{-5}$ for each $k$.

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: 1D, 2D(a), and 2D(b). 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). Similarly to the previous problems, we investigate the final state and dynamic for 1D, 2D(a), and 2D (b) 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.

Factor analysis method has long been applied to analysis of Population Dynamic [27, 2] as well as Ecosystem topic [6, 22]. 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 [16]. 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 $\varepsilon_{k}$, $r_{k}$, $\alpha_{kl}$ with $l\neq k$ and initial conditions $u^{(k)}_{0}=\mathop{\rm const}\nolimits$:

• •

  Two-species: $\varepsilon_{1}$, $\varepsilon_{2}$, $\alpha_{12}$, $\alpha_{21}$ and $u_{0}=[u_{0}^{(1)},u_{0}^{(2)}]$.

• •

  Three-species: $\varepsilon_{1}$, $\varepsilon_{2}$, $\varepsilon_{3}$, $\alpha_{12}$, $\alpha_{13}$, $\alpha_{21}$, $\alpha_{23}$, $\alpha_{31}$, $\alpha_{32}$ and $u_{0}=[u_{0}^{(1)},u_{0}^{(2)}u_{0}^{(3)}]$.

We perform 10,000 simulations for each case with random parameters

and $0.01<u_{0}^{(k)}<0.99$. Two cases of the diffusion coefficient is considered: regular diffusion ($\varepsilon_{k}=D_{k}$) and small diffusion($\varepsilon_{k}=D_{k}/10$).

In investigation we use influence of the following parameters to the system solution: $\varepsilon_{k}$ and $\alpha_{kl}/r_{k}$ for $l\neq k$ and $k=1,2$ for two-species competition and $k=1,2,3$ for three species competition. In Figure 11, we present correlation matrices as well as loading matrix for each factor analysis for all three cases 1D, 2D(a) and 2D(b) with regular and small diffusion. In Table 1, we present a summary for factor analysis of the two-species competition model. From the correlation matrix, we observed a strong correlation between final population density and diffusion rate and between final population density and growth and competition efficiency.

From the loading matrix of factor analysis, we observed that in a two-species system, when the diffusion rates of both species are regular, the dominant factor that causes variation is the diffusion rates of both species. The second most crucial factor is both species’ growth and competition efficiency. However, when the diffusion rates of both species are low, the dominant factor that causes variation becomes the growth and competition efficiency of both species instead of the diffusion rate. 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 effect of boundary is more severe when the diffusion rate is larger.

In Figure 12, we present three - species competition system correlation matrices as well as loading matrix for each factor analysis for all three cases 1D, 2D(a) and 2D(b) with regular and small diffusion. In Table 1, we present a summary for factor analysis of the three-species competition model. We observed that the three-species competing system is similar to the two-species competing system. There is a strong correlation between final population density and diffusion rate and between final population density and growth and competition efficiency. 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 and competition efficiency of the three species.

For two-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 1D, 2D(a), and 2D(b) spatial boundary conditions, and summarized the proportion of different survival status groups, as is in Table 3. We observed that 1D and 2D(b) spatial boundary conditions give similar proportion of survival groups, while 2D(a) gives more different proportion. When there is regular diffusion, for 1D and 2D(b) spatial boundary conditions, there is around 30 percent of simulations that result in either 00(no species survives), 01 (species two survives), or 10 (species one survives). In comparison, only around 10 percent of simulations result in both species surviving together. For 2D(a) spatial boundary conditions, 68.88 percent of simulations result in 00 (both species go extinct), while around 14.87 percent and 14.89 percent result in 01 (species two survives) or 10 (species one survives), and very rare, about 1.36 percent 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 (2D(a) scenario).

When there is small diffusion, all three spatial boundary conditions (1D, 2D(a), and 2D(b)) give a similar proportion of survival status. It is most likely that both species survive together (around 50 percent of simulations). However, it is almost impossible that both species will go extinct. Furthermore, it is equally likely that only one species survives (proportion of simulation both around 20-25 percent).

For 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 1D, 2D(a), and 2D(b) spatial boundary conditions, and summarized the proportion of different survival status groups, as is in Table 4. We observed that 1D and 2D(b) spatial boundary conditions give similar proportion of survival groups, while 2D(a) gives more different proportion. When there is regular diffusion, the probability of one species survives (001, 010, and 100) is similar to one another (around 10 percent for 1D/2D(b) scenario, and around 12 percent for 2D(a) scenario), while the probability of two species survive (011, 110 and 101) is similar to one another(around 12 percent for 1D/2D(b) scenario, and around 1 percent for 2D(a) 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 2D(a) scenario which approximate the lake, around 58 percent of simulation result in 000 (no species survive). When there is small diffusion, all three spatial boundary conditions (1D, 2D(a), and 2D(b)) give a similar proportion of survival status. It’s equally likely that two species survive or all species survive (the proportion for survival group 011/110/101/111 are all around 20 percent), and it’s equally likely that only one species survive (the proportion for survival group 001/010/100 are all around 6-8 percent), and it’s almost impossible that all three species go extinct.

A mathematical model of the spatial-temporal multi-species competition is considered. A discrete system 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 effect 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. Furthermore from general cases where we release holds of diffusion while keeping other parameters fixed, we observed that different combination of diffusion rates of species in the system leads to different final survival status of species. When the combination of both species’ diffusion rates is on the border of different 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 for 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 was 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.