Integrating economic dynamics into ecological networks: The case of fishery sustainability

See allHide authors and affiliations

Science Advances  04 Nov 2020:
Vol. 6, no. 45, eaaz4891
DOI: 10.1126/sciadv.aaz4891


Understanding anthropogenic impacts on ecosystems requires investigating feedback processes between ecological and economic dynamics. While network ecology has advanced our understanding of large-scale communities, it has not robustly coupled economic drivers of anthropogenic impact to ecological outcomes. Leveraging allometric trophic network models, we study such integrated economic-ecological dynamics in the case of fishery sustainability. We incorporate economic drivers of fishing effort into food-web network models, evaluating the dynamics of thousands of single-species fisheries across hundreds of simulated food webs under fixed-effort and open-access management strategies. Analyzing simulation results reveals that harvesting species with high population biomass can initially support fishery persistence but threatens long-term economic and ecological sustainability by indirectly inducing extinction cascades in non-harvested species. This dynamic is exacerbated in open-access fisheries where profit-driven growth in fishing effort increases perturbation strength. Our results demonstrate how network theory provides necessary ecological context when considering the sustainability of economically dynamic fishing effort.


The advent of network theory in ecology and environmental studies has greatly advanced the study of ecological dynamics and complexity (1, 2). These advances have also translated into a growing knowledge of how human-caused disturbances can create far-reaching ecological impacts through indirect effects (36). Often, however, the disturbances in these network studies have been studied as a one-time event or a constant external rate of change, separate from the dynamic elements of the ecological network (3, 6, 7). We argue that developing both sustainable management practices and a fuller understanding of Anthropocene ecological dynamics require recognition that much of the ecological impact of anthropogenic activity is determined by an integrated feedback process between ecological dynamics and socioeconomic conditions as a coupled natural-human system (8, 9). In short, the scope of ecological networks should contain humans as dynamic elements when necessary (10, 11). Here, we expand ecological network theory by incorporating economic dynamics into food web network models to evaluate the coupled natural-human dynamics affecting sustainability in the case of fisheries.

Fisheries are an example of natural-human integrated network systems that are particularly important given their critical role in the economic stability and food security of billions of people (12). Both the available yield and mobilized fishing effort in any fishery are products of a complex set of interacting socioeconomic and ecological factors (13). Understanding the dynamics and consequences of these interactions is particularly pressing, as current management strategies have produced an overexploitation crisis in a multitude of fisheries across the globe (14, 15), threatening both aquatic biodiversity (16) and the aforementioned food security (12). Traditionally, fishery sustainability goals have been implemented into policy on a single-species basis as optimal harvesting strategies (12). While economic considerations have been well integrated into optimization studies (17, 18), these strategies conventionally consider harvested species in isolation from the broader food web and model their dynamics with an idealized growth curve (13, 19, 20). This has limited the ability of research based on these assumptions to address the indirect effects of harvesting on other species within food webs, as seen in network-based approaches (21). Specifically, the population variability induced directly through fishing effort on harvested species (22) can transfer to non-harvested species through trophic interactions (5, 23). This variability has been linked to perturbations, causing reductions in aquatic biodiversity (24, 25) and ecosystem function (26). These reductions cycle back to further affect harvested species (5, 27) and the potential economic returns harvested species provide (14, 28). Any effect on economic returns can affect future fishing effort and the corresponding effect on harvested ecosystems, creating a bio-economic feedback loop. While this bio-economic loop is understood in concept, limited work exists on its actual dynamics, especially in the context of broader ecological networks (29). Until more is known about the complex interactions between ecology and economic factors that drive the bio-economic loop in fisheries, policy-makers risk attempting to optimize a process that we poorly understand.

In accordance with these concerns, we applied network theory to incorporate the ecological complexity of species trophic interactions into the model. This differs from past modeling work using approaches like Ecopath with Ecosim (EwE) (30), which rely on extensive lists of estimated system- and species-specific parameters to manage specific systems rather than study general patterns in fundamental dynamics. This specificity has been essential to the management potential of EwE and allowed it to function as an effective tool in numerous fisheries. However, estimating the full parameter set required for EwE is data intensive (31) and potentially difficult to empirically or theoretically verify (32). This can lead to uncertainty in parameter estimates (13) and prompt some modelers to use “default” values (32). Given that model outcomes can be greatly altered, with multiplicative effects, by erroneous parameter estimates (33, 34), models must be finely tuned to their specific system. EwE includes mechanisms to improve parameter estimates against data from the fisheries in question, but this further delimits investigations to specific systems studied and less toward understanding more universal underlying mechanisms driving the bio-economic feedback loop. While this precision offered by system-specific investigations is undoubtedly important, complementary investigations into the underlying generality of the dynamics seen across fisheries are a necessary part of a multimodel strategy (see appendix S1) (35, 36).

