Effects of stochasticity and social norms on complex dynamics of fisheries
Sukanta Sarkar, Arzoo Narang, Sudipta Kumar Sinha, Partha Sharathi Dutta
aa r X i v : . [ q - b i o . P E ] S e p Effects of stochasticity and social norms on complex dynamics of fisheries
Sukanta Sarkar , Arzoo Narang , Sudipta Kumar Sinha , ∗ and Partha Sharathi Dutta , † Department of Mathematics, Indian Institute of Technology Ropar, Punjab, India Department of Chemistry, Indian Institute of Technology Ropar, Punjab, India (Dated: September 30, 2020; Received :to be included by reviewer)Recreational fishing is a highly socio-ecological process. Although recreational fisheries are self-regulating and resilient, changing anthropogenic pressure drives these fisheries to overharvest andcollapse. Here, we evaluate the effect of demographic and environmental stochasticity for a social-ecological two-species fish model. In the presence of noise, we find that an increase in harvesting ratedrives a critical transition from high-yield/low-price fisheries to low-yield/high-price fisheries. Tocalculate stochastic trajectories for demographic noise, we derive the master equation correspondingto the model and perform Monte-Carlo simulation. Moreover, the analysis of probabilistic potentialand mean first-passage time reveals the resilience of alternative steady states. We also describe theefficacy of a few generic indicators in forecasting sudden transitions. Furthermore, we show thatincorporating social norms on the model allows moderate fish density to maintain despite higherharvesting rates. Overall, our study highlights the occurrence of critical transitions in a stochasticsocial-ecological model and suggests ways to mitigate them.
I. INTRODUCTION
Stochasticity, or only noise, is so prevalent that thestate of a natural system never remains constant. Noiseis often fascinating in its own right and can induce newand unexpected novel phenomena that can not be under-stood alone from the underlying deterministic skeleton[1–3]. Deterministic and stochastic processes interact inexciting ways and share a common origin [4]. Intrinsi-cally noisy and disordered processes generate strikinglyregular patterns similar to those systems with orderedprocesses [2]. Further, in the presence of alternativebasins of attraction, a stochastic event might shift a sys-tem from one attractor to another alternative attractorwith contrasting properties [5, 6]. Ecologists are majorlyconcerned about the temporal fluctuations observed innatural populations, which are the result of intrinsic dis-turbances (demographic stochasticity) and variations inthe global, extrinsic environmental changes (environmen-tal stochasticity)[7–11]. Disturbances can modify ecosys-tems dynamics for various reasons, and predicting ecosys-tems’ response towards these disturbances is a centralobjective in ecology and conservation. Therefore, un-derstanding the influence of natural and anthropogenicdisturbances on ecological variability needs considerableattention. However, failure to elucidate these featuresof stochastic population dynamics is accountable for theextinction of many vulnerable species, including recre-ational fisheries.Increasing anthropogenic pressure on ecosystems haslead to severe biodiversity loss [12, 13]. Further, har-vesting in the absence of adequate fishery managementpolicies has caused overexploitation of many species. Al-though recreational fisheries are, in general, thought to ∗ [email protected] † Corresponding author; [email protected] be self-regulating and highly resilient, a rise in overall de-mand due to human population growth is increasing thevulnerability of these fisheries to sudden collapse [14].According to a recent study [15], vast stocks of com-mercial fisheries in the United States and Europe havewitnessed colossal exploitation. Abundance has beenseverely reduced by harvesting that – reduced mortal-ity may not be sufficient to allow recovery of a popula-tion. Overfishing has driven various species on the vergeof extinction, including Atlantic cod, blue whale, thresh-ers, hammerheads [16, 17]. Also, excessive harvestingdecreases the interaction of overfished populations withother species in the community [18, 19], and eventuallydrives a sudden decline in their abundance. Some pop-ulations might retrieve [20] but mostly exhibit no signsof recovery, which is indeed a matter of concern. In thepresence of natural stochastic fluctuations, sudden tran-sitions from one equilibrium to another due to crossing atipping point (i.e., a threshold or a breakpoint) can makea recovery by species challenging. Therefore, analyzingtipping points at which such sudden shifts or critical tran-sitions occur in marine systems is an important area ofresearch [21–25].There is now growing awareness about the develop-ment of sustainable fishery management policies so thatsudden, irreversible decline of fish stocks can be avoidedwith changing environmental and economic pressures.Since it is well recognized that early warning signals(EWSs) of tipping points in social-ecological systems canadvert sudden transitions, they can be used for fisherymanagement [20, 26]. A lot of efforts has already beenmade to develop EWSs to forewarn critical transitions indiverse stochastic systems [6, 26–29]. Close to a tippingpoint, the recovery rate from perturbation becomes in-creasingly slow, termed critical slowing down, resultingin a loss of resilience in a system. In the presence of smallstochastic fluctuations, such critical slowing down is indi-cated by increasing recovery time, variance, and autocor-relation [6, 28]. The effectiveness of EWSs is supportedby both theoretical and empirical evidence [26, 29, 30].However, after much research on critical transitions andits indicators, the robustness of EWSs is still in question[31, 32]. A thorough understanding of ‘when’ and ‘how’EWSs work reliably is a major research question.This article investigates how a socio-ecological prey-predator fish population model responds to harvestingpressure in the presence of stochastic perturbations. In-creasing social/human influence on ecological systems isknown to exhibit a wide variety of dynamics and criticaltransitions in their structure and function [33]. How-ever, critical transitions and its early warning indicatorsin coupled socio-ecological systems remain less studied.As the harvesting efficacy pushes a system close to a crit-ical point, humans respond to declining fish-yield by in-creasing conservationism and making the system stable.This, in turn, allows moderate fish-yield to be maintaineddespite a higher fishing rate. Mathematically, this can beincorporated in the system by defining social norms andthe imitation dynamics of evolutionary game theory [34].Here, we study the influence of demographic noise intothe social-ecological model of fishing [35] by using themaster equation approach and then extending this modelfor environmental noise. First, we obtain stochastic tra-jectories for both demographic and environmental noises.Then we determine the resilience of alternative steadystates by probabilistic potentials and mean first-passagetime. As increasing harvesting rate in the presence ofstochastic perturbations drives critical transitions in themodel, we examine the efficacy of a few generic earlywarning indicators using simulated time series. When wesimultaneously vary the predation rate and prey harvest-ing rate, it reveals cyclic dynamics in the system. Fur-ther, we also suggest a realizable harvesting strategy byincorporating appropriate social norms into the systemfor the long-term sustainability of recreational fisheries.Our study shows that proposed social norms can effec-tively maintain moderate fish-yield despite a high rate ofharvesting.
II. MODELS AND METHODSA. Deterministic model
To understand the effects of harvesting on prey andpredator fishes, we study a socio-ecological two-speciesmodel where both the prey and the predator fishes areexposed to harvesting by fishers. The two-species deter-ministic model is given as [35]: dxdt = I + r x x (1 − xK x ) − pyx h + x − H x x, (1a) dydt = r y y (1 − yK y ) + cpyx h + x − H y y, (1b)where x and y represent prey and predator fish densities,respectively. I is the stocking rate, r x and r y are the respective growth rates of prey and predator fishes. K x and K y are the maximum potential biomass (or carryingcapacity) of x and y , respectively. p is predation rateand h is value of x where predation is half maximum. c is the consumption rate of the predator fish. H x and H y are harvesting rates of prey and predator fishes, re-spectively. Understanding dynamics of socio-ecologicalmodels is important as humans are responsible for de-cline in biodiversity, in turn decline in biodiversity badlyaffects humans. B. Stochastic model
While examining temporal changes in population dy-namics under natural fluctuations, stochastic modelingcomes into the picture. Noise can be introduced in de-terministic models in terms of intrinsic (demographic)and/or extrinsic (environmental) factors, and are ac-countable for causing population fluctuations in the ma-jority of species [36, 37]. Demographic stochasticityarises due to an individual’s mortality, reproduction,and primarily affects small populations. On the otherhand, environmental stochasticity is the variability im-posed by the environment on a population and mani-fested through random fluctuations in population growthrate. Stochasticity in models can bring new, unexpected,and novel phenomena while interacting with determinis-tic processes.
1. Presence of demographic noise
To study the influence of demographic noise, we con-sider different fundamental processes (e.g., birth anddeath processes) associated with the model (1), presentedin Table I. Using all these birth and death processes (gainand loss probability in Table I), we develop correspond-ing master equation for the grand probability function P ( x, y, t ), and it takes the following form: ∂P ( x, y, t ) ∂t = V IP ( x − , y, t ) − V IP ( x, y, t )+ r x ( x − − x − K x V ) P ( x − , y, t ) − r x x (1 − xK x V ) P ( x, y, t )+ p ( x + 1) y ( x + 1) + h V P ( x +1 , y, t ) − px yx + h V P ( x, y, t ) + ( x + 1) H x P ( x + 1 , y, t ) − xH x P ( x, y, t ) + r y ( y − − y − K y V ) P ( x, y − , t ) − r y y (1 − yK y V ) P ( x, y, t ) + cpx ( y − x + h V P ( x, y − , t ) − cpx yx + h V P ( x, y, t ) + ( y + 1) H y P ( x, y + 1 , t ) − yH y P ( x, y, t ) , (2)where V is the volume of the site in which the processesoccur. Using the van Kampen system’s size expansion[38], the master equation (2) can be written in the fol-lowing operator notation: ∂P ( x, y, t ) ∂t = [ V I ( E − x −
1) + ( E − x − r x x (1 − xK x V )+ ( E x − px yx + h V + H x ( E x − x + ( E − y − r y y (1 − yK y V )+ ( E − y − cpx yx + h V + H y ( E y − y ] P ( x, y, t ) , (3)where E x P ( x, y, t ) = P ( x + 1 , y, t ) = ρ ( x + 1 V , yV , t )= ρ (ˆ c x + 1 V , ˆ c y , t )= ρ (ˆ c x , ˆ c y , t ) + 1 V ∂ρ∂ ˆ c x + 12 V ∂ ρ∂ ˆ c x + ... = (1 + 1 V ∂∂ ˆ c x + 12 V ∂ ∂ ˆ c x ) ρ (ˆ c x , ˆ c y , t ) = E x ρ (ˆ c x , ˆ c y , t )with E x = (1 + V ∂∂ ˆ c x + V ∂ ∂ ˆ c x + ... ), ˆ c x = xV , andˆ c y = yV . Similarly we can compute E − x ρ (ˆ c x , ˆ c y , t ), where E − x = (1 − V ∂∂ ˆ c x + V ∂ ∂ ˆ c x + ... ), and E y ρ (ˆ c x , ˆ c y , t ) and E − y ρ (ˆ c x , ˆ c y , t ). Using the above expansions, (3) can bewritten as: ∂P ( x, y, t ) ∂t = [ ∂∂ ˆ c x ( − I − V r x x (1 − xK x V )+ 1 V px yx + h V + 1 V H x x )+ ∂∂ ˆ c y ( − V r y y (1 − yK y V ) − V cpx yx + h V + 1 V H y y )+ 12 V { ∂ ∂ ˆ c x ( IV + r x x (1 − xK x V ) + px yx + h V + H x x )+ ∂ ∂ ˆ c y ( r y y (1 − yK y V )+ cpx yx + h V + H y y ) } ] P ( x, y, t ) . (4)Equation (4) can also be written in the following com-pact form as: ∂P ( x, y, t ) ∂t = − ∂∂a .v ( a ) P ( x, y, t ) + ∂ T ∂a .B. ∂∂a P ( x, y, t ) , which is the Fokker-Planck equation corresponding to thesystem (1), where v ( a ) = " I + V r x x (1 − xK x V ) − V px yx + h V − V H x x V r y y (1 − yK y V ) + V cpx yx + h V − V H y y , ∂∂a = " ∂∂ ˆ c x ∂∂ ˆ c y , and B = (cid:20) b b (cid:21) with b = IV + r x x (1 − xK x V ) + px yx + h V + H x x and b = r y y (1 − yK y V ) + cpx yx + h V + H y y .
2. Presence of environmental noise
Ecosystems are often influenced by environmental fac-tors such as, climate change, weather conditions, etc.Stochastic fluctuations in population due to such exter-nal variations are termed as environmental noise, and itis known to be one of the major sources affecting theresilience of a system [39]. To study the influence of en-vironmental noise we formulate a stochastic model of (1)with multiplicative noise. The stochastic model is givenas: dXdt = F ( X ) + σXξ t , (5)where X = [ x, y ] T represents prey-predator fish vector. F ( X ) is the interaction between prey and predator fishgiven by the right hand side of (1). σ = [ σ x , σ y ] T is noiseintensity, and ξ t is the Gaussian white noise with zeromean and unit variance. The Fokker-Planck equationcorresponding to (5) can be written in a straightforwardmanner. The final expression of the Fokker-Planck equa-tion takes the following form: ∂P ( x, y, t ) ∂t = − ∂∂a .v ( a ) P ( x, y, t ) + ∂ T ∂a .B. ∂∂a P ( x, y, t ) , where v ( a ) = " I + V r x x (1 − xK x V ) − V px yx + h V − V H x x V r y y (1 − yK y V ) + V cpx yx + h V − V H y y ,∂∂a = " ∂∂ ˆ c x ∂∂ ˆ c y , and B = (cid:20) ˆ b
00 ˆ b (cid:21) , with ˆ b = σ x x and ˆ b = σ y y . TABLE I. Different birth and death processes, change of state vectors, gain and loss probabilities, and their propensity function,associated with the prey-predator fish model (1). V is the system’s volume where all the reactions occur. The symbols (+1) and( −
1) in the column of state vectors represent birth and death processes of the respective chemical species. Here, P stands forthe grand probability function. φ and z are empty state and dummy variables, respectively.Sl.No. Elementaryevents Reaction type Beforereaction Afterreaction Gain probability Loss probability Propensityfunction ( a n )1. φ → x Growth ofprey x − y xy V IP ( x − , y ) V IP ( x, y ) V I x + φ ⇋ x + x Logisticgrowthof prey x − y xy r x ( x − − x − K x V ) × P ( x − , y ) r x x (1 − xK x V ) × P ( x, y ) r x x (1 − xk x V )3. 2 x + y → zz → x Hollingtype IIIresponse forpredation x + 1 y xy p ( x +1) y ( x +1) + h V P ( x + 1 , y ) px yx + h V P ( x, y ) px yx + h V x → φ Mortality ofprey due toharvesting x + 1 y xy ( x + 1) H x P ( x + 1 , y ) xH x P ( x, y ) xH x y + φ ⇋ y + y Logisticgrowthof predator xy − xy r y ( y − − y − K y V ) × P ( x, y − r y y (1 − yK y V ) × P ( x, y ) r y y (1 − yK y V )6. 2 x + y → zz → y + y HollingType IIIprocess ofpredator xy − xy cpx ( y − x + h V P ( x, y − cpx yx + h V P ( x, y ) cpx yx + h V y → φ Mortality ofpredator dueto Harvesting xy + 1 xy ( y + 1) H y P ( x, y + 1) yH y P ( x, y ) yH y C. Stochastic simulation methods
It is a well known fact that presence of noise in abistable system can trigger a tipping [5], which will beevident from its trajectories. In the presence of demo-graphic noise, stochastic trajectories can be generatedby solving the master equation (2). We have used theMonte-Carlo simulation [40, 41] to numerically solve themaster equation (2), and finally generate stochastic tra-jectories. The simulation was performed using the Gille-spie algorithm [40]. The algorithm works in the followingway: it first calculate the propensity functions ( a n ) in Ta-ble I and then draw two uniform random numbers ( r and r ) from the interval (0 , r will be used to update the time ( t ) for the next processstochastically, where τ = a ln( r ), and a = P nj =1 a j .The other number r is used to identify the index n ofthe occurring process and which is given by the smallest integer satisfying: n − X j =1 a j < r a ≤ n X j =1 a j . The system states are updated X ( t + τ ) = X ( t ) + ∆ X ν n ,then the simulation proceeds to the next occurring time.∆ X ν n represents change in species vector in the stochas-tic time interval τ .To calculate trajectories of the model in the presence ofenvironmental noise, we performed stochastic simulationsof (5) in MATLAB (R2018b) using the Euler-Maruyamamethod with an integration step size of 0.01. D. Probabilistic potentials
Here, we try to understand the influence of demo-graphic noise by calculating probabilistic potentials fordifferent prey harvesting rates ( H x ). Probabilistic po-tentials are calculated from numerically simulated tra-jectories. These trajectories are computed for 4 × different initial conditions, which are uniformly drawnfrom [0, 6] × [1 .
5, 2 . x ss and y ss are respective steady state densities of preyand predator fishes, M ssxy is the frequency count to visit asteady state ( x ss , y ss ) and N is the total number of vis-its. Then probability is defined as P ss ( x ss , y ss ) = M ssxy N .The negative logarithm of this probability is the desiredstochastic potential. E. Mean first-passage time
In a bistable region or a bistable potential, two steadystates are separated by an unstable steady state definingthe barrier state or the potential barrier. In the pres-ence of noise, often exit from a potential well is possible.In this regard, one can calculate the first passage time(FPT). The FPT is calculated by considering a family ofinitial points, ( x , y ) that starts somewhere in a poten-tial well of volume V in space, bounded by a surface ∂V .Then the FPT is the time when the point exits V forthe first time. The mean first-passage time (MFPT) isthe average time of the FPT. The permanence of steadystates in a noise-driven system is often characterized bythe MFPT [42]. Ideally, the motion of the family of ini-tial points satisfies both of the master and Fokker-Planckequations and we focus only those points that have notleft V by time t . We remove all the paths that havecrossed the boundary of V before time t , which is usu-ally done by placing an absorbing boundary conditionon ∂V . We impose the absorbing wall at the potentialbarrier of the bistable system. We perform Gillespie sim-ulation for a family of initial points, which allows us tofind the time taken to cross the absorbing wall. Here, themean time is calculated over 5 × realizations, and allthe initial conditions are taken uniformly from the ranges[0 , × [1 . , . F. Early warning signals
The presence of noise in a bistable system can drivetipping, and it can be predicted by critical slowing downbased EWSs. We calculate a few generic EWSs, like vari-ance, lag-1 autocorrelation (AR(1)), and conditional het-eroskedasticity to forewarn upcoming transitions. Vari-ance ( σ ) is the variability of a random variable x aroundthe mean ( µ ) and is given by σ = N P Nn =1 [ x ( t ) − µ ] ,where N is number of time step. We also calculatethe autocorrelation at lag-1, which is given by ρ = E { [ x ( t +1) − µ ][ x ( t ) − µ ] } σ , where E is expectation value of arandom variable, x ( t ) is the value of the state variableat time t . To get robust predictions, we also calculateconditional heteroskedasticity following methods in [43].It uses the Lagrange multiplier test to calculate the con- ACE BDF
MonostabilityBistability
FIG. 1. Bifurcation diagrams of the model (1). Co-dimensionone bifurcation diagram of (A) the prey fish ( x ), and (B) thepredator fish ( y ), for changes in H x . (C-F) Co-dimension twobifurcation diagrams in the following c − H x , H y − H x p − H x and r x − H x planes, respectively. In (A) and (B), solid anddashed curves denote stable and unstable states. In (C- F),the darker regions denote bistability, and other regions denotemonostability. Dashed lines in two-parameter bifurcation di-agrams denote the values used for numerical simulation. ditional heteroskedasticity. III. RESULTS
Depending on the parameter values and the HollingType-III functional response, the model (1) exhibitsmonostability as well as bistability (see Fig. 1). The sys-tem evolves from a high-density state to a low-densitystate of fish via a hysteresis loop with an increase in theharvesting rate. With a low to moderate harvest, there isa high positive equilibrium. With overharvest, the highpositive equilibrium suddenly collapses into a low posi-tive equilibrium. Unless otherwise stated, values of theparameters used throughout this study are: I = 0 . r x = 0 . K x = 6, h = 0 . r y = 0 . K y = 4, c = 0 . p = 0 .
15, and H y = 0 . A. Stochastic potentials
Stochastic potentials of the system (1) for different H x are depicted in Fig. 2. Potentials show a clear effectof stochasticity on the model’s behavior. Increasing theharvesting rate changes the probabilistic potential, eachfor low and high values of H x ; there exists one potentialwell. For an intermediate value of H x , there exist two FIG. 2. Stochastic potential landscapes for different values of the prey harvesting rate ( H x ), obtained by solving the masterequation (2). Stochastic potential for: (A) monostable high density fish state at H x = 0 .
03, (D) bistable high density and lowdensity fish states for H x = 0 .
2, and (G) monostable low density fish state for H x = 0 .
45. The blowup diagrams (B, C, E, F)magnify regions in potential wells. potential wells. The lowest value depicted in the poten-tial color-bar corresponds to the existence of a deep well,i.e., the likelihood of a return to the steady state, afterperturbation, is more in this well. It is evident that, inthe monostable regions, the high-density state’s well isflatter and shallower (Figs. 2A and 2B), and in compar-ison to that the low-density well is deeper (Figs. 2F and2G). However, in the bistable region, the potential wellfor the high-density state is deeper than the low-densitystate (Figs. 2C, 2D and 2G). This deeper well representshigher resilience and indicates that a population near thehigh-density state is more likely to stay there unless it islargely perturbed in the bistable region.
B. Mean first-passage time
The effects of changing H x on MFPT under intrin-sic and extrinsic stochastic fluctuations are shown inFigs. 3A and 3B. The reciprocal of MFPT is the rateof arrival at the potential barrier. The mean time takenfor both the noises to shift from upper to lower steadystate decreases with an increase in the harvesting rateof prey, whereas the mean time decreases with decreas-ing harvesting rate to traverse from a lower to an upperstable state. In other words, the rate of transition from the upper steady state to the lower steady state increaseswith an increase in the prey’s harvesting rate. However,the transition from lower to upper steady state happensquite rarely since the rate is slow. It is also evident thatMFPT for the high-density state is much higher thanthe low-density state, indicating higher resilience of thehigh-density state. The crossing point of MFPTs for twodifferent noises are different. In the case of environmentalnoise, the point is earlier in comparison to demographicnoise (Figs. 3A and 3B). This establishes that MFPTchanges with noise type and the fluctuations observedunder environmental stochasticity are more than that ofdemographic stochasticity. C. Early warning signals
Here we explore the efficacy of EWSs in anticipatingregime shifts by analyzing stochastic time series data(see Fig. 4) for both the demographic and environmentalnoise. To calculate variance and AR(1), we have selecteda time series segment ( Figs. 4A and 4B) before a transi-tion. Indicators of an impending transition are affectedby non-stationarities in the mean of the time series, es-pecially the metrics calculated in a moving windows [28].So, it is always appropriate to remove non-stationarities M F P T M F P T BA FIG. 3. Mean first-passage time with an increase in H x : (A)for demographic noise, and (B) for environmental noise. Solidcurve shows the MFPT taken by the system to switch from theupper steady state to the lower steady state. Similarly dottedcurve shows MFPT taken by the system to switch from thelower steady state to the upper state. Environmental noiseintensity is taken as σ = 0 . as it imposes a formidable correlation structure on timeseries. Hence, for our analysis, all metrics calculatedwithin rolling windows have removed high frequencies us-ing the Gaussian detrending with a bandwidth. Then wecalculate variance (see Figs. 5C and 5K) and AR(1) (seeFigs. 5D and 5L) using rolling windows of 50% and 60%of the size of time series, respectively.We see that under both kinds of stochasticity (de-mographic and environmental), the variance is robustand successfully estimates an impending transition (seeFigs. 5C and 5K). However, we find that AR(1) is onlysuccessful in predicting regime shifts in the case of envi-ronmental stochasticity as depicted in Fig. 5L, but doesnot work well for demographic noise (see Fig. 5D). Thisasserts that the commonly used CSD indicators are notalways robust in predicting an upcoming transition.Therefore we use another indicator called conditionalheteroskedasticity, which minimizes false-positive signalsfrom the time series. It is the conditional error variancein time series, i.e., error variance at a time step is depen-dent or conditional on the error variance of the previoustime step. We consider the time series data for demo-graphic and environmental noise (see Figs. 6A and 6D)to calculate conditional heteroskedasticity. We took arolling window of size 10 of the time series while cal-culating the conditional heteroskedasticity. To removestationarity from the data set, we apply the Gaussian x A Time H x x B Time H x PSfrag replacementsTime λ Budworm DensityWhite NoiseRed − shifted Noise EDCBAFGHI h T x stl → x stu ih T x stu → x stl i τσ = 0 . b > ,c > b > ,c < b < , < c < b / f ( x ) O R O R FIG. 4. Critical transitions in the presence of: (A) demo-graphic noise, and (B) environmental noise, with increasingharvesting rate H x . Cyan solid lines and red dashed lines rep-resent stable and unstable steady states, respectively. Blackdashed vertical line represents the position upto where thetimes series data has taken for statistical analysis. Noise in-tensity for the environmental noise is taken as σ = 0 . kernel fitting to the time series and calculate its residualsquare. Here squared residuals at time step t is plottedwith the previous time step (see Figs. 6B and 6E). Thepositive correlation between these two squared residu-als suggests that both time series are having conditionalheteroskedasticity. Hence, the cumulative number of sig-nificant Lagrange multiplier tests ( C ) is carried out onboth data sets, as shown in Figs. 6C and 6F. We ob-serve an increasing trend in C for both the time series,thus indicating that both the time series are conditionallyheteroskedastic and approach a critical transition. D. Collective influence of predation and harvestingon the dynamics of fisheries
Fish stocks are greatly influenced by two major fac-tors – predation and harvesting. Therefore, illustratingthe criticality of these two processes is significant in reg-ulating fish populations. Here we explore the qualitativedynamics of the two species model (1) with varying thepredation rate p and the harvesting rate H x (see Fig. 7).As can be seen from Figs. 7A and 7B, an increase in H x -0.11-0.050.01 0 300 600 900 A R ( ) Time D V a r i an c e C -0.600.6 R e s i dua l s B x A O cc u rr en c e Kendall τ estimates G AR(1) B and w i d t h Rolling window H K enda ll τ
25 37 50 62 755255075100 -1-0.5 0 0.5 1 B and w i d t h Rolling window F K enda ll τ
25 37 50 62 755255075100 -1-0.5 0 0.5 1 0 130 -1 -0.5 0 0.5 1 O cc u rr en c e Kendall τ estimates E Variance A R ( ) Time L V a r i an c e K -0.500.5 R e s i dua l s J x I O cc u rr en c e Kendall τ estimates O AR(1) B and w i d t h Rolling window P K enda ll τ
25 37 50 62 755255075100 -1-0.5 0 0.5 1 B and w i d t h Rolling window N K enda ll τ
25 37 50 62 755255075100 -1-0.5 0 0.5 1 0 60 -1 -0.5 0 0.5 1 O cc u rr en c e Kendall τ estimates M Variance
PSfrag replacementsTime λ Budworm DensityWhite NoiseRed − shifted Noise EDCBAFGHI h T x stl → x stu ih T x stu → x stl i τσ = 0 . b > , c > b > , c < b < , < c < b / f ( x ) O R O R FIG. 5. Early warning signals were calculated from the time series data for demographic and environmental noises. (A, I)Simulated time series data before a transition. (B, J) Residual time series after applying a Gaussian filter (red curve in (A, I)is the trend used for filtering). (C, D) Variance and AR(1) are calculated using a moving window size of 50% and bandwidth10. Similarly, (K, L) represents variance, and AR(1) calculated using a moving window size 60% and bandwidth 75. (E, G, M,O) Distribution of the Kendall- τ estimates. (F, H, N, P) Sensitivity analysis for the combination of rolling window size andbandwidth. eventually leads to the collapse of prey fish. However,the occurrence of fluctuations in prey densities is gov-erned by the predation level. In Fig 7A, for p = 0 . H x = 0 .
47, whereasfor high predation value i.e. p = 0 .
35 the fluctuationsare evident at lower value of harvesting, i.e. H x = 0 . p and lower values of prey harvesting rate H x ,we observe saddle-node bifurcation (see Fig. 7C). Withan increase in H x and lower p value, near 0 .
29, we getHopf bifurcation. One may note that at lower p value,the system experiences monostable state. Further, fromFig. 7D, we find that a high predator harvesting rate H y and high H x exhibits cyclic dynamics. However, for lowvalues of H y , cyclic dynamics no longer exist and exhibitssaddle-node bifurcation. E. Influence of social norms
Due to over-harvesting, fisheries are prone to suddencollapse, which is evident from our previous analysis (seeFig. 1 and Fig. 5). Therefore, for the future sustainabilityof fisheries, fishing management strategies are required tostrengthen such systems’ resilience. Here, we incorporatesocial norms in the two species model (1). We apply theconcept of opinion dynamics to the prey harvesting rate H x , which now takes the form H x (1 − s ). Here, s isthe proportion of the population that adopts the opinionof conservation. The opinion dynamics is described bythe concept of imitation dynamics of evolutionary gametheory [33]. It captures human behavior, especially howhumans like to imitate successful strategies [34]. Here, weconsider the opinions of either ”adapting conservation”or ”not adapting conservation”, and each one is associ-ated with a payoff. Each individual can sample anotherindividual with a sampling rate. In particular, if the sam-pled person has other opinions with a higher payoff, theindividual switches to the sampled person’s opinion withthe expected gain in the payoff. The following is the re- X A t S qua r ed r e s i dua l s a t t + B C C D t E F PSfrag replacementsTime λ Budworm DensityWhite NoiseRed − shifted Noise EDCBAFGHI h T x stl → x stu ih T x stu → x stl i τσ = 0 . b > , c > b > , c < b < , < c < b / f ( x ) O R O R FIG. 6. Time series before critical transitions for (A) de-mographic noise, and (D) environmental noise. (B, E) Cor-relation between squared residuals at time step t with thenext time step ( t + 1) from an autoregressive lag-1 modeland slanted lines present the linear regression lines. (C, F)Cumulative Lagrange multiplier test ( C ) for conditional het-eroskedasticity. sulting coupled socio-ecological model with social norms: dxdt = I + r x x (1 − xK x ) − pyx h + x − H x (1 − s ) x, (6a) dydt = r y y (1 − yK y ) + cpyx h + x − H y y, (6b) dsdt = ks (1 − s ) (cid:18) d (2 s − − w + 1 x + c (cid:19) , (6c)where k and d are the sampling rate and strength of theinjunctive social norms, respectively. w is the total costof conservation, and c is a control parameter that con-trols the prey fish density with the payoff of conservingfish. In particular, for a small value of x , x + c moti-vates fish conservation, as the utility to conserve preyfish increases when fish densities become very low. Allthe other parameters are the same as in (1).Here, we examine the qualitative dynamics of themodel (6), using one parameter bifurcation diagrams byvarying H x (see Fig. 8). The coupled model exhibitstranscritical and saddle-node bifurcations (see Fig. 8A).We can say that over-harvesting is a significant cause ofthe collapse of fish dynamics (see Fig. 1). However, avery different outcome occurs when small conservationcost ( w = 0 .
3) and weaker social norms ( d = 0 .
1) isapplied. In the vicinity of the tipping point, where H x p H x C MS HBSN x H x A p=0.29 H y H x D MS HBSN x H x B p=0.35 FIG. 7. Codimension-one bifurcation diagrams of the preyfish x with variations in H x , for different values of p : (A) p = 0 .
29, and (B) p = 0 .
35, with H y = 0 .
275 (along thedashed line shown in subfigure (D)). (C-D) Codimension-twobifurcation diagrams in the H x − p plane and H x − H y plane,respectively. In (A-B), red solid curves denote stable steadystate, and black dashed curves denote unstable steady state.Green and blue circles are the stable and unstable limit cycles,respectively. In (C-D), MS, SN, and HB represent monos-table, saddle-node bifurcation, and Hopf bifurcation curves,respectively. AB FIG. 8. (A-B) Codimension-one bifurcation diagrams for theeqn 6. (A) One parameter bifurcation diagram of the preyfish ( x ) at different prey harvesting ( H x ). (B) One parame-ter bifurcation diagram of the social norms at different preyharvesting ( H x ). Blue lines and red dotted lines are denotingstable attractor and unstable attractor respectively. Otherparameters are I = 0 . c = 0 . H y = 0 . r x = 0 . K x = 6, h = 0 . r y = 0 . K y = 4, p = 0 . k = 0 . w = 0 . c = 0 . d = 0 . . s > IV. DISCUSSION
Increasing anthropogenic pressure, in the form of in-creasing population density [48], invoke profound globalchange that has many future repercussions, specificallyincreasing harvesting rates, which may flip some lakesor oceans from a stable state to an uncertain condition[14, 49, 50]. In this paper, we investigate the response ofan ecological model towards harvesting pressures, underthe influence of demographic and environmental stochas-ticity. We proceed by analyzing the behavior of a two-species fish model (1). We find that the model involves ahysteresis loop (see Fig. 1), which primarily instincts thepossibility of sudden critical transitions. Here we haveused the master equation, which follows the Markov pro-cess, to study the effects of demographic stochasticityon the model. For different parameter values, stochasticpotentials are calculated to examine the influence of de-mographic noise on the resilience of steady states whosepersistence is supported by MFPT. MFPT for the con-sidered model reveals that the upper steady state is moreresilient than the lower steady state. Further, for envi-ronmental stochasticity the system experiences more fluc-tuations when compared to demographic stochasticity.In the case of environmental noise, we calculate the sys-tem trajectory by applying the Euler-Maruyama method.Our results show that in the presence of demographic andenvironmental stochasticity, an increase in the harvest-ing rate of prey species induces a critical transition froma high fish density state to a low fish density state.Moreover, we calculate a few generic early warning sig-nals to forewarn the chance of critical transitions. Pop-ulation extinction has always been a major concern, butanticipating species extinction remained challenging [51].This lead to a quest for robust indicators of critical tran-sitions in ecology [27, 52]. We found that the indicatorAR(1) could successfully forewarn an impending transi-tion for environmental stochasticity but fails to provideany warning in the case of demographic noise. However, variance and conditional heteroskedasticity work well inpredicting transitions for both types of stochasticity. Fur-ther, we also emphasized the importance of consideringmanagement steps for the better future of the fisheries.The extended socio-ecological model gives the qualita-tive dynamics of prey fish densities and social norms andtheir relationship with harvesting ( H x ). We observe thatthe application of social norms leads to the growth ofprey species, which otherwise would have resulted in aninevitable collapse of the system.The anticipation of critical transitions could preventsudden catastrophes and thus enhance the stability of asystem. Furthermore, using various early warning indi-cators that could predict an upcoming shift would be in-valuable. For instance, predicting regime shifts in marineecosystems using EWS could impede the switch to lowfish density and result in the persistence of the marinepopulation. However, theories suggest that it is necessaryto analyze specific situations under which generic earlywarning indicators fail to work effectively. Our resultsfirmly outline the need to develop robust early warningindicators with a good understanding of which of themmight be most convenient in terms of perceptibility andreliability.On a final note, our work illustrates the influenceof natural and anthropogenic factors on an ecologicalsystem. Further, understanding critical transitions inhigher-dimensional systems, as in network models, wouldbe a useful research area. Moreover, developing modelsto explore rate-induced tipping and considering the ef-fects of different noise types on the predictability of earlywarning signals is an important future direction. Nev-ertheless, predicting the distance to critical transitionsremains an important research question, together withpredicting the system’s state after these sudden transi-tions. ACKNOWLEDGMENTS
S.S. acknowledges the financial support from DST,India under the scheme DST-Inspire [Grant No.:IF160459]. P.S.D. acknowledges financial support fromthe Science & Engineering Research Board (SERB),Govt. of India [Grant No.: CRG/2019/002402]. [1] O. N. Bjørnstad and B. T. Grenfell, Science , 638(2001).[2] T. Shinbrot and F. J. Muzzio, Nature , 251 (2001).[3] T. Coulson, P. Rohani, and M. Pascual, Trends in Ecol-ogy & Evolution , 359 (2004).[4] A. J. Black and A. J. McKane, Trends in Ecology &Evolution , 337 (2012).[5] M. Scheffer, S. Carpenter, J. A. Foley, C. Folke, andB. Walker, Nature , 591 (2001). [6] M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin,S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes,M. Rietkerk, and G. Sugihara, Nature , 53 (2009).[7] K. Higgins, A. Hastings, J. N. Sarvela, and L. W. Bots-ford, Science , 1431 (1997).[8] R. A. Myers, G. Mertz, J. M. Bridson, and M. J. Brad-ford, Canadian Journal of Fisheries and Aquatic Sciences , 2355 (1998).[9] O. N. Bjørnstad, J.-M. Fromentin, N. C. Stenseth, and J. Gjøsæter, Proceedings of the National Academy of Sci-ences U.S.A. , 5066 (1999).[10] P. Lundberg, E. Ranta, J. Ripa, and V. Kaitala, Trendsin Ecology & Evolution , 460 (2000).[11] J.-M. Fromentin, R. A. Myers, O. N. Bjørnstad, N. C.Stenseth, J. Gjøsæter, and H. Christie, Ecology , 567(2001).[12] P. M. Vitousek, H. A. Mooney, J. Lubchenco, and J. M.Melillo, Science , 494 (1997).[13] J. Lubchenco, Science , 491 (1998).[14] J. M. Fryxell, R. Hilborn, C. Bieg, K. Turgeon,A. Caskenette, and K. S. McCann, Proceedings of theNational Academy of Sciences U.S.A. , 12333 (2017).[15] R. Myers, N. Barrowman, J. Hutchings, and A. Rosen-berg, Science , 1106 (1995).[16] J. K. Baum, R. A. Myers, D. G. Kehler, B. Worm, S. J.Harley, and P. A. Doherty, Science , 389 (2003).[17] J. Gulland, in Dynamics of Populations (Centre forAgricultural Publishing and Documentation Wagenin-gen, 1971) pp. 450–467.[18] D. H. Cushing and D. Cushing,
The provident sea (Cam-bridge University Press, 1988).[19] J. B. Jackson, M. X. Kirby, W. H. Berger, K. A. Bjorn-dal, L. W. Botsford, B. J. Bourque, R. H. Bradbury,R. Cooke, J. Erlandson, J. A. Estes, et al. , Science ,629 (2001).[20] C. F. Clements, M. A. McCarthy, and J. L. Blanchard,Nature Communications , 1 (2019).[21] C. M¨ollmann and R. Diekmann, in Advances in EcologicalResearch , Vol. 47 (Elsevier, 2012) pp. 303–347.[22] T. Oguz and D. Gilbert, Deep Sea Research Part I:Oceanographic Research Papers , 220 (2007).[23] A. C. Kraberg, N. Wasmund, J. Vanaverbeke,D. Schiedek, K. H. Wiltshire, and N. Mieszkowska, Ma-rine Pollution Bulletin , 7 (2011).[24] Y. Jiao, Reviews in Fish Biology and Fisheries , 177(2009).[25] A. J. Pershing, K. E. Mills, N. R. Record, K. Stamieszkin,K. V. Wurtzell, C. J. Byron, D. Fitzpatrick, W. J. Golet,and E. Koob, Philosophical Transactions of the RoyalSociety B: Biological Sciences , 20130265 (2015).[26] C. F. Clements, J. L. Blanchard, K. L. Nash, M. A. Hin-dell, and A. Ozgul, Nature Ecology & Evolution , 0188(2017).[27] J. M. Drake and B. D. Griffen, Nature , 456 (2010).[28] V. Dakos, S. R. Carpenter, W. A. Brock, A. M. Ellison,V. Guttal, A. R. Ives, S. Kefi, V. Livina, D. A. Seekell,E. H. van Nes, et al. , PloS One , e41010 (2012).[29] S. Sarkar, S. K. Sinha, H. Levine, M. K. Jolly, and P. S.Dutta, Proceedings of the National Academy of Sciences U.S.A. , 26343 (2019).[30] L. Dai, D. Vorselen, K. S. Korolev, and J. Gore, Science , 1175 (2012).[31] C. Boettiger and A. Hastings, Journal of the Royal Soci-ety Interface , 2527 (2012).[32] P. S. Dutta, Y. Sharma, and K. C. Abbott, Oikos ,1251 (2018).[33] S. J. Lade, A. Tavoni, S. A. Levin, and M. Schl¨uter,Theoretical Ecology , 359 (2013).[34] J. Hofbauer and K. Sigmund, Evolutionary games andpopulation dynamics (Cambridge university press, 1998).[35] S. R. Carpenter, W. A. Brock, G. J. Hansen, J. F.Hansen, J. M. Hennessy, D. A. Isermann, E. J. Peder-sen, K. M. Perales, A. L. Rypel, G. G. Sass, et al. , Fishand Fisheries , 1150 (2017).[36] R. M. May, Stability and complexity in model ecosystems ,Vol. 1 (Princeton university press, 2019).[37] R. Lande, The American Naturalist , 911 (1993).[38] N. G. Van Kampen, Adv. Chem. Phys , 245 (1976).[39] C. Boettiger, Ecology Letters , 1255 (2018).[40] D. T. Gillespie, Journal of Computational Physics ,403 (1976).[41] N. G. Van Kampen, Stochastic processes in physics andchemistry , Vol. 1 (Elsevier, 1992).[42] C. W. Gardiner et al. , Handbook of stochastic methods ,Vol. 3 (Springer Berlin, 1985).[43] D. A. Seekell, S. R. Carpenter, and M. L. Pace, TheAmerican Naturalist , 442 (2011).[44] J. J. Lever, E. H. van Nes, M. Scheffer, and J. Bas-compte, Ecology Letters , 350 (2014).[45] J. A. Estes, D. O. Duggins, and G. B. Rathbun, Con-servation Biology , 252 (1989).[46] T. P. Hughes, Science , 1547 (1994).[47] H. Lessios, Annual Review of Ecology and Systematics , 371 (1988).[48] S. Barrett, A. Dasgupta, P. Dasgupta, W. N. Adger,J. Anderies, J. van den Bergh, C. Bledsoe, J. Bongaarts,S. Carpenter, F. S. Chapin, et al. , Proceedings of theNational Academy of Sciences USA , 6300 (2020).[49] C. M. Roberts and J. P. Hawkins, Trends in Ecology &Evolution , 241 (1999).[50] Y. Sadovy, Journal of Fish Biology , 90 (2001).[51] D. Ludwig, Ecology , 298 (1999).[52] V. Dakos, M. Scheffer, E. H. van Nes, V. Brovkin,V. Petoukhov, and H. Held, Proceedings of the NationalAcademy of Sciences U.S.A. , 14308 (2008).[53] V. Dakos, E. H. Van Nes, P. D’Odorico, and M. Scheffer,Ecology93