ABSTRACT

The distributions of species in their habitats are constantly changing. This phenomenon is thought to be determined by species’ environmental tolerance and biotic interactions for limited resources and space. Consequently, predicting the future distribution of species is a major challenge in ecology. To address this problem, we use mathematical model to study the combined effects of biotic interactions (e.g. competition) and environmental factors on multiple species community assembly in a heterogenous environment. To gain insights into the dynamics of this ecological system, we perform both analytical and numerical analyses of the range margins of the species. We observe that the range margins of the species can be influenced by biotic interactions combined with environmental factors. Depending on the strength of biotic interactions, our model exhibits coexistence of species and priority effects; mediated by weak and intense biotic interactions respectively. We also show the existence of bifurcation points (i.e. the threshold values of competition coefficient) which lead to the presence—absence of different species. Thus, we suggest that adequate knowledge of biotic interactions and changes in the environments is important for effective maintenance of biodiversity and conservation management.

Keywords:

Abiotic Environments, Biotic Interactions, Coexistence, Priority Effects

1. Introduction

Ecologists have long been interested in answering the fundamental question of what determines the geographical distributions of species [1] and coexistence mechanisms among them. Which species will be present or excluded in a geographical region, where, and why? These are challenging questions in ecology, and it is the main thrust of this study. This is because population occurs in milieu of other species, which often is characterized with multiple and complex species interactions that affect each species performance, population dynamics and then geographic ranges [2] . Darwin [3] is one of the earliest ecologists who demonstrated that biotic interactions such as competition can ultimately influence species spatial distributions and impact range shifts of the species. The novelty of his research initiated several experimental studies [4] [5] [6] [7] which demonstrated biotic interactions, such as competition as one of the significant drivers of species distributions. Today, theoretical evidences exist which also recognized that competition interactions can determine species range margins and as well influence community assembly in geographic environments [8] [9] [10] .

Another crucial force that can affect species distributions and range margins, is environmental factors in the form of climate change. This is because species have distinct responses to environmental changes and biotic interactions [11] . Previous studies have shown that climate change can affect species range margins and biodiversity across a geographical location [12] [13] . Most often, climatic conditions usually have influence on species’ natural population in the form of shifts in their geographic ranges and distributions, changes in abundance, or changes in individual behavior [14] . Thus, climate influence on the distribution of species, often acts as a limiting factor on the extent and location of species’ range margins. Environmental factors have been demonstrated in empirical studies to have been responsible for the distributions of species in their habitats. For example, Perry et al. [15] illustrated that fish distributions have shifted in mean dept (or latitude) due to temperature increase. Similarly, Rowe et al. [11] demonstrated that extreme cold spells have strong impacts on the distributions and abundances of tropical species, particularly at the verge of their individual temperature tolerance.

In addition, the order in which species become established may alter community assembly through priority effects (i.e. the effects one species have on community developments due to early arrival to site). In this case, the outcomes of species interactions depend on the initial abundances of the species; mediated by intense biotic interactions, such that species coexistence is impossible [16] . In the absence of priority effects and when biotic interactions are relatively weak, species can coexist [8] [17] . Loureiro et al. [4] conducted a study using Daphnia species to illustrate that priority effects outcomes can be affected by abiotic factors (e.g. salinity) and then alter community dynamics. Experimental study exists which also, manipulates initial abundance to illustrate occurrence of priority effects [18] . Relative abundance of two species of invasive grasses in California (Bromus diandrus and Bromus madritensis) was empirically demonstrated to have altered depending upon arrival order [19] . Also, empirical study using a microbial community [20] has shown that disparity in the timing of species can lead to alternative community assembly.

However, theoretical studies that analyze the combined effects of biotic and abiotic factors on multispecies dynamics with priority effects are rare. Thus, there is more to know about multispecies dynamics, when competition interactions occur in communities of more than two or three species in a geographically changing environment. According to May and Leonard [21] , when competition interactions involve more than two species, it is likely for things to get complicated and many different community structures can be expected. Therefore, understanding how biotic interactions combined with environmental factors to determine multispecies community assembly is crucial. Thus, we aim to proffer possible explanations for where and when coexistence of species and priority effects are likely to occur using mathematical modelling. This information may be useful for management plans that lessen biodiversity fatalities.