The network approach that we use allows for the development of more fundamental, widely applicable theory because network models can be run across an array of ecological topologies, both empirical (5) and realistically generated (37), by parameterizing metabolic rates and species interactions through allometric scaling (38). Experimental results across these collections of webs can then be searched for patterns and explanatory drivers. For the purposes of this study, we generated realistic food web networks using the niche model (39), each with an initial 30 interacting trophic species (see Materials and Methods). Species available for harvest are labeled “fish” for ease of description (Fig. 1A). Ecological dynamics in each generated web were governed by a series of ordinary differential equations and parameterized through allometrically scaled rates (38), creating allometric trophic network (ATN) models (see Eqs. 1 to 10 in Materials and Methods). Allometric scaling provides a sound theoretical basis from which to realistically parameterize various iterations of trophic interactions and ecological processes based on metabolic theory and life history traits (e.g., trophic level and vertebrate/invertebrate) without exhaustive parameter sweeps (38, 40). Given the high dimensionality of food web models (e.g., coupled differential equations; see Materials and Methods) (30), mathematical analysis is intractable. Instead, general results are obtained through statistical analysis on simulation output. Last, this network approach then also provides a flexible framework that facilitates the integration of economic dynamics.

Fig. 1 Diagrams of the experimental design.

(A to C) Example evolution of a food web structure across complete simulation process with a fixed effort of 5. Node sizes are logarithmically scaled to the biomass at each point in time of the simulation. Edges between nodes represent trophic interactions with arrows indicating the consumer. The vertical axis indicates the trophic level of trophic species. (A) The example food web at initialization (t = 0). Colors of nodes indicate species that will survive or go extinct after initialization. (B) The example food web after initialization and before fishing effort (t = 4000). Colors of nodes indicate species that will survive or go extinct after fishing. Indices (i to iv) indicate the order of extinctions during simulation. (C) The final structure of the food web after the fixed harvesting effort (t = 8000). (D) Example time series from fixed-effort treatment simulations with effort set at 20. (E) Example time series from open-access treatment simulations.

We incorporated two economic models driving fishing effort into the ATN’s ecological network structure with fisheries functioning as an additional node in the networks during simulations (see Materials and Methods for the description of fishing effort). After an initialization period of 4000 time steps, roughly 11 years in model time (see Materials and Methods) (Fig. 1B), each “conserved” food web (see Materials and Methods) is subjected to two fishery treatments (Fig. 1C): (i) fixed effort and (ii) open access. The fixed-effort treatment uses fixed levels of fishing effort starting immediately after the initialization period that does not change within simulations (Fig. 1D). In the open-access treatment, on the other hand, fishing effort is unregulated (41). That is, effort adjusts in response to fishing profits, with effort growing or declining in response to positive or negative net profits, respectively (Fig. 1E). The dynamics of net profits are influenced by both yield and market price, with market price related to yield through a linear pricing model (Eqs. 11 to 13; see Materials and Methods). We illustrate these economic dynamics using a graphical example in fig. S1. The static effort levels of the fixed-effort treatment serves as a control to the variable dynamics of open-access fisheries, but fixed-effort levels also has real-world representations in some subsistence fisheries where strict permitting rules produce stable effort levels over the years (42). The open-access treatments simulate fisheries from their initialization and therefore start from a low initial effort at t = 4000. All open-access results presented in the main text pertain to initial effort levels of 1, while initial effort levels of 0.5 and 2 were also considered for sensitivity analysis (see Materials and Methods). Effects of different economic conditions are studied by parameter sweeps across levels of price sensitivity to yield (b), effort’s sensitivity to changes in profit (μ), and maximum price (a) paid for the harvested species. See Materials and Methods for more. We simulated single-species fisheries, labeling the single harvested fish in each simulation, H. The remaining non-harvested species in each simulation are then labeled, N-H. While single-species fisheries are not necessarily the dominant form of fishery in terms of quantity, they make up some of the largest fisheries globally and have a large economic and ecological influence (43, 44), making them a useful starting point for our network-based studies of the bio-economic feedback loop.

We use our dynamic model to evaluate the economic and ecological factors that determine (i) the impact of fishing on harvested and non-harvested species (ecological impacts, i.e., changes in abundance and persistence), (ii) the conditions for fishery “success” (i.e., a sustained nonzero fishing effort from start to finish), and (iii) the different ecological impacts of fisheries within fixed and open-access regimes (i.e., the persistence of harvested species and degree of species loss).


