Long-Term Conservation Effects of Protected Areas in Stochastic Population Dynamics

Terrestrial and marine protected areas are essential tools in mitigating anthropogenic impacts and promoting population persistence and resource sustainability. Adequately implemented protected areas (PAs) aim to promote conservation by increasing population size and reducing its variability. To resolve how these effects depend on PA features, I develop and analyze new models of stochastic processes that encompass the fluctuations generated by demographic or environmental stochasticity in PAs management. The stochastic model is built upon individual processes. In the model, density-independent mortality, migration between PAs and non-PAs, organism preference for PAs, and size characterize the features of the PA. The effect of PAs size is also examined. The long-term conservation effects are quantified using the coefficient of variation (CV) of population size in PAs, where a lower CV indicates higher robustness in stochastic variations. The results from this study demonstrate that sufficiently reduced density-independent mortality in PAs and high site preference for PAs and immigration rate into PAs are likely to decrease the CV. However, different types of stochasticity induce rather different consequences: under demographic stochasticity, the CV is always reduced because PAs increase the population size therein, but an increased population size by PAs does not always decrease the CV under environmental stochasticity. The deterministic dynamics of the model are investigated, facilitating effective management decisions.


INTRODUCTION
Terrestrial and marine protected areas are being expanded worldwide in response to increasing concern about species loss (Watson et al., 2014). These protected areas (PAs) have become essential tools for mitigating anthropogenic impacts, promoting population persistence and resource sustainability, and enhancing ecological resilience (IUCN, 2008;Venter et al., 2014;Watson et al., 2014;Sala and Giakoumi, 2018). The establishment of PAs is not in itself a goal, but PAs are assumed to promote long-term conservation (Margules and Pressey, 2000). Existing strategic PA site selection methods face difficulties in envisioning their long-term conservation effects, because site selection often involves a snapshot of optimality, rather than a long-term consideration of the optimum (Possingham et al., 2000).
In the literature, several long-term benefits of PAs arise in deterministic models under the assumption that populations approach an equilibrium state after the PAs are implemented, such as improvement of yields (Takashina, 2016) and mitigating bycatch (Hastings et al., 2017) in fisheries management. This concept is ubiquitous and characterizes an expected long-term effect of conservation practice and management (e.g., Clark, 1990;Holden et al., 2018). In nature, however, population size fluctuates according to demographic and/or environmental stochasticity, with the former attributed to the probabilistic nature of (intrinsic) demographic events and the latter attributed to the noise induced by external factors (Nisbet and Gurney, 1982;Lande et al., 2003). Consequently, ecological status (and, hence, PA effects) fluctuates over time (Figure 1A), perhaps in different ways under different types of stochasticity. In such situations, effects of PAs are no longer persistent, and a population may undergo a period during which a given conservation effect is weak, potentially enhancing the risk of extinction. Therefore, variations in conservation effects, along with their average effects, are critically important to understand effects of PAs. However, prevailing equilibrium discussions cannot deal with such variations in conservation effects. Nonequilibrium discussions, however, can inform long-term PA effects (Figure 1B), and are the focus of this study.
In practice, long-term field observations are very costly; hence, there are limited time-series data available (Geldmann et al., 2013). Instead, modeling has power to effectively explore such long-term effects. Most modeling frameworks integrating stochastic events based on PA effects are in the context of marine protected areas (MPAs) (e.g., Mangel, 2000a;Aiken and Navarrete, 2011;Hopf et al., 2019), which are often associated with fisheries management and focus on optimal harvesting strategies (Costello and Polasky, 2008), tradeoffs between conservation effects and fisheries profits (Mangel, 2000b;Grafton et al., 2005;De Leo and Micheli, 2015), population persistence (Aiken and Navarrete, 2011;White et al., 2020), and resilience (West et al., 2009;Barnett and Baskett, 2015; FIGURE 1 | (A) Population size in protected areas (PAs) show variations over time; and (B) such variations characterizes the long-term effect of PAs. A large variation indicates that the effect is variable over time, while a small variation suggests a robust long-term conservation effect. Aalto et al., 2019). Few studies have investigated the stochastic influences on long-term MPA effects in a fisheries context, or how MPAs affect variability in population size or catch (Mangel, 2000a;Grafton et al., 2005;Fryxell et al., 2006). Barnett and Baskett (Barnett and Baskett, 2015) demonstrated that, by assuming stochastic recruitment in predatory fish species, the coefficient of variation in the catch could be reduced in the case of fisheries management using MPAs. In a study of metapopulation dynamics in the Great Barrier Reef Marine Park, Hopf et al. (2019) demonstrated that marine reserves could promote the stability of populations and fishery yields, regardless of fishing intensity. However, Hopf et al. (2019) also showed that this conclusion depends on the location of reserves: more variable biomass can occur if disturbed reefs are protected, compared to protecting undisturbed sites. This suggests a need to understand the underlying mechanisms in determining of conservation effects. The models used in these studies target marine species and often have complex structures to describe the life histories of species. Hence, our understanding of stochasticity and PA effects is very limited, despite their broad applicability to terrestrial and marine ecosystems.
This paper develops general theoretical insights regarding the long-term effects of stochasticity on PAs (and MPAs), without restricting the discussion to fisheries management or marine environments. I focus on (i) whether PAs suppress stochastic fluctuations (i.e., providing a long-term conservation effect) and under what conditions, if any, this is achieved; and (ii) whether demographic and environmental stochasticity affect the longterm conservation effects of PAs in different ways. The developed master equation allows existing analytical methods to be used, and I develop analytical insights into the stochastic population model. This would be a difficult task using the existing complex models. This is not merely for mathematical understanding, but provides a more explicit underlying mechanism of variability in the effects and parameter dependence of PAs.
In this paper, I address these questions using a spatially explicit stochastic population model (e.g., McKane and Newman, 2004;Hakoyama and Iwasa, 2005;Gardiner, 2009). I begin by looking at individual processes and formulating a master equation, and then obtain the corresponding stochastic differential equations (SDEs). While the derived stochastic population models with PAs are new, it turns out that they are stochastic analogous of existing deterministic models used to analyze equilibrium properties in MPA management. I measure the time variations of PA effects via the coefficient of variation (CV), where a smaller value indicates a more robust long-term conservation effect against stochasticity, and vice versa. I use the CV to quantify time-variable PA effects, which are scaled by the mean population size (note that, in general, the variance increases with the mean), but I present the CV along with the mean value.
My approach enables the long-term conservation effect of PAs under stochasticity to be discussed. Multiple parameters determine the quality of PAs, and I can examine various biological and management scenarios. Additionally, by defining "inappropriate" PAs as sites with a higher mortality than non-PAs (e.g., caused by illegal use of protected species within PAs Razafimanahaka et al., 2012;Harasti et al., 2019), we discuss how inappropriately enforced/managed PAs affect conclusions. This approach provides the opportunity to discuss effective PA implementations under stochastic population dynamics.