Thus, we extend the deterministic model of Godsoe et al. [10] to model biotic interactions of multispecies community dynamics across heterogenous environments. The model is a spatially explicit extension of the Lotka-Volterra competition model, which is a system of ordinary differential equations (ODEs). To gain insight, we analyze the model with the assumption that species competition strengths are symmetrical. Thus, the remaining part of the study is organized as follows: the next section presents and describes the deterministic model which incorporates biotic interactions and environmental factors. This is followed by derivations of analytical results on the range margins of the species. Thereafter, we present the results of numerical simulations on the range margins of species. We also, generate summary plot and bifurcation analysis as a model’s parameter varies; in order to show the presence-absence of species across a geographical region. Then, we discuss the implications of our results.

2. Materials and Methods

2.1. The Models

We consider a multispecies deterministic model which is an extension of Lotka-Voltera competition model. The model is a system of ordinary differential equations for the densities N_i(t, x) of n species extended along one dimensional environmental gradient x, as in Equation (1) where 0 ≤ x ≤ 1.

$\frac{\text{d}{N}_i}{\text{d}t}=\frac{r_i{N}_i}{k_i\left(x\right)}\left(k_i\left(x\right)-{\displaystyle \sum _{j=1}^n{\alpha }_{ij}{N}_j}\right);\left(i=1,2,\cdots ,n\right)$ (1)

Here r_i is the intrinsic growth rate, k_i(x) is the carrying capacity, α_ij is the competitive strength of species j on species i, α_ii is the intraspecific competition coefficients (i.e. a measure of the strength of competition within the same species) and N_i is the densities of species i at time t. To make the model simple, we set intraspecific competition coefficients α_ii to 1. Thus, interspecific competition coefficients, α_ij represent the ratio of interspecific to intraspecific competition coefficients. In this study, we consider competition of four species (i.e. n = 4) and assume that competition strengths of the species are symmetrical (i.e. α_ji = α_ij = α). Then, Equation (1) becomes a system of four ordinary differential equations, one for each species with competition coefficient α.

The environmental suitability is modelled into the carrying capacity, k_i(x) of the species as a spatial dependence on the locations x (i.e. the carrying capacity, k_i(x) vary with the locations x). In this case, the environmental location is represented by x; which is used as a proxy to represent environmental factors like temperature, moisture, salinity and any other environmental factors that may affect the species. Therefore, the effects of biotic interactions on the species may depend on how the species respond to environmental factors. To show these effects, the carrying capacity is modelled such that it varies linearly with x.

$k_i\left(x\right)=m_{i}x+b_i$ (2)

In Equation (2), b_i is the species i carrying capacity at x = 0; m_i denotes the changes in the suitability of the environments with respect to abiotic component x (i.e. the gradient of the carrying capacity) and k_i(x) represents the carrying capacity of the species i at each location x. In this case, the maximum density that species i can attain is at x = 1 if m_i is positive.

2.2. Numerical Methods

To understand the dynamics of the system, $\left(d{N}_i\right)/dt$ in Equation (1) is set to zeros and we solved numerically for the steady states using MAPLE package. The stability analysis of the steady states is also carried out using MAPLE package. Thus, at a location x, the steady state who’s all the real parts of the eigenvalues are negative is stable. Based on the steady states and the numerical simulations results on the range margins, we used the techniques of invasion analysis to derive analytical results on the species’ range margins of the Equation (1). We also show from the invasion points of the species, the threshold values of the competition coefficient which can lead to presence-absence of different species across the locations x. The numerical simulation results on the range margins of the species are obtained by employing MATLAB ode15 solver for t = 1000 to solve Equation (1) until steady states are achieved. We also generate numerical simulation for the summary plot, using MATLAB ode15 solver in order to show different species present and their range margins across the geographical locations x as a model parameter, α changes. MAPLE package is also used to verify that the simulation results are stable. To further cross check the simulation results, numerical simulation package XPPAUT is used. Thus, the steady states of Equation (1) is computed with the aid of cvode solver for t = 1000. We then continued the steady states in AUTO; where the stable and unstable steady states as well as the bifurcation points are tracked as α changes. A parameter value, 10⁻¹/10⁻⁶ is used as the maximum/minimum allowable step size. Other parameter values that are used in the computation of the results are given in Table 1. Otherwise, the parameter values are written below the respective Figures.