When we implement fishing through the fixed-effort treatment, higher effort levels increase the mortality of H, intuitively causing more biomass depletion (fig. S2), more extinctions (fig. S3), and quicker times to extinction (Fig. 2A) for the harvested species, H. Among the hundreds of ecological factors analyzed (see Materials and Methods), we found those H extinctions to be more prevalent at higher trophic levels (Fig. 2B and fig. S4), reflecting numerous empirical examples (4547). The best ecological predictor of H extinctions, the population biomass of H at the start of fishing (BH0; table S1), has a nonlinear effect (Fig. 2B and fig. S5A). Compared to populations with the lowest BH0, moderate increases in the starting population biomass of the harvested species decreased the prevalence of H extinctions because more abundant harvested populations are more resistant to the direct extraction-induced mortality. However, for all but the highest trophic levels, we saw that further increases in BH0 escalate extinction risk of H. This nonlinearity initially seems counterintuitive, as one could expect larger populations to better withstand harvesting pressure.

Fig. 2 Fixed-effort results.

(A) Scatterplot showing the effect of log(BH0) on time to H extinction. (B) Generalized additive model (GAM) binomial regression of the effects of log(BH0) on the probability of fishing resulting in H extinction across four trophic level groupings of H. Lines represent fit, and shaded regions show 95% confidence interval. Including log(BH0), trophic level of H, and fixed-effort level gives the binomial GAM results as follows: 51.1% of variance exampled at R2 = 0.57, UBRE (Un-Biased Risk Estimator) = −0.32. (C) GAM regressions showing the percent of the community lost based on log(BH0). Lines represent fit, and shaded regions show 95% confidence interval. Higher effort levels (e.g., E = 15, 20) do present a nonlinear trend at high BH0. This is due to the expedited H extinctions at high effort levels ending the active fishing disturbance on the rest of the community quicker than lower effort levels. See fig. S6 for more. (D) Scatterplot showing positive correlative relationship between the percent of community lost to N-H extinctions and the percent of H’s prey items lost. Colors represent fixed-effort values similar to (A), and the green line represents GAM regression fit with 95% confidence interval.

This nonlinearity occurs because fishing higher BH0 generally induces greater levels of variability in the rest of the populations within the food web (fig. S6; see Materials and Methods), mirroring past work that finds greater population variability when removing species with higher initial biomasses (see table S2) (25, 37). Higher levels of variability generate extinction cascades of non-harvested species (N-H extinctions; fig. S7), which threaten H when extinction cascades negatively affect its prey items. This causation chain of harvesting larger BH0 inducing higher population variability causing more N-H extinctions functions as the mechanism behind BH0, generating higher levels of N-H extinctions shown in Fig. 2C (also see figs S5B and S8). BH0 is also the best single pre-fishing predictor of N-H extinctions (table S3), and its effects are exacerbated when higher effort levels induce stronger perturbations in the community (Fig. 2C, fig. S6, and table S3). The majority of these N-H extinctions occur downstream from H (fig. S9A) and are trophically close to H (fig. S9B), with the average distance becoming closer with a higher number of trophic links to H (fig. S10). Consequently, we can show that the loss of H’s prey options increases proportionally with overall N-H extinctions (β = 1.5, P < 0.0001, R2 = 0.75; Fig. 2D). The higher amount of downstream extinctions relative to upstream extinctions is due to the generally higher trophic level of harvested fish (average of 3.36) compared to the overall species trophic level (average of 2.33). That is, there are proportionally more non-harvested species downstream than upstream of the harvested species. This is partially due to the size restrictions on our generated food webs (30 species maximum), although empirical food webs are also generally more species rich at lower trophic levels. Regardless, our results demonstrate that perturbations due to fishing can cycle through the food web and threaten even highly abundant harvested species through their prey items. The nonlinear effect of BH0 on its own survival indicates that one cannot necessarily assume linear relationships in ecological responses to perturbations.

In the open-access treatment, BH0 also played a critical role in fishery sustainability, as it was a principal driver of market dynamics. Open-access effort is dynamic, capable of both declines and growth (Fig. 3A). Growth in open-access effort is a function of BH0 and the maximum price of H (i.e., parameter a in Eq. 13). Higher max price reflects higher base demand for H, and higher BH0 provides potential yield to meet demand, thereby driving greater profits and effort (fig. S11A). Therefore, dependent on demand and yield, effort levels can range from low to high values. However, with sufficiently low demand or yield, net profit is consistently negative and effort declines to zero, meaning that the fishery fails to sustain itself (fig. S12). Consequently, harvesting the most abundant fish per food web sustained more fisheries (78% of simulations sustained) than harvesting randomly chosen fish populations (26% of simulations sustained). However, higher prices (i.e., demand) can sustain effort on rare species by supporting higher profits on lower yield (fig. S13). This potential growth in effort as a result of combined max price (a) and harvestable initial biomass (BH0) strongly predicts the peaks in effort early in open-access fisheries’ time series (Fig. 3, A and B, and tables S4 and S5).