METHODS
In the model, the focal region has two categories: PAs and non-PAs (Figure 2). I examine the CV in PAs, non-PA, and the whole region for various sizes of PAs, different levels of quality (measured by the degree of decline of mortality rate therein), and several parameters such as site preference, and the degree of stochasticity. The immigration and emigration of individuals to FIGURE 2 | The model scheme. The concerned region is subdivided by subdivisions and these have two categories: PA (fraction R; green region) and non-PA (fraction 1-R; white region). The emigration from PAs and the immigration to PAs exchange individuals between the areas at constant rates of m 1 and m 2 , respectively. Each area is characterized by an intrinsic growth rate r i and a carrying capacity K i (i = 1, 2), and site preference affects the likelihood of individual migration (see the main text). and from PAs connect these areas. The sizes and site preferences affect the likelihood of individual migration. Inhomogeneous mortality rates and sizes in the two regions induce different growth rates r i and carrying capacities K i , which characterize the quality of PAs-higher-quality PAs offer lower mortality in the region, leading to a larger conservation effect. First, I discuss a simple situation in which the population size remains fixed, and introduce some key aspects of the model and analytical results. Then, I discuss a more general situation in which demographic stochasticity occurs in birth, death, and migration events. I also make the explicit connection to existing deterministic models, and thereby introduce SDEs with environmental stochasticity.
When PAs have higher growth rates than non-PAs, it is possible to swap the definition of the two areas. That is, one can regard areas having lower growth rates as PAs, and discuss the effect of "inappropriate" PAs. In the following, I will provide results for both PAs and non-PAs, and the discussion of "good" PAs encompasses "inappropriate" PAs.