3. Results

3.1. Analytical Results on Range Margins of Species

The analytical results of Equation (1) are presented in this section to show the range margins and threshold values of competition coefficient ( ${\alpha }_{T_i}$ ) of the species across locations x. The analytical results of the Equation (1) are derived using the method of invasion analysis to obtain the invasion points, x_i (i.e. the positions across the locations x where species i can invade when rare). The method uses the criterion that a species that can invade at a location x, must be rare at that point and its growth rate must be greater than zero (i.e. $\left(d{N}_i\right)/dt>0$ ) [22] . The potential of the species to invade depends on the carrying capacities and competition coefficients of the species. Godsoe et al. [10] used the idea of invasion analysis to find the range margins of two species of Equation (1) with n = 2, by obtaining the invasion points of the species; and suggested two ecological scenarios that can shift the range margins of the species. We extend the techniques of analyzing the two species invasion points of the Equation (1) to the case of four species that is considered in this study. The method of the analysis is by tracking the species that can invade at a location x when rare, and the species that are present at that location. The analysis requires the steady states (computed with the aid of MAPLE), competition coefficient (α) and carrying capacities of the species. The potential of a species to invade depends on its carrying capacity and competition coefficient. As stated in Table 1, species 1 - 4 carrying capacities are k₁(x) = m₁x, k₂(x) = m₂x, k₃(x) = b₃ and k₄(x) = b₄ respectively.

Based on Figure 1(a), the presence of species 3 and 4 corresponds to the invasion of species 1. However, for species 1 to invade, it requires that the right-hand side of $\left(d{N}_1\right)/dt$ in Equation (1) to be greater than zero. In this case, species1 is considered rare at its invasion point and so, its density (i.e. N₁) is taking to be zero. In a similar way, the density of species 2 (i.e. N₂) is also considered to be zero; since at the invasion point of species 1 only species 3 and 4 are present. Thus, species 1 can invade if $k_1\left(x\right)>\alpha N_3^{*}+\alpha N_4^{*}$ , where $N_3=N_3^{*}$ and $N_4=N_4^{*}$ . Therefore, the point x where:

$k_1\left(x\right)=\alpha N_3^{*}+\alpha N_4^{*},$ (3)

satisfies species 1 invasion point (x₁). But k₁(x) = m₁x and for stable steady state, $\left(0,0,N_3^{*},N_4^{*}\right)=\left(0,0,\left(\alpha k_4-k_3\right)/\left({\alpha }^2-1\right),\left(\alpha k_3-k_4\right)/\left({\alpha }^2-1\right)\right)$ , where k₃ = b₃ and k₄ = b₄, which on substitution into Equation (3), species 1 invasion point is given as:

$x_1=\frac{{\alpha }^2\left(b_3+b_4\right)-\alpha \left(b_3+b_4\right)}{m_1\left({\alpha }^2-1\right)}$ (4)

Equation (4) implies that there exists an asymptote (i.e. x₁→∞) in the range margin of species 1 whenever α = 1 (i.e. the range margin of species 1 tends to infinity at α = 1). Thus, based on the Equation (4), the threshold value of competition coefficient of species 1 in Equation (6) can be established through Equation (5).

$m_1{x}_1\left({\alpha }^2-1\right)=\left({\alpha }^2-\alpha \right)\left(b_3+b_4\right)$ (5)

Thus,

${\alpha }_{T_1}=\frac{m_{1}x}{\left(b_3+b_4\right)-m_{1}x}$ (6)