Fig. 3 Open-access results.

(A) Example of effort dynamics where the red dot represents a local peak effort. (B) Peak efforts averaged across 2 years per product of starting biomasses of harvested species (BH0) and max price (a). Colors indicate the value of μ, economic sensitivity/reactivity of the effort per simulation. (C) Box plot showing BH0 of failed and sustained fisheries across different max price (a) when fishing the max BH0 fish population per food web. Increasing a supports sustained fishing effort in less abundant fisheries until spiking effort peaks produced by higher profits when fishing abundant species cause the H extinctions and fishery failure. (D) Three economic outcomes characterized by efforts averaged across the first fishing year, showing the importance of transient dynamics in the cycle of adjustment. Sustained (green) fisheries exhibit intermediate efforts, producing enough profit without collapsing the resource. Unsustained (orange) fisheries failed before causing species extinctions because their profits were lower than harvesting costs. Extinct (red) fisheries failed because high profits produce high effort peaks, causing the harvested species extinction. Significant differences are indicated by different letters (Tukey post hoc, P < 2 × 10−16). In boxplots of (C and D), boxes represent the interquartile range with the horizontal line showing the median, the lower box representing the 25th percentile, and the upper box showing the 75th percentile. Upper and lower lines extending from the boxes show the most extreme values within 1.5 times the 75th and 25th percentile, respectively. Outliers are shown as single dots.

However, growth itself, if unchecked, can also lead to effort declines during the “cycle of adjustment” (48), an empirically detected bio-economic process (49) reproduced in our simulation results. That is, effort declined (Fig. 3A) either when past fishing effort oversupplied H causing market saturation (reductions through price sensitivity to yield, via parameter b) or when past effort overfished H (reductions through a lack of biomass from which to profit). These mechanisms of reactionary effort reduction can function as “self-corrections” that potentially protect the harvested population (H) against excessive fishing effort, allowing their regrowth and recovery. These self-corrections, coupled with the potential for open-access fisheries to fail to sustain effort at low BH0, explain the relative lack of direct extinctions of the harvested species (H extinctions) in the open-access treatments (fig. S11B). However, combinations of sufficiently high max price (a) and BH0 can increase effort to levels (Fig. 3B) that cause H extinctions before self-correction occurs (Fig. 3C), especially in fisheries with high growth potential (high prices and BH0) and high sensitivity to profit (high μ in Eq. 11). After analyzing the range of this bio-economic feedback, we found three general categories in the outcomes of open-access fisheries that emerged from the early dynamics in the cycle of adjustment: (i) failed fisheries due to low BH0 and max price failing to maintain or grow effort, (ii) failed fisheries due to high BH0 and price driving H extinctions through excessive effort, and (iii) sustained fisheries existing in a middle ground between the two (Fig. 3D).

These three general economic outcomes underpin much of the relationship between the harvested species (H) and fishing effort in open-access fisheries. In particular, these results indicate that H extinctions were actually most common when harvesting highly abundant H (species with high BH0), as unregulated growth in effort depleted fish reserves (fig. S14). This also formed the basis of the nonlinear effect that BH0 had on open-access fishery persistence (Fig. 4A) and the level of effort sustained by the end of the simulation (Fig. 4B). In addition, while H extinctions were relatively less common in open-access fisheries, the association between harvestable biomass and effort growth reversed the positive relationship between BH0 and the time to H extinction seen in fixed-effort fisheries (compare Fig. 2A to Fig. 4C). Last, the dynamics of effort in response to yield also drive the different effects of fishing on non-harvested (N-H) species between open-access and fixed-effort fisheries.

Fig. 4 Open-access fishery’s economic and ecological sustainability results.

