Critical patch sizes in spatial population dynamics with two habitat fragments

Abstract

This work starts from the Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) equation to model diffusive populations in fragmented regions characterized by patches with different environmental conditions. Using this model, we found an expression for the general case of critical patch size, i.e. sizes below which the population goes extinct in a system of two patches (life-beneficial regions) surrounded by sinks (life-detrimental regions). This expression alone is an interesting result because it allows for the study of phenomena and interpretations presented in this work, as well as other features not included in this study. From the analysis of this expression, we extracted the analytical prediction that in a system formed by a patch and two semi-infinite sinks, if we replace one semi-infinite sink with a patch plus another semi-infinite sink, more lethal than the original, the added patch will only be beneficial to the original system if it is large enough to compensate for the added sink. This interpretation arises from a case study done with a specific set of parameters, but will appear for other values that meet the presented condition.

Keywords:
Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) equation; Minimum size of fragmented regions; Relation between neighbor patches

1. Introduction

Habitat fragmentation is a problem for the survival of different species living in several environments, from mammals scattered around the globe [¹] to bacteria in the laboratory [², ³]. Concerns about the relation between the fragment size and population density have also been discussed in the literature [⁴], concluding that larger and more intact areas are more beneficial for species preservation [⁵]. However, small fragments should not be overlooked, as they have a fundamental contribution to biodiversity and species conservation [⁶], sometimes preventing some of them from extinction [⁷].

Mathematical models developed decades ago [⁸, ⁹] and their successors [¹⁰–¹³] have sought to provide insights into finding the minimum size of a fragment that can sustain a stable population within it, with experimental confirmation already documented [³]. Experimental results serve as important guides for extending or refining models to accurately describe reality.

From the literature [¹², ¹⁴, ¹⁵], in a system of two communicating patches, both need a smaller size to provide life inside them, compared with a single patch. The decrease in the minimum size of each fragment comes from the communication between them through a region where the diffusion of the species in question is possible, but the region does not have resources to support the population indefinitely. However, these studies only considered identical fragments, raising questions about whether this behavior would remain the same for asymmetric fragmentation.

In this sense, this work is a natural continuation of previous studies [¹¹, ¹⁴, ¹⁵]. The main focus of this research is to determine whether the presence of a small fragment for the species in question has any impact on another nearby fragment. This other fragment is connected by an unfavorable region for the species’ life, but it allows the diffusion of individuals within the population between the fragments. The study explores whether the insertion of a very small fragment near an existing fragment, that already meets the minimum size to sustain a stable population, will increase, decrease, or not affect the minimum size of the original fragment.

To answer this question, it is necessary to predict the behavior of the patch size in a system of two totally asymmetric patches. Mathematically, this implies to find an analytic expression for the minimum size, as found in the literature for particular cases [¹¹, ¹², ¹⁴, ¹⁵]. Here, in this paper, the objective is to present the general case.

The paper is disposed as it follows: section ² presents the problem and discusses known techniques and results that will be useful to obtain new results of this work and their comparison with the literature. In section ³ a brief mathematical discussion of the main steps to obtain an analytical prediction of the minimum patch sizes for the general asymmetric case of homogeneous regions is presented. Section ⁴ will address results as well as their discussion. The concluding remarks close the paper.

2. Model

The Fisher-Kolmogorov-Petrovskii-Piskunov (FKPP) equation has been used to model population dynamics from microscopic entities as genes [⁸], bacteria [², ³] and other cells [¹⁶], to macroscopic ones, such as spruce budworm [¹³], such that its application in the study of diffusive populations is relevant.

In population dynamics, one feature is the existence of a critical size L_c, below which a stable population cannot exist within the fragment [⁸, ⁹, ¹³]. This critical size can be understood as the minimum size for the preservation of a species, or as the maximum size that guarantees the extinction of this species, i.e. if the patch size L is smaller than L_c, the species will be extinct from this patch.

The convenient form of the one-dimensional FKPP equation found in the literature [², ⁹, ¹³] as a mathematical model for describing a diffusive population is:

(1) $$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+a(x)u-b{u}^2,$$

where the variables related to population are: u = u(x, t) the population density, a(x) the growth rate, b ≥ 0 the saturation rate or intraspecific competition, and D the dispersion coefficient. The x and t variables represent space and time, respectively.