Moving along the environmental gradient x to the right of x₁ (see Figure 1(a)), we observed the presence of species 1, 3 and 4 with species 2 absent. Thus, species 2 is the one that can invade in the presence of the three species. The invasion point (x₂) of species 2 can be derived like species 1. In this case, $\left(d{N}_2\right)/dt$ in Equation (1) must be greater than zero for species 2 to invade. Also, at the invasion point of species 2, its density (i.e. N₂) is also, considered equivalently zero. Thus, species 2 can invade if $k_2\left(x\right)>\alpha N_1^{*}+\alpha N_3^{*}+\alpha N_4^{*}$ , where = $N_1=N_1^{*}$ , $N_3=N_3^{*}$ and $N_4=N_4^{*}$ . Then, species 2 invasion point satisfies that:

$k_2\left(x\right)=\alpha N_1^{*}+\alpha N_3^{*}+\alpha N_4^{*}$ , (7)

where k₂(x) = m₂x, k₃ = b₃, k₄ = b₄ and the steady state is given as:

$\begin{array}{l}\left(N_1^{*},0,N_3^{*},N_4^{*}\right)=(\left({\alpha }^2\left(k_3+k_4-k_1\right)-\alpha \left(k_3+k_4\right)-k_1\right)/\left(2{\alpha }^3-3{\alpha }^2+1\right),\\ 0,\left({\alpha }^2\left(k_1+k_4-k_3\right)-\alpha \left(k_1+k_4\right)-k_3\right)/\left(2{\alpha }^3-3{\alpha }^2+1\right),\\ {\alpha }^2\left(k_1+k_3-k_4\right)-\alpha \left(k_1+k_3-k_4\right)/\left(2{\alpha }^3-3{\alpha }^2+1\right)).\end{array}$

Thus, substituting into Equation (7), gives the invasion point of species 2 as in Equation (8).

$x_2=\frac{{\alpha }^3\left(b_3+b_4\right)+\left(1-2\alpha \right)\alpha \left(b_3+b_4\right)}{m_2\left(2{\alpha }^3-3{\alpha }^2+1\right)-\alpha m_1\left({\alpha }^2-2\alpha +1\right)}$ (8)

Like species 1, the threshold value of competition coefficient of species 2 becomes:

${\alpha }_{T_2}=\frac{m_{2}x}{\left(b_3+b_4\right)+m_{1}x-2m_{2}x}$ (9)

Based on Equation (8), the first scenarios by which species 2 range margin is likely to increase depends on strong interspecific competition from species 3 and 4 at the boundary of species 2 fundamental niche. For instance, larger values of α(b₃ + b₄) can shift species 2 invasion point (x₂) and then, increase the range of x for which species 2 can be absent. Therefore, a modest change in the model parameters will also cause a modest change in the range margin of species 2 [10] . Secondly, species 2 range margin can also increase when the denominator in Equation (8) gets smaller. In this case, the range margin of species 2 can shift if there is similarity between the changes in the combined effects of species 1, 2, 3 and 4 (i.e. $m_2\left(2{\alpha }^3-3{\alpha }^2+1\right)$ ) and the changes in the combined effects from species 1, 3, and 4 (i.e. $\alpha m_1\left({\alpha }^2-2\alpha +1\right)$ ); along the environmental gradients (αm₁). As the denominator tends to zero, there exists an asymptote (i.e. x₂ →∞). In this scenario, small changes in the model parameters will result to high changes in the invasion point (x₂) of species 2 [10] .

In a similar way, species 4 invasion point can be computed in the presence of species 1, 2 and 3. Thus, invasion point of species 4 is given as:

$x_4=\frac{b_4\left(2{\alpha }^3-3{\alpha }^2+1\right)-\alpha b_3\left({\alpha }^2-2\alpha +1\right)}{\left({\alpha }^2-2\alpha +1\right)\alpha \left(m_1+m_2\right)}$ (10)

The threshold value of competition coefficient of species 4 is given as in Equation (11).