(A) Probability of producing a sustainable open-access fishery across the BH0 levels of H. The graph depicts the output of a binomial GAM regression with 95% confidence interval to demonstrate the nonlinear effect of BH0 on probability of fishery success. Full GAM regression model including economic factors log(a), μ, and b gives R2 = 0.66 with UBRE = −0.39. (B) Effect of BH0 of H on sustained open-access effort measured as average of final 400 time steps approaching t = 8000. The graph depicts output of a gamma GAM regression (log link) with 95% confidence interval to demonstrate the nonlinear effect of BH0 on final effort levels. Full GAM regression model including economic factors and square root–transformed effort levels results in R2 = 0.54 [GCV (Generalized Cross Validation) = 0.65] when using a Gaussian and a 61.5% deviance explained (GCV = 1.7) when using gamma (log link) distributions. (C) Scatterplot showing the effect of log(BH0) on time to H extinction in open-access simulations. Note the qualitatively reversed relationship from fixed-effort fisheries. The color gradient represents the max price effort (a) × the market reactivity (μ). Notice the reversal of the positive trend shown in Fig. 2A. (D) Bar graph with SE comparing the percent of the community lost through N-H extinctions in fixed-effort fisheries to open-access fisheries with comparable effort levels at t = 8000. (E) Hedge’s G effect size comparison of the two fishery treatments (fixed and open access) on N-H extinctions across different max prices (a). Dots represent effect size, and lines show 95% confidence intervals. Comparisons are made within webs at comparable effort levels (see Materials and Methods).

The profit-driven growth in effort caused by high BH0 means that (given sufficient demand) economic factors incentivize subjecting food webs to high levels of harvesting pressure on their more abundant species. This induces more variability in the biomass of the rest of the community than fixed-effort fisheries (fig. S15). Given the connection between this variability and non-harvested (N-H) extinctions (see fig. S7), market-generated variability expectedly induced more N-H extinctions than fixed-effort simulations. This was the case whether comparisons between fixed-effort and open-access results were made with attainable effort levels across all webs (Fig. 4D and figs. S16 and S17) or strictly within each web at the most comparable effort levels between fishery treatments (Fig. 4E).


Our study uses the ATN framework to study integrated feedback processes between ecological and economic dynamics by subjecting a wide array of food web topologies to two different harvesting treatments. Results demonstrate the notable nonlinear effects of the initial population biomass of harvested species (BH0) on both the ecological (e.g., Fig. 2B) and economic (e.g., Figs. 3C and 4B) sustainability of fishing. It also reveals the role of BH0 in driving variations in effort (Fig. 3B) and exacerbating the ecological costs of harvesting in open-access fisheries (Fig. 4, D and E). Temporary fluctuations in effort above sustainable levels are a known process in the cycles of adjustment seen in open-access fisheries (41), but we show here that the indirect effects of those fluctuations can drastically change food web structure through local species extinctions. Overall, our results indicate that considering the sustainability of an ecosystem service requires accounting for its surrounding ecology and its interface with the economy that it provides. In addition, the study, as a whole, demonstrates the flexibility of network approaches to integrate external anthropogenic variables into dynamic ecological processes.

Our study succeeds in expanding the purview of ecological networks and addressing the need for generalizability in fisheries studies, but further considerations are necessary in future developments. Important details in fishery output can result from spatial features of the fishery or from the ontogeny of the species in the aquatic community (5). Developing network ecology to the point of direct application in fisheries requires consideration of these ecological components. In addition, while single-species fisheries can have a large ecological and economic impact (43, 44), multispecies fisheries make up the majority of unique fisheries across the globe. Our model’s current form serves as a useful first step, but its flexibility makes it well poised in expanding toward multiple harvested species and other crucial developments.

Fisheries management has come a long way from focusing on single harvested species in isolation to considering the interconnected surrounding social, economic, and ecological conditions in which fishing occurs, commonly called ecosystem-based fisheries management (EBFM) (29). Such a holistic collection of considerations requires a multifaceted array of models and techniques to study (35). Our current model qualitatively reproduces past theoretical results (table S2) (37), well-documented empirical patterns (fig. S4), and output from models specifically trained on empirical data (fig. S18) (25). Open-access results were also qualitatively robust to different starting effort levels initiating the fishery (figs. S19 to S22). This gives us confidence in vetting our results and further developing this framework to contribute toward the creation of operational EBFM options for fisheries.


Experimental design

The objective of this study was to expand ecological network theory to include dynamic economic components. This Materials and Methods section describes the creation of ecological-economic networks and how they were analyzed. The ecological network models are composed of two major parts, the underlying structure of the network of interactions between species and the population dynamics of the interacting species (6, 39). In our model framework, we expand upon preliminary work (11, 41) and add economic dynamics to the established population dynamics.

Network/food web structure