The a(x) profile is commonly used to represent space fragmentation, because if a(x) ＜ 0, in Eq. (1), the population density u(x, t) goes to zero for large times (t → ∞), which does not happen if a(x) ＞ 0. However, if a(x) is a composition of regions where a(x) positive and regions where a(x) negative, it is possible a particular solution that satisfies Eq. (1) in the considered domain. A common interpretation is to assume as patches (or habitat fragments), the regions where a(x) ＞ 0, surrounded by unsuitable regions, where a(x) ＜ 0, this last one will be labeled here as sinks, see Fig. (1).

Figure 1
Representation of growth rate profiles used to describe the general case of (a) a two patches system and (b) a single patch.

In the proposed model, the simplest possible fragmentation is to merge homogeneous regions (where a(x) is constant) favorable and unfavorable to life. For this type of fragmentation, the general form of a two fragment system is represented by Fig. (1a), which is composed of two patches and three sinks, as described in Table 1.

Tabele 1
Definition of regions and their properties.

Short label	Region label	Space region
S 1	sink 1	x ＜ −L1
P 1	patch 1	−L1 ≤ x ≤ 0
Ss	sink s	0 ＜ x ＜ s
P 2	patch 2	s ≤ x ≤ L2 + s
S 2	sink 2	x ＞ L2 + s

In Fig. (1b), the general case of a single patch is presented. Its critical size was recently described by Pamplona da Silva [¹¹] in explicit form:

(2) $$L_{sph}=\sqrt{\frac{D}{a_i}}\left\{\arctan \left(\sqrt{\frac{h_i}{a_i}}\right)+\arctan \left(\sqrt{\frac{p}{a_i}}\right)\right\}.$$

From Eq. (2) the critical size for a semi-isolated fragment is easily obtained at the limit h_i → ∞, i.e.,

(3) $$L_{spi}=\sqrt{\frac{D}{a_i}}\left\{\frac{\pi }{2}+\arctan \left(\sqrt{\frac{p}{a_i}}\right)\right\}.$$

Figure (1b) can be obtained from Fig. (1a) taking the limit s → ∞. The variable h_i is used to represent h₁ and h₂ and a_i to represent a₁ and a₂, alluding now to fragments 1 and 2, respectively.

3. Mathematical Analysis

To describe a time-varying population, it is necessary a model with temporal dependence, in which the population cannot grow infinitely. In this sense, Eq. (1) is minimally qualified to model a diffusive population dynamics at any point (x, t) of spacetime. However, as already discussed by Pamplona da Silva et al. [¹¹, ¹⁴, ¹⁵] and Kenkre and Kumar [¹²], following the ideas of Ludwig et al. [¹³]. The focus of this study is on stable solutions over time, that is, after the transient period. Thus, the temporal evolution term (∂u/∂t) in Eq. (1) can be suppressed. The intended results also focus on a limit condition for the fragment size. This condition corresponds to a size below which the fragment cannot sustain a stable population within it. The quadratic term can be neglected (when compared to the linear term) since u tends to zero at the extinction-survival transition, with u² tending to zero faster. Therefore, for the purposes of this study, it is possible to work only with the linear stationary equation, namely:

(4) $$D\frac{{\partial }^2\text{Φ}}{\partial x^2}+a(x)\text{Φ}=0,$$

where Φ = Φ(x) = u(x, t) for t → ∞.

To determine the critical sizes L₁ and L₂ of fragments 1 and 2, respectively, Eq. (4) is discussed, region by region, (see Table 1), enforcing the continuity of the function and its derivative at the boundaries. Thus, the functions, which are solutions of Eq. (4) in five regions, are initially identified:

(5) $$\mathrm{\Phi }{s}_1(x)=A{e}^{k_{1}x}\text{+}A1e^{-k_{1}x},$$

(6) $$\mathrm{\Phi }{P}_1(x)=C\text{sin (}{\alpha }_{\text{2}}x)+D\cos ({\alpha }_{2}x),$$

(7) $$\mathrm{\Phi }{S}_e(x)=G{e}^{k_{3}x}+H{e}^{-k3x},$$