${\alpha }_{T_4}=\frac{b_4}{\left(m_1+m_2\right)x-2b_4+b_3}$ (11)

At the location x = 0.5 for instance, the threshold value of competition coefficient of species 4is equal to the threshold value of competition coefficient of species 2. This implies simultaneous extinction of species 4 and 2 at the same competition strength and location.

Also, the invasion point of species 3 is computed to give:

$x_3=\frac{b_3\left({\alpha }^2-1\right)}{\left(\alpha -1\right)\alpha \left(m_1+m_2\right)}$ (12)

Like species 1, the range margin of species 3 tends to infinity if α = 1 in Equation (12). This implies that the range margins of species 1 and 3 both tend to infinity at α = 1; since x₃ = x₁ at the value of α = 1. Based on the Equation (12), the threshold value of competition coefficient of species 3 is given as in Equation (13).

${\alpha }_{T_3}=\frac{b_3}{\left(m_1+m_2\right)x-b_3}$ (13)

The threshold value of competition coefficient of species 3 is also, equal to the threshold value of competition coefficient of species 1 at the location x = 0.5 and α = 1.25. See Equation (6) and parameter values given in Table 1.

3.2. Numerical Results on Range Margins of Species

Numerical results are presented in this section to illustrate the influence of biotic interactions and environmental gradients on the range margins of species across locations x. Both the numerical and analytical results agreed with each other (compare Figure 1(a) with analytical results). As stated in Table 1, we study competitions of two pairs of ecologically similar species. In each pair, the two species are ecologically similar based on their carrying capacities. For instance, species 3 and 4 are both homogenously distributed across the locations x, which means that both species have the same environmental tolerance. In the same way, species 1 and 2 are also ecologically similar species, since the two species have the same environmental requirements (i.e. their carrying capacities increase from left to right as x increases from 0 to 1). The carrying capacities of the species are also stated below the Figures.

To show the impacts of competitive strengths on multispecies community assembly, the numerical results are obtained separately for different competitive strengths of the species. For example, when competitive strengths of the species are relatively weak (see Figure 1), this situation leads to multispecies coexistence at the center. Contrarily, when competition strengths of the species are relatively strong (see Figure 2), it leads to exclusion of species such that only one species occupies a location x. In this case, the initial abundances of the species determine the outcomes of the competitions. Summary plot is also presented in this section, which illustrate the presence-absence of the species across the locations as α varies with respect to locations x (see Figure 3(a)). We also presented the bifurcations analysis results, showing the bifurcation points and the densities of the species present at a location x as α varies (see Figure 3(b)). The bifurcation analysis and the summary plot result also agreed with each other on the

presence-absence of the species. A detection threshold value of 0.5% of the maximum observed density of the species is employed for the numerical results presented in this section. This means that a species is considered absent if its density is below this expected value [23] [24] .

3.2.1. Multispecies Range Margins Due to Weak Biotic Interactions (i.e. α <1)

In Figure 1, the solid lines indicate steady states of the species, dotted lines represent the carrying capacities and circles on the horizontal axis denote the invasion points of the species. Figure 1(a) illustrate competition outcomes of the species when interspecific competition is weaker (i.e. α < 1) than intraspecific competition coefficient. In this case, the impact of competitions of one species on other species is relatively weak. Consequently, coexistence of the species at one location is possible. We observed multispecies species coexistence at the central location x. However, due to competition, only species 3 and 4 coexisted on the left side of the locations x and then exclude species 1 and 2 from the locations x. Thus, species 1 and 2 range margins are shifted from the origin to a point x₁ = 0.3375 and x₂ = 0.4655 respectively. Similarly, species 3 and 4 are shifted by species 1 and 2 on the right side of the locations x;such that only species 1 and 2 coexisted in the domain. Hence, species 4 and 3 range margins are shifted to a point x₄ = 0.5370 and x₃ = 0.7407 respectively. The multispecies coexistence at the central region, illustrates the most favorable locations in the environments. The range margins shown in Figure 1(a) correspond to the invasion points of our analytical results (compare the circles in Figure 1(a) with analytical results in Equations (4, 8, 10, 12)). We observed that biotic interactions and environmental factors combined to determine the range margins of the species.