Population Dynamics With Fixed Population Sizes
I begin with a simple situation in which the population dynamics of a focal species are driven by the immigration and emigration of individuals to and from PAs, and no birth or death events occur. Each area has a site preference, which affects the realized migration rate. The realized migration rate determines the degree of mixture between PAs and non-PAs. Detailed technical discussions can be found in Appendix A1. The number of individuals remains fixed at N, and I write the population sizes in PAs and non-PAs as n and N − n, respectively. Let p(n, t) be the probability of n individuals located in PAs at time t. Then, the population dynamics can be described by a simple gain-loss process (Appendix A1): dp(n, t) dt = g n−1 p(n − 1, t) + l n+1 p(n + 1, t) − (g n + l n )p(n, t), (1) where g n and l n are the gain and loss rates corresponding to one individual gain and loss in PAs, respectively. Let N s and n s be the total number of subdivisions of the concerned region and the number of subdivisions categorized as PAs, respectively, and let R = n s /N s be the fraction of PAs in the region of interest (Figure 2). Only the immigration and emigration of individuals drive population changes in PAs, and the gain and loss terms become where m 1 and m 2 are the emigration and immigration rates, respectively. The parameter α controls the preference of PAs, accounting for preferred/non-preferred/neutral sites by the species (Figure A1). Neutral preference (α = 1) means that the destination of an individual on the move is determined at random and weighted by the sizes of the PAs and non-PAs. When PAs are preferred (α < 1) or non-PAs are preferred (α > 1) by the species, the probability of choosing PAs is higher or lower than the neutral choice, respectively.

Population Dynamics Under Demographic Stochasticity
When birth and death occur in a population, the total population number N is no longer constant, but is the sum of the population sizes in PAs n 1 and non-PAs n 2 : N = n 1 + n 2 . When the population size changes dynamically as a result of births, deaths, and migrations, Equation (1) becomes (see Appendix A2 for full details) where r is a state and W(n | m) is the transition rate from state m to state n. These transition rates are where b is the birth rate, d i is the density-independent mortality rate in site i (i = 1, 2), and d is the density-dependent mortality rate. The factors R and 1 − R in Equations (4c, 4d) denote the fractions of PAs and non-PAs, respectively. For example, this accounts for a larger density-dependent mortality in a smaller region. The intrinsic growth rate r i and the carrying capacity K i in site i (i = 1, 2) are defined as r i = b − d i and K i = r i /d, respectively (see Appendix A2 for details). The connections to existing deterministic models can be observed by deriving a deterministic representation of Equation (5). Multiplying both sides by n 1 and summing over all probabilities, I recover a deterministic two-patch dynamics equation (Appendix A2): where E[x] represents the average of x. This model is discussed, for example, in Takashina et al. (2017) and Takashina (2020) with some arrangements.

Population Dynamics Under Environmental Stochasticity
I next introduce a stochastic model incorporating environmental stochasticity. This can be obtained from the representation of SDEs (Turelli, 1977;Lande et al., 2003). Under the deterministic dynamics of Equation (5) and with continuous variables X i representing the population size in PAs (i = 1) and non-PAs (i = 2), the SDEs can be written as Gardiner (2009) dX 1 = M 1 (X 1 , X 2 )dt + σ e X 1 dW(t), where M i corresponds to the deterministic part of the population dynamics in region i (Equation A19) and σ e is the intensity of environmental stochasticity.