Network structure describes the trophic links between resource, prey, and predator populations, irrespective of the strength of the link. Initially, 1100 trophic food webs were created as networks using the niche model (39). All food webs contained 30 trophic species with a connectance of 0.15 (37) within an error of 0.025 (39). Trophic species (S) define a population of individuals with similar resources and consumers. Connectance (C) is the fraction of realized trophic links (LS2), where L is the actual number of realized links. Trophic species are assigned a niche value across a single trait axis (considered body size in this case; see the “Ecological and economic dynamics” section below) defined by a feeding range and center along the trait axis. Species P eats species V if V’s niche value lies within P’s consumption range. This process does allow for the species’ consumption range to be at a higher niche value that allows for cannibalism. The consumption range of the lowest species in the niche axis is fixed to 0, assuming at least one primary producer/basal species. Each of the 30 trophic species is assigned values iteratively until (i) the web is connected (cannot be divided into two independent webs), (ii) every consumer species is linked to at least one basal species through a trophic chain, and (iii) the realized connectance is within the error range set by C. Trophic species that can be potentially harvested are labeled fish for ease of discussion. The fish species are chosen among the consumer species of each web with a Bernoulli’s trial (P = 0.6) and modeled as vertebrates, while the remaining species are modeled as invertebrates (11). Niche model food webs have been shown to exhibit empirically observed patterns in field webs (distribution of trophic species across different trophic levels, mean trophic chain length, etc.), especially in aquatic systems (39, 50, 51).

Ecological and economic dynamics

Population dynamics of the trophic species on the food web are regulated through the allometric trophic scaling (38), extended to a multispecies food web to function as an ATN (52). The ATN has advanced our understanding of ecosystems by providing flexible ways of expanding ecological models to incorporate dozens to hundreds of interacting species while maintaining tractable methods of analysis. The ATN models the biomass (μgCL) dynamics of trophic species through a system of consumer-resource (Eqs. 1 and 2)dBidtProducerbiomass=ri(1ΣkproducersBkK)Net primaryproductionΣjpredatorsxjyjejiFjiBjLoss bypredation(1)dBidtConsumerbiomass=xiBiMaintenanceloss+ΣjpreyxiyiFijBiGains bypredationΣjpredatorsxjyjejiFjiBjLoss bypredationqiEiBiLoss byharvesting(2)

Parameter values shown in Eqs. 1 and 2 are described in Table 1 with differences based on life history traits labeled in consecutive entries. The functional response of consumers to prey biomass is given in Eq. 3. Note that B0, j denotes the half saturation biomass of species j and that ωij denotes the preference of consumer i for prey species j. Each preference, ωij, equals the inverse of the total number of i’s prey species and changes through time when prey go extinctFij=ωijBjhB0,jh+Σkprey speciesωikBkh(3)

Table 1 Parameter/function definitions, values, and sources.

N/A, not applicable.

View this table:

The applicability of the ATN framework to models with such a large number of species stems from the ability to create plausible parameters for the physiological rates through a negative-quarter power law relationship with species’ body masses, the single trait axis described in the “Network/food web structure” section (38, 53). Part of the analytical power of this approach is that ecological interactions are reasonably parameterized, thereby allowing researchers to focus analysis and sensitivity analysis on other model aspects, in our case, food web structure and fishery parameters. This is realized through the three rates: reproduction, R; metabolism, X; and maximum consumption, YRp=arMp0.25(4)XC=axMC0.25(5)YC=ayMC0.25(6)where ar, ax, and ay are allometric constants that determine rates based on the body size Mi. The subscripts C and P indicate consumer and producer parameters, respectively (38).

The time scale of the dynamics in each food web is defined by setting the mass-specific growth rate of the basal species (the producers) to one. With this as a reference and the assumption that basal species share the same body size, the mass-specific metabolic rates of all species are normalized by the time scale (RP), and the maximum consumption rates are normalized by the metabolic ratesri=1(7)xi=XCRp=axar(MCMP)0.25(8)yi=YCXC=ayax(9)

The body masses of predators are expressed relative to basal species, and the body size ratio between predators and prey is considered to be a constant (Z), a reasonable approximation in an aquatic system (52). This allows xi to be expressed as belowxi=axar(ZTi1)0.25(10)where Ti is the prey averaged trophic level of species i calculated from network topology (50). Model time is set similar to past work (54) where one model time step equals one real-time day.