To illustrate the impact of the environmental gradients, we computed Figure 1(b) by changing the carrying capacity of species 2 from k₂(x) = 0.8x to k₂(x) = 2x – 0.8. Other parameter values remain as in Figure 1(a). Consequently, the community assembly illustrated in Figure 1(a) is altered due to changes in environmental gradients. Thus, we observed that only two and three species coexistence are possible across the location x. The order of invasion points (x_i) of the species changed to x₁, x₄, x₂ and x₃ in contrast to community assembly observed in Figure 1(a).

3.2.2. Multispecies Range Margins Due Intense Biotic Interactions (i.e. α > 1)

In this section, we investigate the effects of intense biotic interactions (i.e. α > 1) on competition outcomes among the species. In this case, we observed that coexistence of species in one location is impossible due to aggressive interactions. As a result, this situation leads to occurrence of priority effects; and the dynamical model’s behavior is such that species’ initial abundances determine the outcomes of the competitions. Ecologically, priority effects can be referred to as alternative stable states; a situation where the presence-absence of the species depends on the order of arrival to site [25] . For example, Figure 2 illustrates the range margins of the species resulting from α = 1.4 with initial abundances favoring species 2 and 4 in the first row (i.e. Figure 2(a), Figure 2(b)). To illustrate the impacts of initial abundance on competition outcomes, the initial abundance of species 2 is reduced in the second row (Figure 2(c), Figure 2(d)) to N₂(x) = 0.75k₂(x); and in the third row (Figure 2(e), Figure 2(f)), the initial abundance of species 4 is reduced to N₄(x) = 0.75k₄(x). This implies that in Figure 2, the initial abundances of the species in Figure 2(a) and Figure 2(b) are given as: N₁(x) = 0.1k₁(x), N₂(x) = 0.9k₂(x), N₃(x) = 0.1k₃(x), N₄(x) = 0.9k₄(x); in Figure 2(c) and Figure 2(d), the initial abundances of the species are given as: N₁(x) = 0.1k₁(x), N₂(x) = 0.75k₂(x), N₃(x) = 0.1k₃(x), N₄(x) = 0.9k₄(x) and in Figure 2(e) and Figure 2(fe), the initial abundances of the species are: N₁(x) = 0.1k₁(x); N₂(x) = 0.9k₂(x); N₃(x) = 0.1k₃(x); N₄(x) = 0.75k₄(x). In each case, the spatial domain is partitioned into four regions, depending on the initial abundances of the species (see Figure 2). Each region is occupied by a single-species, existing to the maximum density of the carrying capacity

Throughout, species 2 and 4 are more abundant and so, have higher potential to exclude species 1 and 3 from some locations and occupied a larger domain. However, species 3 and 4 dominated the left side of the domain (see Figures 2(a)-(f)) and excluded species 1 and 2 in the region. Species 4 with higher initial abundance compared to species 3 dominated the larger part of the left region and shifted species 3 to a narrower central region. Similarly, species 1 and 2 dominated the right side of the domain (see Figures 2(a)-(f)) and exclude species 3 and 4. Species 2 with initial abundance advantage, also dominated the larger part of the right domain and then shifted species 1 towards a smaller right center. We observed that the decrease in species 2 initial abundance (see second row) tipped a balance more for species 1 and increase its region of dominance. Similarly, the reduction of initial abundance of species 4 (see third row) tipped a balance more for only species 3. These observations show that in a multispecies community dynamic, initial abundance can determine the range margins of the species; but ecologically similar species may likely have more impacts on one another than the reverse.

3.3. Summary and Bifurcation Plots of Species Presence-Absence as Competition Coefficient Varies