Long-Term Effects of PAs With a Fixed Population Size
When population changes in PAs are solely due to migrations between PAs and non-PAs, I can obtain analytical insights that are not possible for more general situations. Later, I will show that analytical insights from this simple situation can be applied to these more general situations. In the model, I can derive explicit forms of the mean and the CV in PAs as follows (Appendix A1): The parameter dependence on the CV in PAs is analyzed by differentiating (Equation 7b) with respect to the parameter of interest. For example, the relationship ∂CV/∂m 2 < 0 indicates that an increase in the immigration rate to PAs m 2 decreases the CV, and ∂CV/∂α > 0 indicates that an increase of the site preference for non-PAs α increases the CV, respectively. In addition, I conclude that an increase in the PA fraction always decreases the CV in PAs because ∂CV/∂R < 0. Figure 3 represents the relationships among the expressions in Equation (7), and confirms the above discussions. I can further simplify (Equation 7b) by setting E[n] = n in Equation (7a), solving for an arbitrary parameter, and plugging the result into Equation (7b). This cancels R, α, m 1 , and m 2 in Equation (7b), and the CV in PAs can be described only by the population size in PAs, n, and the total population size, N (Equation A14 in Appendix A1). This indicates that a larger population size in PAs results in a smaller CV. Hence, I conclude that any parameters that improve the population size in PAs reduce its CV.

Long-Term Effects of PAs Under Demographic Stochasticity
When birth and death events occur in addition to migrations, the total population size changes dynamically. Although this situation is not amenable to mathematical analysis because of the nonlinearity in the demographic rate terms, numerical simulations imply that I can still discuss this situation following a similar line to the analytical discussions above. For instance, Figures 4A,B show that the site preference α, immigration rate to PAs m 2 , and PA fraction have a qualitatively similar dependence on the CV in PAs to that of the fixed-population scenario (Figure 3).
Regarding the dependence of demographic parameters, I can use the relationships obtained in the preceding analysis. That is, changes in parameter values that increase the population size in PAs will reduce its CV (Figures 4C,D). This can be heuristically seen using the result in the preceding analysis (Appendix A1), whereby a larger population size in PAs, n 1 , leads to a smaller CV in PAs, while the magnitude of the population size in non-PAs, n 2 , scales with the CV: Generally speaking, expanding the PAs increases their population size and decreases that in non-PAs. If a conservation effect of PAs is sufficient, it also improves the total population size N. This leads to a smaller CV in PAs as the PA fraction increases (Figure 4). Figure A2 shows the results for non-PAs and the whole region. The opposite trend from that in PAs can be observed:  Frontiers in Ecology and Evolution | www.frontiersin.org increasing the PA fraction increases the CV in non-PAs, which is associated with a population decline in the region. Additionally, the total population size determines the CV in the whole region.