Loss of a population’s biomass due to harvesting (labeled “loss by harvesting” in Eq. 2) is measured as the rate of fishing effort, Ei (Eq. 11), on species i (55). Effort is a broadly applicable index used to measure the amount of fishing/harvesting taking place in a fishery, including capital and labor (41). Depending on the specific fishery, effort can track the number of fishing lines, boats, workers, work hours dedicated to harvesting, etc. The results presented in this work would not qualitatively change on the basis of the specific details of the effort metric. This study focuses on single-species fisheries. Therefore, Ei = 0 for all species that are not the harvested fish species, H. The fishing effort on H is greater than or equal to 0 and changes in E derive from the product of net profit and an adjustment parameter, μ, representing the economic sensitivity of fishing effort to changes in net profit (55). The economic sensitivity of the fishery’s effort (μ ∈ [0,1]) describes the sensitivity of the effort to changes in profit or loss. Net profit is defined as the product of price per unit of biomass harvested and the actual yield (Y) caught at the current E level (Eq. 12) minus the costs per unit effort (55). The yield at any given time in the model is a product of the current level of effort, the available biomass of the harvested species (H), and q, the “catchability” of H per unit effort. Yield translates into supply and informs the market price, p, through the linear Eq. 13, where a represents the maximum price and price sensitivity to yield is labeled b. Equation 13 is incorporated as a piecewise equation such that p = 0 when Y1b. Last, costs are removed from gross profit to reach net profit by subtracting coE from pY, where co is the fixed value of cost per unit effort. All open-access fisheries start with E(0) = 1 to model a fishery from its initiationdEidtEffortlevel=μSensitivityto net profit*pYGrossprofitcoEiCostsper effortNet profit(11)Y=qEBH(12)p=a(1bY)(13)

Unlike optimization studies that can analytically seek out effort levels that maximize profits, the open-access equations (1113) move toward equilibrium (dEdt=0), as total revenue approaches total costs and economic profits reach zero.

Experimental setup and treatment design

All 1100 food webs created with the niche model were randomly assigned initial biomass conditions per species using a uniform distribution, U ∈ (0,1). The ATN framework shown in Eqs. 1 and 2 was then used to simulate each food web for 4000 time steps in a fishing-free stage (effort = 0). This initial fishing-free period limits possible effects of transient dynamics on results of fishery treatments. After the initial fishing-free period, only conserved webs that met fishery criteria are chosen to be subjected to the two different fishery treatments. A food web is considered conserved if (i) it is connected, (ii) every consumer species is linked to at least one basal species through a trophic chain, (iii) the number of remaining trophic species is higher than or equal to 20, and (iv) it has at least one fish species.

These criteria resulted in 480 unique conserved food webs that were subjected to an additional 4000 time steps of simulation in the fishing stage. Subjecting the initial food webs to these criteria eliminated the potential for analytical comparisons of fishery effects between categorically different types of food webs. For example, criterion (i) not being met in a food web would potentially protect certain species from any perturbations caused by fishing, making it an incompatible comparison to connected food webs perturbed by fishing.

Fishing effort was applied as one of two fishery treatments, fixed effort and open access. In the fixed-effort treatment, dEdt=0 and fishing effort is held constant from t = 4000 to 8000 during the entire fishing treatment. The tested fixed-effort values are E = [0,1,2,3,4,5,6,7,8,9,10,12,15,20]. In the open-access treatment, effort varies on the basis of the interactions between economics and ecology through profit on yield described above. Open-access results are tested across the values of effort sensitivity to profit μ = [0.05,0.2,0.5,0.8,1.0], price sensitivity to yield b = [0.01,0.05,0.1,0.5], and maximum price a = [10,20,40,70,100,150]. Initial effort levels in all open-access simulations are set at 0.5, 1, or 2 to model a fishery from its beginning. Initial efforts of 1 are presented in the main text with initial efforts of 0.5 and 2 presented in the Supplementary Materials (figs. S19 to S22). Any larger initial efforts would take simulations away from modeling fisheries from their inception. Roughly 99.5% of our open-access simulations result in final effort levels ≤20. Therefore, exploring the range of fixed-effort values between 1 and 20 is the exact range that we want to cover so that we can compare management strategies at different effort levels.

Regardless of fishery treatment, fishing effort was applied to a single fish species per simulation, H, chosen at the end of the initial stage in conserved webs, at t = 4000. We do not model bycatch or multiple-species fisheries. In the fixed-effort treatments, every single unique fish species was harvested in every web. In the open-access treatment, given the importance of population biomass seen in the fixed-effort results, the identity of the harvested fish species was selected either (i) randomly across all available fish in the food web, labeled “Random” in figures and results, or (ii) as the fish species with the highest biomass after the initial stage (BH0 at t = 4000), labeled “Max” in figures and results. Webs with only one harvestable fish species were only used once per parameter combination per fishery treatment. These open-access results could then be directly compared to their fixed-effort counterparts. All simulations were completed using MATLAB 2018a and the solver ode45 for numerical integration (relative and absolute error tolerances both equal to 10−8).

Data collection and categorization

Ecological, economic, and network data were compiled initially at t = 0, at the end of the initial fishing-free stage at t = 4000, and at the end of the fishing stage t = 8000. Species extinctions during simulations were considered at the threshold, Bext, listed in Table 1 and are displayed in figs. S3 and S11B. Beyond the time series data of each species in each food web simulation, a large number of variables describing pre- and post-fishing overall network attributes (e.g., connectance) and local network structure around the harvested species (e.g., number of direct trophic links) were also considered. A full listing of these variables appears in table S6.

