StrathE2E2: An r package for modelling the dynamics of marine food webs and fisheries

We present StrathE2E2, an r package for modelling the whole ecosystem, or ‘big‐picture’ effects of hydrodynamics, temperature, nutrient additions and fishing on continental shelf marine food webs. StrathE2E2 has two linked parts—a fishing fleet model and an ecology model. The fishing model combines harvesting, discarding and seabed disturbance rates across a range of gears and passes the results into the ecology model. The ecology model is a network of coupled ordinary differential equations representing the rates of change in nitrogen mass of organic detritus, dissolved inorganic nutrient and coarse guilds of living biomass spanning microbes to megafauna. The equations include representations of feeding, metabolism, reproduction, active migrations, advection and mixing. Environmental driving data include temperature, irradiance, hydrodynamics and nutrient inputs from rivers, atmosphere and ocean boundaries. The package includes functions for parameter optimization, global sensitivity analysis and Monte Carlo estimation of credible intervals for model outputs. A fully developed and documented implementation for the North Sea is included with the package.


| INTRODUC TI ON
The effects of anthropogenic or natural pressures applied to any part of an ecosystem are eventually felt everywhere to some extent through the phenomenon known as a 'trophic cascade' (Pace et al., 1999). Cascading effects are attenuated or amplified as they propagate through the food web, depending on the nature of the pressure and details of the ecology (Heath et al., 2014). Diagnosing the type and magnitude of pressures that an ecosystem can sustain before being fundamentally altered requires simulation with mathematical models that aim to represent the key ecological components and processes which govern cascades.
We present the package StrathE2E2 for the r statistical environment (R Development Core Team, 2014), which models both bottom-up and top-down trophic cascades in shelf-sea ecosystems, spanning inorganic and organic nutrients through to birds and mammals. The model takes a macroscopic view of ecology, aggregating over the many microscopic details of taxonomy, demography and spatial structure (Giricheva, 2015). The aim is to represent the gross dynamics with a tolerable parameter count and fast run-time, so as to enable 'big-picture' strategic scenario analyses. The basic model is supported by functions for computational parameter optimization, sensitivity analysis, estimation of credible intervals of model outputs and network analysis.

| Ecology model general description
The ecology model, developed from an earlier prototype (Heath, 2012), is a network of mass conserving coupled ordinary differential equations (ODEs) describing spatially averaged rates of change in state variables representing organic detritus, dissolved inorganic nutrient and living biomass. To simplify the description we can think of the variables as being divided between two coupled sub-networks: a predator-prey network-the food web-and a nutrient recycling network. Between the two, all marine lifeforms are explicitly or implicitly accounted for, but aggregated into coarse groups or 'guilds' defined mainly by feeding characteristics and diet preferences ( Figure 1). All state variables, except macrophytes, are expressed solely in terms of nitrogen mass, since this element is the most commonly limiting in temperate shelf seas.
Macrophytes are expressed in terms of both nitrogen and carbon mass with dynamic stoichiometry since these organisms have an exceptional capacity to seasonally absorb and store nitrogen (Appendix S1).
Each ODE comprises a set of rate-of-change terms representing a variety of biological and physical processes (Appendix S1).
Biological terms describe the balance between gains due to assimilation of food, and losses due to mortality and metabolism.
Some components of the food web (planktivorous and demersal fish; suspension/deposit feeding and carnivore/scavenge feeding benthos) are resolved into life stages, and for these the equations F I G U R E 1 Schematic of the food web compartments of the StrathE2E2 model. Green arrows represent advection, mixing and migration; orange arrows represent fishery-related fluxes; black arrows represent biological fluxes. Red labelled components are active migrators while blue are subject to passive advection and mixing and black are anchored. Pale blue boxes represent quantities that are exported from the model while yellow are imported. The model also includes fluxes from living components to ammonia, detritus and corpses due to excretion, defecation and death but these are not shown for clarity. Also for clarity, birds, pinnipeds and cetaceans are combined as a single box but in the model are separate entities.