Figure 3(a) is the summary plot which illustrate Equation (1) predictions on the presence-absence of species as competition coefficient, α changes with respect to environmental locations x. The plot is generated using the carrying capacities of species as described in Table 1. Four sets of initial abundances are used; in each set, the initial abundances favor one of the four species at a time. Colors are used to represent the range margins of the species present at a location x. Different combinations of species present are illustrated by changes in color across the geographical locations x. The location where one-color changes to another, signifies critical value of the competition coefficient; where a stable combination of species exchange its stability for another one. In Figure 3(a), the regions labelled bistable, tristable and tetrastable indicate two, three and four stable combinations of single-species presence respectively at the same location x.Each single-species exist to the maximum density of the carrying capacity. In this case, the single-species it converges to depend on the initial abundances of the species.

When competition coefficient, α < 1, multispecies coexistence (red colored region) is observed at the center of the location x. Due to changes in the environments and biotic interactions, the number of species coexistence decreases as one move away from the center. These observations are also, illustrated in Figure 1(a). However, increasing competition coefficient α beyond species 2 and 4 threshold value of competition coefficient of 0.6667 at the location x = 0.5 (see Figure 3(a) and refer to Equations (9 and 11)), this situation leads to species 2 and 4 simultaneously excluded from the four species coexistence. Thus, only species 1 and 3 are observed to coexist (blue colored region), which also disintegrate at the location x = 0.5 and α = 1 to give rise to bistable single-species (lime colored region) of species 1 and 3. Beyond α > 1.25 (see Figure 3(a) and refer to Equations (6 and 13)), the dynamical model’s behavior illustrates tetrastable single-species (gray colored region); each at the maximum density of the carrying capacity. The single-species it converges to depend on the initial abundances of the species.

To further understand the different species presence-absence in the summary plot, we employed numerical continuation to compute Figure 3(b). The stable and unstable steady states of the species in the Equation (1) are tracked as α varies. In Figure 3(b), it illustrates the stable steady states densities of species 1 at the location x = 0.5. Any other species other than species 1 can be used for the plot. The Figure is computed with the same carrying capacities of the species as in Figure 3(a); and as specified in Table 1. The numerical continuation results illustrate several stable and unstable steady states of species presence-absence and existence of threshold values of competition coefficient (α). The threshold values of competition coefficient (α) in the numerical continuation results, correspond to the critical values, α in the summary plot. The threshold values of competition coefficient α correspond to transcritical bifurcation (i.e. α_t1) and saddle node bifurcation (α_T2) points in Figure 3(b). These bifurcation points lead to different branches of stable (red line) and unstable (black line) steady states of the species. In Figure 3(b), there exist four species coexistence stable steady state branch when α < α_t1 (compare Figure 3(b) with Figure 3(a)). The bifurcationpoint, α_t1 is also, consistent with species 2 and 4 threshold values of competition coefficients at the location x = 0.5 (see Figure 3(b) and refer to Equations (9 and 11)). This is followed by another stable steady-state branch of species 1 and 3 coexistence (i.e. α_t1 < α < α_T2); with simultaneous exclusion of species 2 and 4. The vertical red line at α = 1 indicate an asymptote, where the range margins of species 1 and 3 tend to infinity (refer to Equations 4 and 12). Beyond α > α_T2, are another stable steady state branch of species 1 and 3, each existing as a single-species to the maximum density of the carrying capacity. This is separated by unstable steady state (black line) of two species coexistence. Based on our analytical results (refer to Equations (6 and 13)) and the summary plot (Figure 3(a)), the bistable single-species of species 1 and 3 undergoes bifurcation at α = 1.25 to exchange its stability with tetrastable single-species for α > 1.25. The results in this study are also cross checked in MAPLE and found to be consistent.

4. Discussion

We investigate the dynamical outcomes of multispecies competition in an environmentally changing habitat. One of the significant results we observed is that when biotic interactions are relatively weak, species can coexist, with multispecies coexistence near the center. This form of community assembly has earlier been observed in an empirical study of small mammal species along elevational gradients [26] [27] . This may be attributed to locations where the environments are moderately suitability with low competition intensity on all the species. Consequently, the exclusion of the species at the lower and upper environments of the locations, as illustrated in Figure 1(a) is expected, due to unfavorable environments on some species; coupled with competition interactions from the environmentally favored species. Connell [28] reported in his empirical research using Balanus and Chthamalus species, that competition interactions and environmental factors can combine to determine the presence-absence of species. Therefore, the conservation and diversity of species may depend on relatively weak biotic interactions and moderate environmental components. In this way, both species can favorably compete for space and resources without anyone being eliminated from the community.