Per simulation, 408 attributes were considered in both the fixed-effort treatment and open-access treatment. Additional variables were added to the open-access treatment, for a total of 432, to detail the dynamics of effort through the simulation and to make direct comparisons to the fixed-effort treatments (see Fig. 4E). While effort levels were directly recorded in open-access treatments, simulation results were also categorized as sustained or failed at the end of simulations. To avoid considering temporary troughs in effort at the end of simulations as complete failures, effort levels were averaged over the last 400 time steps. Fisheries that did not maintain effort levels above 1% of starting effort values were labeled failures, while those that did were labeled sustained.

In each food web simulation, dynamic variability of each population was measured using a population’s biomass time series’ coefficient of variance (C.V.). A similar process was used to measure the variability in effort in open-access fisheries. For the community’s dynamic variability, we used the mean of each population’s biomass time series’ C.V. as a proxy. In other words, we took the C.V. of each species’ biomass time series and averaged them to get the mean variability of a community. This was done in the 400 time steps before the start of fishing, in the first year after the start of fishing, across the first 3 years after the start of fishing, and during the last year of fishing before each simulation ended. See table S6 for the full list of analyzed factors.

Statistical analysis

Given the large number of variables/factors, we used classification and regression tree analysis (56) through JMP Pro 13 to obtain variables of interest in the search for drivers of fishing-induced changes across our webs and simulations. Once obtained, potential drivers were further explored using R 3.5.1 statistical software. Continuous variables were initially explored using generalized linear models. Binomial regressions were used for binary response data, while gamma or Gaussian distributions were used for continuous response data. Models were trained and vetted on various resampled subsets of data to assess consistency of the model results and avoid the pull of outliers. In the case of biomass of the harvested species (BH0), resampling indicated potential nonlinear responses. As a result, we used generalized additive models (GAMs) to account for and visualize BH0’s effect on response variables using the mgcv package in R (57). When response variables were continuous and limited between 0 and 1 (fig. S7), beta regressions were used (beta package R) (58). Categorical variables can be incorporated into GAMs as factors (table S3). In situations where GAMs were not appropriate, comparisons across categorical variables were done using Tukey’s post hoc test (Fig. 4D).

There are two major comparisons of ecological cost of harvesting through the non-harvested (N-H) extinction prevalence between the two fishery treatments. The first (Fig. 4D and figs. S16 and S17) compares the average N-H extinctions associated with comparable effort levels between the fixed-effort and open-access fishery treatments. Because of inherent variability in effort levels in open-access treatments, not every fish in every web can be guaranteed to be fished at every effort level, as was done in the fixed-effort treatment. Therefore, open-access effort per simulation was averaged across three time periods—the first year (fig. S16A), the third year (fig. S16B), and across the last year (Fig. 4D)—and then grouped into one unit effort block. For example, 0 < Effort < 1, 1 < Effort < 2, 2 < Effort < 3, and so on along the tested fixed-effort values. N-H extinctions from those open-access simulations were then compared against these from the fixed-effort simulation.

The first method’s focus on comparisons across effort requires comparisons across different webs. To constrain the N-H extinction assessment to within web comparisons, a second method was implemented. For each max price value (a), every single open-access simulation’s N-H extinction level was compared to the N-H extinctions from the fixed-effort treatment of the same fish in the same web with the closest matching fixed-effort level when compared to the open-access effort level averaged across the third year of open-access fishing. The third year was chosen because it captured the important transient behavior (when most changes to the food web occur) and was generally the beginning of simulations reaching asymptotic behavior. The effect of fishery treatment was ascertained using Hedge’s G (effsize R package) (59) to compare the effect size of fishing with open access with the fixed-effort treatment as the control (Fig. 4E). The same process was also used to compare other metrics, such as mean community biomass variability (fig. S15).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to thank K. Kroetz, M. I. Avila, the Valdovinos Lab, and the researchers at the “Socio-Environmental Networks of Common Pool Resources” Working Group at the National Socio-Environmental Synthesis Center (SESYNC) for feedback on study design, results, and interdisciplinary communication. We would also like to thank J. Megahan for assistance in producing the figures for publication. Funding: No outside funding was applied to this study. Author contributions: F.S.V. conceived the study. F.S.V., P.G., and V.C. designed the study. V.C. and P.G. wrote and implemented the simulation code. P.G. and F.S.V. completed the analysis and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Simulation code and data is available at the repository Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article