Cetacean landings
also include the balance between gains due to recruitment and losses due to developmental progression or spawning. In addition, all components of the model are replicated across homogeneous spatial compartments, so each ODE also includes terms representing sinking, advection, mixing and migration flows through the system.
The model has a coarse spatial structure consistent with the coarse guild definitions of the living and chemical components of the system. The spatial domain is first divided horizontally into two bathymetric/hydrographic zones ( Figure 2) which represent the most basic spatial division in shelf seas-a shallow, vertically mixed zone mostly influenced by tides and freshwater inputs, and a deeper, potentially seasonally stratified zone mostly influenced by exchange with an external ocean (Pingree & Griffiths, 1977).
For convenience we refer to these as the inshore and offshore zones respectively, though there is no necessity for the inshore zone to be adjacent to the coast-in principle it could represent a shallow offshore bank. The water column in the offshore zone is divided vertically into two compartments or layers, while the inshore zone is represented by a single compartment. The inshore and offshore zones are each divided into discrete seabed habitat types comprising exposed rock and up to three compartments of different sediment properties per zone. The sediment habitats are notionally mud, sand and gravel, but defined by median grain size and natural disturbance rates. It is not necessary or computationally efficient to represent each state variable in every spatial compartment; for example, cetaceans do not need to be resolved vertically. State variables are therefore resolved hierarchically to spatial compartments with the largest (in terms of body size) and/ or most mobile guilds being represented at the coarsest spatial resolution (Appendix S1).

| Predator-prey connections, demography and mortality
Ingestion of prey by a predator is governed by a preference matrix and a standard Type II functional response in which per-unit-biomass predator consumption rates increase asymptotically towards a Q 10 temperature-dependent maximum with increasing prey concentration (Appendix S1). A proportion of ingested food becomes new body mass in the predator. The remainder is divided equally between fluxes to organic detritus and ammonia, to represent defecation of undigested material and food-dependent metabolism. Background (non-feeding) metabolism increases with temperature but with a higher Q 10 than maximum uptake rates, so the net result is that productivity, that is, production rate per unit biomass, will exhibit a dome-shaped response to temperature.

| Nutrient recycling network
Six forms of organic detritus are represented in the recycling network: suspended material, labile and refractory sediment material, 'macrophyte debris', 'corpses' and 'discards'. Both the suspended and sediment fractions implicitly include dissolved and particulate organic matter and associated bacterial flora, and are formed in the living food web by defecation and density-dependent mortality fluxes from plankton and the larval stages of fish and benthos.
Corpses are produced by density-dependent mortality of fish, benthos and top predators, and the decay of discards. The latter are a short-lived, special form of detritus generated as a byproduct of fishery harvesting. Macrophyte debris is created by wave and density-dependent destruction of living macrophyte forest biomass. All forms of detritus are regarded as a potential food source for detritivorous and scavenge feeding guilds of living organisms.
The dynamics of each detritus and dissolved nutrient category are governed by an ODE in which the rate-of-change terms correspond to the production and consumption rates elsewhere in the food web, plus physical flows between spatial compartments. Q 10 temperature-dependent coefficients govern transformations between different forms of detritus, conversions of detritus into ammonia (mineralization); ammonia to nitrate (nitrification), and nitrate to nitrogen gas (denitrification). To complete the biogeochemical cycle, nitrate and ammonia are re-absorbed into the food web by phytoplankton and macrophytes, governed by light and temperature-dependent uptake responses.

F I G U R E 2
Schematic showing the horizontal and vertical spatial structure of the model. The compartments S0-S3 and D0-D3 refer to inshore/shallow and offshore/deep seabed habitats respectively. S0 and D0 are rock habitats which reflect, rather than absorb, settling material back into the water column. S1 and D1 represent fine (muddy), S2 and D2 medium (sandy) and S3 and D3 coarse (gravelly) sediments respectively

