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):

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₁ and m₂ 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.

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₁ and non-PAs n₂: N = n₁ + n₂. 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₁ 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.

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)

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.

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₂ < 0 indicates that an increase in the immigration rate to PAs m₂ 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₁, and m₂ 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.

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₂, 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₁, leads to a smaller CV in PAs, while the magnitude of the population size in non-PAs, n₂, 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: 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.

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₂=0.1 and 10). However, these actions may reduce the CV when m₂ = 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₁ (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₁)/(b − d₂) 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₁, m₂ ≫ 1). In this condition, the CV in PAs is described as (Appendix A3)

where r^ 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₁ > b − d₂), 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₂ = 10).

When PAs provide a higher conservation effect (d₁ = 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 σe2=0.3, and obtained qualitatively similar results (Figure A5).

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 long-term 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).

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 long-term 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₂ = 10).

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.

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₁) and to promote the probability of remaining in PAs (i.e., site preference α, emigration m₁ and immigration m₂ 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.