Long-Term Effects of PAs Under Environmental Stochasticity
Environmental stochasticity affects the whole population, and its influence does not vanish even with large population sizes, unlike demographic stochasticity. Under environmental stochasticity, the relationships discussed above are no longer useful, and the interpretation of results becomes more intricate. For example, increasing the PA fraction and population size in PAs does not guarantee a reduction in the CV, but can instead increase the CV (Figure 5B; immigration rates m 2 =0.1 and 10). However, these actions may reduce the CV when m 2 = 1.0. Similarly, the site preference α exhibits a complex response in terms of the effect on the CV in PAs (Figure 5A), and for PAs with a higher preference, increasing the PA fraction may increase the CV (Figure 5A; α = 0.8 and PA fraction R < 0.2).
The establishment of PAs is likely to decrease the corresponding CV if a conservation effect of PAs is increased by reducing the density-independent mortality rate in PAs, d 1 (Figure 5D). Increasing the birth rate b (a focal species has a high fecundity rate) shows a similar effect, and tends to decrease the CV regardless of the PA fraction ( Figure 5C). However, a larger birth rate decreases the conservation effects of PAs (e.g., the relative difference between the net growth rate of PAs and non-PAs (b − d 1 )/(b − d 2 ) decreases when b increases), and the effect of reducing the CV in PAs becomes smaller. These numerical observations can be verified analytically when the immigration and emigration rates are sufficiently large (m 1 , m 2 ≫ 1). In this condition, the CV in PAs is described as (Appendix A3) wherer can be interpreted as an average population growth rate in the whole region weighted by the site preference (α). From Equation (9), I conclude that an increase in the PA fraction decreases the CV in PAs if the growth rate in PAs is greater than that in non-PAs (b − d 1 > b − d 2 ), and vice versa. I also state that a small PA size has a large effect on reducing the CV when the site preference for PAs is higher than that for non-PAs (α < 1), and vice versa (Appendix A3).
Results for non-PAs and the whole region are provided in Figure A3, but the CV trends in these cases are not easy to distinguish from those in PAs. This highlights the difference from the situations under demographic stochasticity (Figure 4, A2; bottom), where the opposite trends in population sizes and CVs occur between the two areas. This may be because environmental stochasticity affects all individuals regardless of the population size, and the stochastic influence tends to be consistent across the whole region. Additionally, the CV in the whole area is more likely to decrease with increasing total population size ( Figure A3). However, this is not the case when the immigration rate is large (Figure A3F; m 2 = 10).
When PAs provide a higher conservation effect (d 1 = 0.1), the trend for an increasing CV in PAs would be mitigated ( Figure  A4B; m = 0.1, 10), or even reversed ( Figure A4A; α = 0.8). I also examined a situation with a higher environmental stochasticity of σ 2 e = 0.3, and obtained qualitatively similar results ( Figure A5).

General Findings
A multitude of indices can characterize the effect of PAs (Leverington et al., 2010). Here, I have analyzed the long-term conservation effects of PAs under stochasticity, which is not simply an equilibrium discussion. I measured the long-term PA effects via the CV of the population size, where PAs with a small CV offer less variation, and hence have a larger longterm conservation effect. Loosely speaking, the CV in PAs is suppressed by increasing their area if the PAs provide a sufficient conservation effect. Hence, the implementation of effective PAs promotes long-term conservation effects. This general finding agrees with previous studies focused on the marine environment (Barnett and Baskett, 2015;Mellin et al., 2016;Hopf et al., 2019).

Effects of Demographic/Environmental Stochasticity
Demographic stochasticity is suppressed by a large population size (Lande et al., 2003). Therefore, if PAs increase the population size by, for example, expanding in size or increasing the immigration probability (Figures 3, 4A,B), or reducing the mortality rate (Figure 4D), the population fluctuations are suppressed; this is a mechanism for achieving longterm conservation effects. Moreover, if PAs increase the total population size in the focal region, fluctuations in the population are suppressed by the same mechanism (Figures A2E-H). The developed master equation approach allows us to derive a concrete analytical insight for a fixed population size, e.g., under the assumption that the target species has a low growth rate compared to the timescale of a conservation program. In fact, many conservation practices act over much shorter timescales than ecological timescales (Willis et al., 2005;Froyd and Willis, 2008).
In contrast, environmental stochasticity affects all individuals, and population size plays a minor role in determining the long-term conservation effects (Figures 5A,B,D). However, if the conservation effect of PAs is high (e.g., mortality is sufficiently reduced in PAs), expanding the PAs tends to reduce the CV. Equivalently, the CV in PAs is increased by introducing ineffective PAs (e.g., mortality is not sufficiently reduced or increased in PAs). In addition, if PAs increase the total population size in the whole region, the CV tends to decrease in that region (Figure A3), except for situations with a high immigration rate, where the opposite trend is observed (e.g., Figure A3F; m 2 = 10).