| External boundary fluxes
Passive hydrodynamic influxes of dissolved nutrient, suspended detritus and phytoplankton into the model domain from adjacent sea areas are defined by driving datasets (Appendix S1). These comprise water volume influxes and associated boundary concentrations of advected material, the products of which (volume flux × concentration) represent mass fluxes. These are incorporated into the rate terms in the relevant ODEs describing changes in the receiving spatial compartments (Appendix S1). Additional inflows due to dry and wet deposition of atmospheric nutrient to the sea surface, nutrient inputs from river discharges, and other unspecified sources (e.g. aquaculture) are accounted for in the same way. proportional loss rate from the guild to an ocean sink.

| Interior fluxes between spatial compartments
The vertical and horizontal links between sediment and water column compartments within the model domain are represented by combinations of passive advection, diffusion and sinking and active migration fluxes (Appendix S1). Dissolved nutrients, detritus and phytoplankton are subject to passive vertical and horizontal exchanges. Zooplankton and larvae of fish and benthos are subject to passive horizontal exchanges, but active vertical. Fish and top predators (pinnipeds, cetaceans and birds) undertake active horizontal migrations, but their vertical distributions are not resolved.
Passive exchanges are the product of dynamic differences in concentrations between vertical or horizontal spatial compartments, scaled by hydrodynamic mixing coefficients supplied as time-varying parameters (Appendix S1). Mixing coefficients between sediment habitats and the water column are defined by sediment permeability, modified by representations of bioturbation by deposit feeding benthos, natural erosion by bed shear stress and fishing-related abrasion.
In reality, plankton and larvae undertake active diel vertical migrations but we do not resolve these sub-daily behaviours. Rather,

| Representation of fishing in the ecology model
Living biomass guilds considered vulnerable to targeted capture or incidental by-catch by fishing gears are the top-predators (birds, pinnipeds and cetaceans); planktivorous, demersal and migratory fish; carnivorous/scavenge and suspension/deposit feeding benthos, carnivorous zooplankton and macrophytes. In each case, the fishing process is represented in the ODE for each guild by a 'harvest ratio' (proportion of instantaneous biomass captured per unit time).
A proportion of the catch is directed to the discards class, comprising whole-animal rejects and viscera arising from at-sea processing of the remaining catch. Only the residual fraction of the catch weight is exported from the model as landings.
The collateral effects of fishing activity-release of sediment pore-water nutrients, resuspension of sediment detritus and damage mortality of benthos-are driven by the area-proportion of each seabed sediment habitat abraded per unit time by fishing gears.

| Fishing fleet model description
The fishing fleet model is a static, matrix-based scheme which generates harvest ratios and discarding rates for each guild in the ecology model, and abrasion rates for each seabed habitat, due to the combined actions of up to 12 different fishing gears. (Appendix S2).
Key inputs are, for each gear type, the spatial distribution of activity density, catching power, selectivity, discards and at-sea processing rates for each ecology model guild, and contact rate with the seabed (Appendix S2). Activity density is defined as the deployment duration of a given gear per unit sea surface area in a given time interval, integrated across all vessels (units: m −2 ). The power of a gear is a measure of its efficiency at catching biomass of a given resource guild. The product of activity density and power is a quantity that we refer to as fishing effort. For a given resource guild, effort is proportional to the harvest ratio and so can be summed across gears. The core network of differential equations is coded in C for solving with the lsoda function as implemented in the 'deSolve' package for r (Soetaert et al., 2010). lsoda swiches automatically between stiff and non-stiff methods-stiff methods are required when high rates of change are encountered to avoid numerical instabilities. All The package documentation includes guidance on how R users can extract data from the returned structures so as to conduct their own analyses. In addition, the.csv outputs offer language-neutral analysis options.