It is also observed that when biotic interactions are severe (α_j > 1), this situation can lead to occurrence of priority effects. In this case, the initial abundances of the species determine the presence-absence or the range margins of the species. Ecologically, these qualities could be implemented in biocontrol management; for preserving some species of interest or to reduce the excesses of a species whose activities are undesirable in the habitat. As illustrated in Figure 2, species 1 for instance, can be sustained in the habitat by introducing more of it against species 2. Similarly, if species 4 activities are not desirable in the habitat, the activities can be limited by introducing more of the species 3. According to Godsoe et al. [10] , if a focal species is ecologically similar to its competitor, a small change in the biology of either species can radically change the range margins of the species. It is also, previously observed that initial abundance can be used in biocontrol management, to regulate the presence-absence of species [29] [30] .

Also, environmental factors combined with biotic interactions are illustrated in our model results (compare Figure 1(a) with Figure 1(b)) to affect the range margins of the species. This observation agreed with [31] who reported that, species distributions can be influenced by combined effects of biotic interactions and environmental factors. Therefore, if the outcomes of species’ competitions depend on combined influence of biotic and environmental factors, this calls for caution in the model’s parameterization; as different magnitude of carry capacity, k_i(x) may engender different dynamical behaviors. However, the choice of k_i(x) which affects only species 3 and 4 carrying capacities, will qualitatively correspond to the dynamical behaviors illustrated in this study. For instance, changing the slope of the carrying capacities of species 3 and 4 from m₃ = 0 to m₃ < 0 (respectively from m₄ = 0 to m₄ < 0), will not change the invasion points of the species as much from the one illustrated in this study. This is because, we still have two pairs of ecologically similar species interacting with one another. In this case, a small variation in the steepness of species 3 and 4 carrying capacities, will lead to a modest change in the invasion point of the species [10] . Consequently, the dynamical outcomes are expected to be consistent with our results, as the model’s parameter α changes across the locations x.

Also, our numerical continuation results which illustrate both stable and unstable steady states and bifurcation points of the models, proffer detail explanation for the different species presence-absence observed in the numerical simulation results. Our model exhibits existence of threshold values of competition coefficient, α. The threshold values correspond to critical values (or color change) in the summary plot, where one combination of species presence exchanged its stability for another combination of species. The bifurcation points therefore, give rise to different dynamical behaviors of the model; such as coexistence, exclusion of species and priority effects.

5. Conclusion

In conclusion, we used deterministic model, which is a system of ordinary differential equations, to investigate the combined effects of competition interactions and environmental factors on multispecies community dynamics. The model is analyzed for the range margins of the species using analytical and numerical methods. Both analytical and numerical simulation results are found to agree with each other. The findings of this study have shown that biotic interactions and environmental factors can combine to determine the range margins of species. The results revealed different dynamical outcomes, which depending on the species’ competition strengths can be coexistence of species or priority effects. Thus, based on the findings, we suggest that adequate knowledge of biotic interactions and changes in the environments is essential for successful maintenance of biodiversity and conservation management. Also, ecological factors such as dispersal may change the outcomes of the competition dynamics presented in this study. We therefore, recommend dispersal inclusion in the deterministic model in this study; as an interesting extension which can lead to robust predictions of the range margins of species across a geographical region.

Acknowledgements

Thanks to the School of Mathematical Sciences and the Universitti Sains Malaysia for the support and encouragements.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Omaiye, O.J. and Mohd, M.H. (2019) The Roles of Biotic Interactions and Environmental Factors on Multispecies Dynamics. Open Journal of Ecology, 9, 426-442. https://doi.org/10.4236/oje.2019.910028

References