(8) $$\mathrm{\Phi }{P}_2(x)=E\text{sin}\text{(}{\alpha }_{\text{4}}x\text{)}\text{+}F\text{cos}({\alpha }_{\text{4}}x),$$

(9) $$\mathrm{\Phi }{S}_2(x)=B_1{e}^{k_{5}x}+B{e}^{-k_{5}x},$$

where $k_1^2=h_1/D$, ${\alpha }_2^2=a_1/D$, $k_3^2=p/D$, ${\alpha }_4^2=a_2/D$ e $k_V^2=h_2/D$.

From the boundary conditions ${\text{Φ}}_{S_1}(-\infty )=0$ and ${\text{Φ}}_{S_2}(\infty )=0$, the constants A₁ = 0 and B₁ = 0 were obtained respectively. With the functions (solutions of Eq. (4) in their respective regions) given by Eq. (5), Eq. 6), Eq. (7), Eq. (8) and Eq. (9), it is enforced continuity at borders: x = −L₁, x = 0, x = s and x = L₂ + s, respectively obtaining the expressions:

(10) $$A{e}^{-k_{\text{1}}{L}_{\text{1}}}\text{=}-C\text{sin (}{\alpha }_{\text{2}}{L}_1)+D\cos ({\alpha }_2{L}_1),$$

(11) $$D=G+H,$$

(12) $$G{e}^{k\text{3}s}\text{+}H{e}^{-k\text{3}s}\text{=}E\text{sin (}{\alpha }_{\text{4}}s)+F\cos \text{(}{\alpha }_{\text{4}}s),$$

(13) $$E\text{sin}{\alpha }_4\text{(}{L}_{\text{2}}\text{+}s\text{)}\text{+}F\cos {\alpha }_4(L_2+s)=B{e}^{-k_5(L_2+s)},$$

as well as the continuity of the first derivative at the same points provides:

(14) $$k_{\text{1}}{\text{Ae}}^{-k_{\text{1}}{L}_{\text{1}}}\text{=}{\alpha }_{2}C\text{cos (}{\alpha }_{\text{2}}{L}_1)+{\alpha }_{2}D\sin ({\alpha }_2{L}_1),$$

(15) $${\alpha }_{2}C=k_3(G-H),$$

(16) $$k_{3}G{e}^{k_{3}s}-k_{3}H{e}^{-k_{3}s}={\alpha }_{4}E\cos ({\alpha }_{4}s)-{\alpha }_{4}F\sin ({\alpha }_{4}s),$$

(17) $${\alpha }_{4}E\cos {\alpha }_4(L_2+s)-{\alpha }_{4}F\sin {\alpha }_4(L_2+s)=-k_{5}Be-k_5(L_2+s).$$

To solve the algebraic system composed by Eq. (10) to Eq. (17), in order to determine the constants, several paths can be followed to eliminate redundant equations and to obtain a linearly independent system. If it is desired to get exactly and directly the same format shown in this paper, just follow these steps: from Eq. (10) and Eq. (14), obtain the relation C = R₁ D (see Eq. (19), below) and replace it in linear combinations of Eq. (11) and Eq. (15) resulting in one relation between G and D and other between H and D. Similarly, with Eq. (12), Eq. (13), Eq. (16) and Eq. (17) find a relation between D and F and other between H and F. Finally eliminate the constants G and H by obtaining a system in the variables D and F, which in matrix form can be expressed by:

(18) $$\left(\begin{array}{cc}\left(m_1{R}_1+m_2\right)e^{-k_{3}s}& -\left({\alpha }_2{R}_1+k_3\right)\\ \left(m_3{R}_2+m_4\right)e^{k_{3}s}& {\alpha }_2{R}_1-k_3\end{array}\right)\left(\begin{array}{c}F\\ D\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$$

where

(19) $$\begin{array}{lll}{R}_1\hfill & =\hfill & \frac{k_1\cos ({\alpha }_2{L}_1)-{\alpha }_2\sin ({\alpha }_2{L}_1)}{{\alpha }_2\cos ({\alpha }_2{L}_1)+k_1\sin ({\alpha }_2{L}_1)},\hfill \\ R_2\hfill & =\hfill & \frac{k_5\cos ({\alpha }_4(L_2+s))-{\alpha }_4\sin ({\alpha }_4(L_2+s)}{{\alpha }_4\cos ({\alpha }_4(L_2+s))+k_5\sin ({\alpha }_4(L_2+s))},\hfill \\ m_1\hfill & =\hfill & -[{\alpha }_4\cos ({\alpha }_{4}s)+k_3\sin ({\alpha }_{4}s)],\hfill \\ m_2\hfill & =\hfill & k_3\cos ({\alpha }_{4}s)+{\alpha }_4\sin ({\alpha }_{4}s),\hfill \\ m_3\hfill & =\hfill & {\alpha }_4\cos ({\alpha }_{4}s)-k_3\sin ({\alpha }_{4}s),\hfill \\ m_4\hfill & =\hfill & k_3\cos ({\alpha }_{4}s)+{\alpha }_4\sin ({\alpha }_{4}s).\hfill \end{array}$$

The nontrivial solution for the system of Eq. (18) generates the secular equation Eq. (20), which will be presented in Section ⁴.

A non-trivial solution means that the total population is finite. This solution provides the critical sizes of the patches. The critical size is the smallest size at which a fragment can sustain life within it. In this work, the critical sizes are represented by L_i, where the subindices i identify the different patches. In a system of two fragments, when the population of one is zero, the population of the other must also be zero. Otherwise, both populations could become non-null due to diffusion.

4. Results and Discussion

This work has a mathematical result, which is an interesting outcome by itself, since it generalizes, in the critical condition, the relation between the problem parameters (L₁, L₂, a₁, a₂, h₁, h₂, D, p, and s). This result arises from the secular equation associated with the system of Eq. (18),

(20) $$(m_1{R}_2+m_2)({\alpha }_2{R}_1-k_3)e^{-k_{3}s}+(m_3{R}_2+m_4)({\alpha }_2{R}_1+k_3)e^{k_{3}s}=0.$$

Equation (20) provides many possibilities of analysis, including a particular one, which has been a subject of interest in the literature [⁹, ¹¹, ¹³] and is the main phenomenological focus of this paper, namely, the relation between fragment sizes (L₁ and L₂), in critical condition. Specifically, how the presence of a patch affects its neighbor critical size, as discussed below.

In particular, Fig. (2) illustrates a phenomenon that arises from the approximation of two patches connected to each other by a sink with life difficulty (p). To observe this phenomenon, it is not necessary for the condition p to have a value equal to the difficulties of the external sides of the system (h₁ and h₂), as shown in Fig. (1a).

Figure 2
Critical patch size for parameters a₁ = 4.8, a₂ = 2.5, h₁ = 1.7, h₂ = 3.9, D = 1.1, p = 0.9 and s = 0.4.

Figure (2) shows, on the oblique continuous line, the relation between L₁ and L₂ for the parameter set, a₁ = 4.8, a₂ = 2.5, h₁ = 1.7, h₂ = 3.9, D = 1.1, p = 0.9 and s = 0.4. In this one, it is observed that, in the presence of a very small patch 2, the fragment 1 has a critical size (L₁) larger than if it was alone (L_1sph). Similarly, fragment 2 has its critical size increased by the proximity of a very small patch 1 (see regions L₁ ＞ L_1sph and L₂ ＞ L_2sph). At first glance, this result is counterintuitive, as the insertion of a second fragment, even if very small, should not be able to hinder the original fragment.

In the literature [¹², ¹⁴], the presence of a fragment always benefited its neighbor, even if only to a small extent. In other words, as we brought two patches closer together, their critical size would become smaller or equal to that of each patch when isolated. This assistance could be so minimal that it might be negligible. It was not expected to hinder its neighbor, as discussed earlier. However, in the cited literature, the fragments had symmetries or system configurations that always led to this positive interaction between them.

This counterintuition can be clarified by noting that in the example presented in Fig. (2), we had a system 1 (initially alone) that originally could sustain stable life within it. This system had a patch of size L₁ with parameters h₁ = 1.7, a₁ = 4.8, and p = 0.9, where p = 0.9 represented the condition of one of the sinks (let’s assume it was on the right side), see Fig. (1). When we replaced the sink on the right side of system 1 with a system 2, which originally could not sustain stable life within it, the dynamics changed. Although this new system had a “patch 2" with size L₂ (and parameters p = 0.9 (on the left), a₂ = 2.5, and h₂ = 3.9), this system 2 began to act as a more lethal sink than the original sink with condition p = 0.9 on the right of system 1.

5. Concluding Remarks

The main result of this work arises when, in a system initially alone, formed by a patch and two semi-infinite sinks, we replace one semi-infinite sink with a patch plus another semi-infinite sink. In this new dynamic, the critical size of the fragment in the original system increases in the presence of the second fragment. This result can be heuristically explained by observing that the new semi-infinite sink is more lethal than the original semi-infinite sink, so that adding a new fragment will be sufficient to compensate for this new more lethal sink only if the new patch is large enough.

The presented result was achievable solely through the analytic expression, Eq. (20), for the general case of the critical sizes for two patches with homogeneous conditions inside and outside them, i.e. spatial heterogeneities occur only abruptly and at the fragments borders. From this expression, only one parameter was evaluated, and there are still possibilities of this equation to contribute for future works related to this subject.

As known in the literature [¹², ¹⁴], the interaction between two identical patches connected by a sink, as long as it is not impenetrable (p ＜ ∞), is always beneficial to both, because they can be smaller than when they are alone. However, in cases where a very small patch separates one sink from another sink that is more lethal than the first, what the system experiences is this more lethal sink, not the small patch that separates them.

Acknowledgments

The authors (DJPS and TAB) thank the “Coordenação de Aperfeiçoamento de Pessoal de Nível Superior” (CAPES) and the (DJPS and ECO) “Departamento de Matemática Aplicada” (DMA-UNICAMP). The authors deeply thank the reviewer selected by RBEF for the careful and qualified review, which undoubtedly enhanced the quality of this article significantly.

References

• ^[1]
  K.R. Crooks, C.L. Burdett, D.M. Theobald, S.R.B. King, M. Di Marco, C. Rondinini, and Boitani L., Proceedings of the National Academy of Sciences 114, 7635 (2017).
• ^[2]
  V.M. Kenkre and M.N. Kuperman, Physical Review E 67, 051921 (2003).
• ^[3]
  N. Perry, Journal of The Royal Society Interface 2, 379 (2005).
• ^[4]
  M.A. Bowers and S. F. Matter, Journal of Mammalogy 78, 999 (1997).
• ^[5]
  E.F.Connor, A.C. Courtney and J.M. Yoder, Ecology 81, 734 (2000).
• ^[6]
  D. Lindenmayer, Proceedings of the National Academy of Sciences 116, 717 (2018).
• ^[7]
  B.A. Wintle, H. Kujala, A. Whitehead, A. Cameron, S. Veloz, A. Kukkala, A. Moilanen, A. Gordon, P.E. Lentini, N.C.R. Cadenhead et al., Proceedings of the National Academy of Sciences 116, 909 (2018).
• ^[8]
  R.A. Fisher, Annals of Eugenics 7, 355 (1937).
• ^[9]
  J.G. Skellam, Biometrika 38, 196 (1951).
• ^[10]
  R.A. Kraenkel and D.J. Pamplona da Silva, Physica A: Statistical Mechanics and its Applications 389, 60 (2010).
• ^[11]
  D.J. Pamplona da Silva, Ecological Modelling 384, 168 (2018).
• ^[12]
  V.M. Kenkre and N. Kumar, Proceedings of the National Academy of Sciences 105, 18752 (2008).
• ^[13]
  D. Ludwig, D.G. Aronson and H.F. Weinberger, Journal of Mathematical Biology 8, 217 (1979).
• ^[14]
  D.J. Pamplona da Silva and Kraenkel R.A., Physica A: Statistical Mechanics and its Applications 391, 142 (2012).
• ^[15]
  D.J. Pamplona da Silva, Villar R.P. and Ramos L.C., Applied Mathematics and Computation 315, 494 (2017).
• ^[16]
  K. Takamizawa, S. Niu and T. Matsuda, Journal of Biomaterials Science, Polymer Edition 8, 323 (1997).

Publication Dates

• Publication in this collection
  19 July 2024
• Date of issue
  2024

History

• Received
  09 Apr 2024
• Reviewed
  28 May 2024
• Accepted
  04 June 2024