| E X AMPLE MODEL OF THE NORTH SE A
The geographic domain for the North Sea example model provided with the package is divided between inshore and offshore zones at around the 30-m isobath (Appendix S4). The model was optimized for two contrasting periods of fishing, temperature, hydrodynamics and river nutrient emissions (1970-1999, cool temperatures and high harvest ratios for finfish; 2003-2013 warmer temperatures and low finfish harvest ratios), with a common set of biological parameters for the ecology model ( Figure 3). Full details of the physical configuration, assembly of driving data and the observational data for fitting are provided in Appendix S4.
As an example of the many possible ways of using the model to explore scenario situations, Figure 4 shows the difference between a baseline model representing the stationary state with 1970-1999 environment and fishing activity, and the stationary state for a scenario in which all towed fishing gears are prohibited in the inshore zone and displaced offshore. Note that the overall activity rates of the gears remain unchanged. This scenario directly affects the fish and benthos guilds which are the main targets of the mobile fishing gears. However, the ecosystem response also extends to cetaceans, pinnipeds, birds, zooplankton, phytoplankton and nutrients through indirect 'who-eats-whom' network effects which cascade through the food web. Other illustrative scenario experiments are presented in Appendix S4.

| APPLYING THE MODEL TO NEW REGIONS
The ease of applying StrathE2E2 to a new region will depend on the availability of driving data, and observational data to form  Development, assembly and refining of these data for the North Sea model was a substantial task, but copying and adapting this setup to represent a different region is well within the scope of, for example, a Masters or PhD project. Our experience is that a credible draft model can be produced in a few days by roughly adapting the North Sea example given basic knowledge of the environmental driving conditions such as annual cycles of light and temperature, bathymetry, seabed sediment properties and fishing fleet compositions.

| D ISCUSS I ON AND FUTURE DIREC TIONS
There is no simple answer to the question of an appropriate spatial, taxonomic and biological granularity for an ecosystem/food web model (Fulton, 2010;Giricheva, 2015;Iwasa et al., 1987). In StrathE2E2, our aim was to provide a 'big picture' marine ecosystem modelling tool amenable to computational parameter optimization, sensitivity analysis and exploration of uncertainty in model outputs, on ordinary desktop computing hardware. Achieving this required us to sacrifice some spatial, taxonomic and biological granularity in order to span the ecosystem and food web from physics, nutrients and microbes through to megafauna and fishing fleets. Other end-to-end models which aim to Partial models, which dynamically represent only a subset of taxa of interest at greater spatial, taxonomic or demographic resolution (e.g. Ecopath with Ecosim (which does not dynamically represent the physics, nutrient or microbial system, Christensen & Walters, 2004);Speirs et al., 2010;Scott et al., 2014) can be problematic for simulating cascading effects because propagation is sensitive to the way in which boundary conditions are represented, especially for short food chains (Heath et al., 2014).

Inclusion of tools for exploring uncertainty in model outputs
and parameter sensitivity is a relatively new development for ecological modelling packages (Steenbeek et al., 2018;Wu et al., 2013).
However, interfaces to these techniques are among the essential requirements for future generations of models, allowing uncertainty in the outputs to be attributed to different sources of uncertainty in the inputs (Pianosi et al., 2016). This is increasingly required for corroboration, quality assurance and the defensibility of model-based scenario analyses. In the StrathE2E2 package, we have provided integrated analysis functions based on some established methods (Simulated Annealing for optimization, Morris Method for sensitivity analysis), but there are many other methods available (e.g. Sobol Methods; Sobol, 2001), and Bayesian methods using Markov chain Monte Carlo algorithms). Users can explore these by simply embedding our e2e_run() function (which executes a single model run with a given parameter set) into their own scripts to encode different methodologies.

ACK N OWLED G EM ENTS
The main financial support for the development of StrathE2E2

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.

DATA AVA I L A B I L I T Y S TAT E M E N T
The package, complete with documentation, the North Sea model and all data shown in this paper, is available to download from CRAN (https://CRAN.R-proje ct.org/packa ge=Strat hE2E2) or the GitLab repository https://marin ereso urcem odell ing.gitlab.io/index.html.