Inappropriate PAs May Amplify Fluctuations
In the proposed model, the density-independent mortality rate characterizes the quality of PAs, and I implicitly assumed that PAs reduced the mortality rate in those areas. In fact, this definition is arbitrary, and one can argue that PAs have a higher mortality than non-PAs, a situation that could be described as "inappropriate" PAs. The results of this case are already discussed as "non-PAs" in this paper (e.g., Figures A2A-D). Alternatively, an increase in mortality inside PAs would reverse patterns identified here (i.e., proceeding right to left along x-axes in Figures 3-5): increasing PA fraction would increase overall mortality, thereby increasing the CV.
In practice, poorly implemented PAs can arise due to inappropriate planning or management process, such as lack of sufficient commutations between stakeholders (Agardy et al., 2011). Existing unfairness affects compliance of management and can increase conflict between user groups (Agardy et al., 2011). Illegal use of protected species is one of the potential risks of PAs management that can increase the mortality of a target species (Razafimanahaka et al., 2012). For instance, Harasti et al. (2019) reported that illegal fishing activities potentially reduced the abundance of a fish species by 55% from 2011 to 2017 in the Seal Rocks no-take area in Australia. Similarly, Hopf et al. (2019) demonstrated certain PAs designs can further destabilize a system. Hence, effective enforcement of PAs is necessary to promote its benefit.
Strategies to Achieve Robust PA Management I have identified different mechanisms for enhancing the population size in PAs, such as those to improve the demographic rate (i.e., birth b and death d 1 ) and to promote the probability of remaining in PAs (i.e., site preference α, emigration m 1 and immigration m 2 rates). In practice, if PAs adequately regulate anthropogenic activities and are monitored (e.g., strict nature reserve IUCN, 2008), the demographic rate may be improved in PAs, providing a long-term conservation effect under demographic/environmental stochasticity. However, ineffective PA management is often associated with a failure to reduce human activities (Craigie et al., 2010;Schulze et al., 2018), which may amplify population fluctuations. The immigration and emigration rates are not purely biological parameters, but can be controlled by the configuration of the PAs (Takashina, 2020). Protecting the favored sites of a target species is desirable as a means of increasing the conservation effects of PAs (Hunt et al., 2020). However, when MPAs are used as a sustainable fisheries management tool, a moderate spillover effect is necessary to promote fishing yields (McClanahan and Mangi, 2000;Stobart et al., 2009). Under environmental stochasticity, this can lead to a reduced long-term effect of PAs, and careful assessment is necessary.
MPAs can improve the resilience of populations under existing multiple stable states (Takashina and Mougi, 2014;Barnett and Baskett, 2015;Aalto et al., 2019). The findings in this paper further complement this insight. Namely, while a catastrophic shift of population is incurred by population perturbations (Scheffer et al., 2001), adequately established PAs tend to suppress population fluctuations (i.e., the population is more robust against perturbations). Therefore, populations within PAs are more likely to remain in the basin of attraction of a current stable state. A more explicit discussion addressing the relationships between the CV of the population dynamics, the strength of perturbations, and the basin of attraction will further improve our understanding.
Likewise, there are multiple directions for further extending this analysis. For example, here, I have assumed that environmental stochasticity affects PAs and non-PAs equally. While this is reasonable when the concerned region is small, and the environment in each type of area is not different, two areas may be subject to other environmental fluctuations (Hakoyama and Iwasa, 2005). In fact, the size of an MPA sometimes becomes significant (Dulvy, 2013), and the model developed needs to incorporate heterogeneous environmental stochasticities. Age and metapopulation structures have been used in the context of fisheries management (Hopf et al., 2019), and these investigations are also relevant to the context of this paper. While the study discussed the expected population size along with the CV provided by a fraction of PAs established, an expected timescale to observe such a population size is an important management consideration (Kaplan et al., 2019;Barceló et al., 2021). Kaplan et al. (2019) demonstrated recovery is fast at small population sizes where variations of population fluctuations are also small. Revealing the role of the CV in population recovery will further complement the current knowledge. However, it should be noted that applying complex models has a large cost in terms of reduced generality and analytical intractability. With this in mind, the approach of the study will guide the development of a general framework to discuss the long-term conservation effects of PAs.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author. And I did not detect any particular expressions.

AUTHOR CONTRIBUTIONS
NT conceived the idea, conducted the analyses, and wrote the manuscript.

FUNDING
The University of Tokyo provided funding for this project. Also, this work was partially supported by JSPS KAKENHI Grant no. 21K17913.