Addressing initialisation uncertainty for end-to-end ecosystem models: application to the Chatham Rise Atlantis model

Introduction

Atlantis models are end-to-end ecosystem models that can be used to model everything in an ecosystem from sunlight to market, including feed-back loops and flow-on effects (Audzijonyte et al., 2017a). There are many published Atlantis models (Link, Fulton & Gamble, 2010; Ainsworth et al., 2011; Ainsworth, Schirripa & Morzaria-Luna, 2015; Fulton et al., 2011; Weijerman et al., 2014; Sturludottir et al., 2018; Ortega-Cisneros, Cochrane & Fulton, 2017; Porobic et al., 2019; McGregor et al., 2019b), and sensitivities and uncertainties are beginning to feature in the more recently published work (Ortega-Cisneros, Cochrane & Fulton, 2017; Sturludottir et al., 2018, Ortega-Cisneros et al., 2018; Hansen et al., 2019a, 2019b; McGregor et al., 2019b).

An important step that has not yet been addressed for Atlantis models, or any end-to-end ecosystem models, is to address initialisation uncertainty. In a review paper, Payne et al. (2015) found generally marine ecosystem models have not explicitly addressed uncertainty of initialisation, and more recently, Hansen et al. (2019a) noted it is not something that has been done for Atlantis models. Payne et al. (2015) speculated as to the likely effects of initialisation uncertainty in end-to-end models such as Atlantis, noting long-lived species might dampen effects, and short-lived species may amplify them. While accounting for uncertainties of such complex models is no small task, it is possible to focus on a sub-component of model dynamics rather than the entire model. Atlantis models are particularly well suited for this approach as they are structured using components. These components dynamically interact, and they define functionality for the physical, biological, fisheries, harvest and economic components of the system (Audzijonyte et al., 2017b). Rather than vary the initial conditions for thousands of spatially explicit parameters for the entire ecosystem, we focus this study on the biological sub-component and explore sensitivity to the initial conditions of species groups. These will have flow-on effects to the whole system, but the analyses focus on the effects as they relate to population dynamics of the modelled species.

Ecological theory has featured stability, chaos, the importance of initial conditions, and how these relate, although definitions seem to vary. The most commonly known attribute of chaotic dynamics is sensitivity to initial conditions, often referred to as the butterfly effect. In a system with chaotic dynamics, a very small change in its initial state will result in a different trajectory (Baker, Baker & Gollub, 1996). Stability often relates to variability, persistence and robustness of a system (Tilman, 1996; Ives, Gross & Klug, 1999). Multiple studies have investigated characteristics of an ecosystem that are linked with stability, both from a theoretical perspective, and from observation. May (1972) showed mathematically, that large complex systems with high levels of diversity are unstable. However, there seem to be exceptions to this rule, as later studies have shown. Roberts (1974) argued that most systems in practice appear to be more stable with more connections—contrary to the mathematical analysis of May (1972). Roberts (1974) showed if only feasible solutions are included in the analyses, such that no species may have a negative population, larger systems are actually more stable. May (1975) examined non-linear difference equations with respect to chaotic, cyclic and stable biological dynamics. Other aspects subsequently shown to increase stability of ecosystems include negative pairwise correlations (Tang & Allesina, 2014), species dispersal (Allesina & Tang, 2012), modularity (subsets of closely connected components) (Grilli, Rogers & Allesina, 2016), predator-prey relationships (Tregonning & Roberts, 1979), a high proportion of weak interactions (Olsen et al., 2016) and spatial structure (Fulton, 2001).

The Chatham Rise Atlantis model was selected for this study as it has been shown to have sound dynamics with respect to our current knowledge of the system (McGregor et al., 2019b), has been tested for sensitivity to its ocean forcing variables (McGregor et al., 2019b) and spawning stock recruitment assumptions (McGregor, Fulton & Dunn, 2019a), but is yet to be tested with respect to its initial conditions. The initial conditions for the Chatham Rise Atlantis model were specified to reflect the ecosystem in its unfished, or virgin state. We have varying levels of understanding of the components of this ecosystem in its unfished state, and as such, there are varying levels of confidence around the estimates for the initial conditions, with all components having some level of error.

This study presents an approach for addressing and exploring sensitivity to initial conditions of the biological sub-component of an end-to-end ecosystem model, with application to the Chatham Rise Atlantis model. We go beyond the question of whether the model is sensitive to its initial conditions, and analyse what features of the Chatham Rise ecosystem, and how we have modelled it, affect this answer. We discuss the likely impacts of our findings for future use of this ecosystem model. We highlight areas of potential future research with respect to model development, and to support decisions relating to the sustainable use of the Chatham Rise marine ecosystem resources.

Methods

The Chatham Rise Atlantis model is spatially defined using 24 dynamic polygons, and five water column depth layers. Species are modelled using 55 species groups, which include species of bacteria, detritus, phyto-plankton, invertebrates, fish, sharks, cetaceans and birds. Some species groups were modelled as biomass pools, and some with age-structure, using numbers-at-age and mean weight-at-age. For many of the species, we have estimates of biomass, growth rates, age of maturity, natural mortality, spatial distributions and diets, although some species have more knowledge gaps than others. McGregor et al. (2019b) characterised the species groups by keystoneness, responsiveness, and informance, and brought these together to highlight which data gaps are likely to influence model results. We can use these attributes to perturb the initial conditions in a meaningful way based on likely uncertainties, and they may add context when analysing the results. Keystoneness measures the effect changes in biomass of a species group has on the rest of the system; responsiveness measures how responsive a species group is to changes in biomass of other species groups within the system. Informance was a qualitative measure used to reflect both how well informed each species group was, and how well it performed in the model. It considered whether key dynamics such as growth, mortality rates, diets, and responses to fishing were all realistic based on current knowledge of these dynamics.

Sensitivity to initial conditions is characterised by two criteria in this study: (1) persisting variability between simulations, assessed at the species group level; (2) consistency of responses to fishing pressure, both at the species group level and at the system level. Variability of species biomass between simulations with perturbed initial conditions is important to consider as it may affect conclusions drawn from the model, and can provide error bounds around results. The consistency of responses to fishing pressure has also been considered as it is often the direction of a response that is of interest when exploring scenarios using ecosystem models. For example a conclusion that increasing fishing by 20% reduced diversity by 5–15% (with the range presented perhaps reflecting the initialisation uncertainty) is more informative than whether simulated biomasses of individual species groups varied across runs by less than 20%.

The analyses presented here were carried out in three main sections: (1) components of the modelled system were characterised with respect to attributes that may affect sensitivity to initial conditions; (2) the initial conditions were perturbed, and the resulting model simulations were compared; (3) correlations between component attributes and responses to perturbations of the initial conditions were analysed, thus linking the first two sections.

The models were run on a Linux cluster using Moab job scheduler, and used Atlantis version 6262M (trunk). We created a new initialisation.nc file for each set of perturbed initial conditions, and these were created using R (version 3.6.1) to read in the base initialisation file, scale the values, then create the new initialisation file. R was also used to create run files that interacted with Moab to schedule the model runs. Scripts used for this work are available on GitHub (McGregor, 2018) (https://github.com/mcgregorv/CRAM_chaos). The R scripts were all developed in a Windows environment.

Varying initial conditions

We varied the initial conditions for the number-at-age variables of age-structured species groups, and the biomass of biomass-pool species groups (Tables 1 and 2). For the biomass-pool species groups, biomass is the only option to perturb. Age-structured species groups could have errors in the specification of numbers-at-age and/or size-at-age, both of which affect the biomass-at-age. In a stock assessment model, size-at-age (or growth rates) are generally the same with respect to time, whereas there is more likely a difference in numbers with respect to time (especially before fishing compared to after fishing). Hence, a different virgin biomass in a stock assessment model would generally be made up of a different number of fish, rather than the same number of fish but a different size. To align with this, we perturbed numbers rather than size for the age-structured species groups. The resulting number of variables to perturb was 361, of which 341 were numbers-at-age (number of age-classes for a species group ranged from 2 to 10, with 10 being most common), 18 were the nitrogen content of biomass-pool species groups, and two were the silicate content of biomass-pool species (diatoms and microphytobenthos). The numbers-at-age of all age-classes for a given species group were scaled by the same amount for each simulation, such that the proportions-at-age were preserved. Again, this aligns most closely with the structure of a population’s virgin state in a stock assessment. The assumed M (instantaneous natural mortality) is preserved through the proportions-at-age based on an exponential decay curve, and the total numbers are what change when varying the initial conditions. Age-structured species groups were modelled with between 2 and 10 age-classes, and were perturbed by applying one scalar for all age-classes of a given species group. This significantly reduced the number of scalars required to 57, of which 37 were for scaling numbers-at-age for age-structured species groups.

Table 1:

List of age-structured species groups for CRAM (McGregor et al., 2019b).

‘Keystone’ is the keystone ranking from McGregor et al. (2019b), where 1 is the highest. ‘Informance’ is the informance rating from McGregor et al. (2019b) where 1 is the highest (‘No data gaps, performed well, abundance index available’) through to 4 which is the lowest (‘Poorly specified’). Ben: benthic; Dem: demersal; invert: in-vertivore; pisc: piscivore; Invert comm: commercial invertebrates; herb: herbivore; scav: scavenger.

Specices group	Main species	Keystone	Informance
Arrow squid	Arrow squid	19	3
Baleen whales	Southern right whales (Eubalaena australis)	25	4
Basketwork eel	Basketwork eels (Diastobranchus capensis)	18	3
Baxters dogfish	Baxter’s dogfish (Etmopterus baxteri)	28	2
Ben fish deep	Four-rayed rattail (Coryphaenoides subserrulatus)	35	2
Ben fish shal	Oblique banded rattail (Coelorinchus aspercephalus)	5	1
Black oreo	Black oreo (Allocyttus niger)	26	2
Bollons rattail	Bollons’ rattail (Caelorinchus bollonsi)	16	1
Cephalopod other	Squid & octopus	14	2
Cetacean other	Primarily sperm & pilot whales & dolphins	9	4
Dem fish pisc	Giant stargazer (Kathetostoma giganteum)	31	3
Elasmobranch invert	Primarily skates & dogfish	32	1
Elasmobranch pisc	Primarily semi-pelagic sharks	17	2
Epiben fish deep	Spiky oreo (Neocyttus rhomboidalis)	22	2
Epiben fish shal	Common roughy (Hoplostethus atlanticus)	8	2
Ghost shark	Dark ghost shark (Hydrolagus novaezealandiae)	30	2
Hake	Hake (Merlucciidae australis)	10	1
Hoki	Hoki (Macruronus novaezelandiae)	1	1
Invert comm herb	Paua (Haliotidae) & kina (Evechinus chloroticus)	21	2
Invert comm scav	Primarily scampi & crabs	33	2
Javelinfish	Javelinfish (Coelorinchus australis)	27	1
Ling	Ling (Genypterus blacodes)	11	2
Lookdown dory	Lookdown dory (Cyttus traversi)	20	1
Mackerels	Slender jack mackerel (Trachurus murphyi)	15	2
Orange roughy	Orange roughy (Hoplostethus atlanticus)	2	1
Pelagic fish lge	Southern bluefin tuna (Thunnus thynnus)	29	2
Pelagic fish med	Barracouta (Thyrsites atun)	7	3
Pelagic fish sml	Myctophids (Myctophidae)	4	3
Pinniped	NZ fur seal (Arctocephalus forsteri)	36	4
Reef fish	Blue cod (Parapercis colias)	24	2
Rock lobster	Rock lobster (Jasus edwardsii)	37	3
Seabird	Seabirds & shorebirds	6	4
Seaperch	Seaperch (Helicolenus spp.)	34	2
Shovelnosed dogfish	Shovelnosed dogfish (Deania calcea)	12	2
Smooth oreo	Smooth oreo (Pseudocyttus maculatus)	23	2
Spiny dogfish	Spiny dogfish (Squalus acanthias)	3	2
Warehou	Silver, white & blue warehou (Seriolella spp.)	13	3

DOI: 10.7717/peerj.9254/table-1

Table 2:

List of species groups modelled as biomass-pools for CRAM (McGregor et al., 2019b).

Zoo: zooplankton.

Species group	Description
Benthic Carniv	Benthic carnivores
Carniv Zoo	Planktonic animals (size 2–20 cm)
Carrion	Dead and decaying flesh
Deposit Feeder	Detritivores and benthic grazers
Diatoms	Diatoms (large phytoplankton)
DinoFlag	Dinoflagellates
Filter Other	Non-commercial benthic filter feeders
Gelat Zoo	Salps, ctenophores, jellyfish
Labile detritus	Organic matter that decomposes at a fast rate
Macroalgae	Macroalgae
Meiobenth	Benthic organisms (size 0.1–1 mm)
MesoZoo	Planktonic animals (size 0.2–20 mm)
Microphytobenthos	Unicellular benthic algae
MicroZoo	Heterotrophic plankton (size 20–200 μm)
Pelagic bacteria	Pelagic bacteria
Pico-phytoplankton	Small phytoplankton
Refractory detritus	Organic matter that decomposes at a slow rate
Sediment bacteria	Sediment bacteria

DOI: 10.7717/peerj.9254/table-2

Initially, we perturbed all species group initial conditions using the same scalar for all variables within each model run. The scalars we used were 0.5, 0.8, 0.9, 0.95, 1.05, 1.1, 1.2, 1.5. These scalars were chosen to cover a range from slight (±5%) to extreme (±50%) errors in the initial conditions.

As shifting initial conditions by the same amount may not give an indication as to how robust or sensitive the model is to mis-specification of the initial conditions where changes could vary in direction and magnitude, we next simulated multiple model runs, with the initial conditions scaled with some random variability. We scaled the initial conditions of each variable, sampling the scalar for each from a normal distribution, N (0, σ) with σ chosen based on how large we assumed a plausible change could be.

In total, we ran three sets of simulations, and repeated each set with and without fishing.

Set 1: All up or down. All species group initial conditions were scaled (numbers for age-structured, biomass for biomass-pool) with the same scalar for each run;

scalars ∈ {0.5, 0.8, 0.9, 0.95, 1.05, 1.1, 1.2, 1.5}

Set 2: High uncertainty. All initial conditions were scaled (numbers for age-structured, biomass for biomass-pool), with the scalars sampled from normal distributions with μ = 0 and σ set based on the informance ratings defined in McGregor et al. (2019b) (Fig. 1). Biomass-pool species groups were assumed poorly specified as these were not ranked in McGregor et al. (2019b).

Set 3: High keystone species. These runs only scaled the initial conditions of species groups likely to be most influential on the system. The species groups that ranked in the top 10 for keystoneness in McGregor et al. (2019b), and all biomass pool species groups were scaled using normally distributed scalars sampled with μ = 0 and σ = 0.25, giving 95% confidence intervals of ≈±0.5. All other species groups were unchanged (Fig. 2).

Characterising the species groups

If certain species groups appear to be more stable than others, we wanted to be in a position to investigate whether the more stable species groups have shared characteristics—for example are there links between sensitivities to changes in the initial conditions and how connected each species group is in the system, how abundant they are, how long they live, or some combination of these.

We characterised the species groups based on the base model presented in McGregor et al. (2019b) so we could test for links between these attributes of the species groups and sensitivities to changes in the initial conditions. We considered keystoneness, trophic level, biomass, animal size, lifespan, additional (background) mortality, number of trophic connections, and proportion of most dominant (‘top’) prey. All but the final three of these indices were available from McGregor et al. (2019b). The proportion of diet made up by most dominant prey, number of trophic connections, and the proportion of natural mortality that was made up of additional mortality were calculated for this study using R version 3.4.3.

Proportion of top prey

For each species group, we calculated the contribution to a predators diet from the single most dominant prey in their diet, as a proportion of biomass consumed, using averaged diets from the base model. The intent was to classify the extent to which each species group was eating as a specialist or generalist as they are modelled. It is possible for a predator to perform in the model as more of a specialist due to aggregation of species into groups—they could predate on several prey species that are modelled in the same species group. For each species group, we summed the prey eaten over the entire model region and all modelled years 1900–2015, then selected the largest proportion. Due to this averaging, a predator switching between two prey (for example) would appear similar to a predator consistently consuming those two prey. Hence, a predator would need to consistently predate on a given prey species group to be considered a specialist feeder for this metric.

Number of trophic connections

Trophic connections were calculated as ‘primary connections’ (predators and prey of the species group), ‘secondary connections’ (predators and prey of the predators and prey of the species group) and ‘tertiary connections’ (similarly) (illustration in Fig. 3). As sometimes a predator may eat a very small, negligible amount of a prey, we included a cut-off at 1%, such that a prey or predator was not included in the connections count if they made up less than 1% of the total prey consumed or predation pressure applied, respectively.

Additional natural mortality

There is the option in Atlantis to apply additional natural mortality either as a quadratic term, which is density dependent, or as a linear term (Audzijonyte et al., 2017a). The balance between additional natural mortality and mortality coming from dynamics within the model may affect the model’s stability. Higher levels of additional mortality reduce the strength of connections in the model, with 100% additional mortality effectively resulting in parallel single species models. Additional mortality was required for some species groups in the Chatham Rise Atlantis model that did not suffer sufficient natural mortality through predation, starvation or disease in the model to match estimates of mortality from the literature. For all age-structured species groups in this model, linear rather than quadratic mortality was applied as this is a close approximation to instantaneous natural mortality (M). When M is small, as it is when applied at small time-steps, e^−M can be approximated by 1 − x using the first two terms of its Taylor series expansion. Hence, if we take N_t to be the number of individuals at timestep t and N_{t + δ} to be the number at timestep t + δ where δ is small we get (1) $$N_{t+\delta }=N_t{e}^{-M_{\delta }}\approx N_t(1-M_{\delta })$$

As linear mortality, m_L, is applied at every timestep (12 h for this model), we can use M_δ to approximate m_L. This is, however, complicated by a temperature effect which is applied to m_L in Atlantis. Additional mortality in Atlantis are assumed related to metabolic rates, and hence are temperature dependent. The temperature effect is applied as a scalar (T_corr) calculated as a function of the current water temperature (T) (in a given cell at a given time) relative to a base temperature, set at 15 °C (Eq. (2)).

(2) $$T_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}=2^{(T-15)/10}$$

As temperature varies spatially and temporally, so does the scaled m_L. We calculated the additional applied mortality for each species group based on their spatial distribution, m_L values, and temperature corrections, using the median, upper and lower quartile, and 95% confidence intervals for the applied additional mortality to reflect the variability of temperature spatially and temporally. These were calculated for both juveniles and adults as m_L and spatial distributions were defined separately for these life stages.

Total realised mortality rates were estimated from the model by fitting an exponential decay curve to the proportions-at-age. By running the model with no fishing, the realised mortality consisted entirely of natural mortality, including sources within the model such as predation, as well as additional mortality from m_L. We then compared the total realised natural mortality with the range of additional mortality to estimate what proportion of natural mortality was coming from dynamics within the model, and what proportion was forced. We produced a weighted average for each species group that combined the proportions for adults and juveniles, weighted by the numbers of adults and juveniles respectively.

Results

Variability from initial conditions

Variability between simulations with perturbed initial conditions remained high for some species groups, while others converged to almost identical outputs. Generally, the biomass-pool species groups were more likely to have persistent high CVs between model runs (Fig. 4). Fishing sometimes reduced the between-model CVs for age-structured species groups, such as for hoki, but the effects of fishing were not apparent in any biomass-pool species groups—in biomass trajectories or CVs between model runs (Fig. 5 for two examples; Appendix A for the full set of figures). Responses of age-structured species groups to fishing were generally consistent across model runs. This included direct effects of fishing on a species (such as hoki, hake, orange roughy and ling), and predation-release responses (such as cephalopods and pelagic fish). Exceptions were invert comm herb (primarily paua and kina), invert comm scav (primarily scampi), dem fish pisc (primarily giant stargazer) and seaperch, which all gave varied responses with fishing included in the model.

Ecological indicators

Ecological indicators demonstrated variability from the perturbed initial conditions that generally persisted throughout the model simulations. However, the responses to heavy fishing from the mid-1970s were consistent across runs, with a decline in mean trophic level, a slight increase in diversity, and an increase in the ratio of pelagic biomass over total biomass (Fig. 6). There was a slight decline in mean trophic level from 1900 to 2015 in some of the unfished models, although the decline was approximately 0.02 of a trophic level over 100 years, so rather small. The ecological indicators present with more overlap between simulations than the individual species biomass outputs (Figs. 5 and 6).

Characterising the species groups

Proportion of top prey

Some diets consisted almost entirely of one species group, but many others did not have a dominant species with the ‘top’ prey making up less than 50% of the diet, and there was quite an even spread in between, with top prey making up around 50–70% of many diets (Fig. 7). Not all species groups predate, which is why some species groups (such as sediment bacteria, macroalgae) do not have a highest proportion of prey and are not included in the plot.

Number of trophic connections

The number of primary connections ranged from 1 through to 30, and with fairly even spread in between (Fig. 8). Most species groups were almost fully connected by the third level, and all species groups had at least 43 tertiary connections, of the 55 available species groups. Many of the species groups had more than 20 secondary connections, and those with fewer secondary connections generally had fewer primary connections. The number of secondary and tertiary connections are unlikely to be informative for stability between runs as there is little contrast.

Additional natural mortality

The proportion of natural mortality forced with additional mortality through the m_L term ranged from just over 0.8 for spiny dogfish down to zero for several species (Fig. 9). While baleen whales, cetacean other, pinnipeds and seabirds all have zero additional mortality through m_L, this does not mean their populations are entirely constrained due to mortality within the model, as these groups all migrate out of the model and their populations are restrained on re-entry into the model domain. Pelagic fish small (primarily myctophids), arrow squid, cephalopods other, and invert comm scav (primarily scampi) have all their natural mortality from sources such as predation within the model. Just over half (19/37) of the age-structured species groups had more than 80% of natural mortality forced as adults. Most age-structured species groups (31/37) had less than 50% of natural mortality forced for juveniles.

Modelling stability

GLMs fitted at each timestep

The models fitted at each year with all species groups combined (ALL) selected the interaction term ChaosAlt (the way in which the initial conditions were perturbed) and trophic level, and explained 43–53% of the null deviance (r²) (Table 6). The ChaosAlt:trophic level interaction term was also the most important explanatory variable for biomass-pool (BP) only species group models, explaining 40–41% of the null deviance, and age-structured (AS) only species group models explaining 31% of the null deviance (Table 6). For all models, ChaosAlt:trophic level was selected as an explanatory variable for most years (at least 102 out of 116) (Table 6). BP models consistently selected a second term, which was generally the interaction of the number of primary trophic connections and virgin biomass (B₀) and explained an additional 16–17% of the null deviance. The AS models had different explanatory variables selected at different timesteps, and these were also influenced by whether fishing was included in the models (Table 6; Fig. 10). The interaction between ChaosAlt and the number of primary trophic connections was the most consistently selected second explanatory variable for AS models (Table 6; Fig. 10). The AS models seemed to have a shift at around 1910, and explanatory variables selected prior to 1910 did not include ChaosAlt, but trophic level and informance were important (Fig. 10).

Table 6:

Ranges of r² values for selected explanatory variables from GLMs fitted to between-run CVs at each year 1900–2015 for each subset of species groups (BP: Biomass Pool species groups; AS: Age-structured species groups, and ALL: all species groups).

Each column reflects three subsets of model runs (without fishing; with fishing; both with and with-out fishing). The range in r² values reflects the year variability as well as whether or not fishing was included in the model. The numbers in brackets are the number of years (out of the 116 years for which GLMs were fitted) the explanatory variable was selected, with the range re-flecting the variability from whether or not fishing was included in the model. The AS results are also presented in Fig. 10.

Explanatory variable	BP	AS	ALL
ChaosAlt:Linf		0–0.01 (4–6)
ChaosAlt:NumL1cons		0.04–0.06 (38–59)
ChaosAlt:TL	0.4–0.41 (103–107)	0.31–0.31 (102–105)	0.43–0.53 (116–116)
Keystone:propAdM		0.01–0.01 (6–7)
Lifespan:Linf		0–0.03 (3–25)
NumL1cons:B0	0.16–0.17 (103–107)
propAdM:Linf		0–0.04 (4–29)
TL:B0	0.04–0.05 (9–13)
TL:Informance		0.01–0.02 (7–10)
TL:NumL1cons	0.02–0.02 (9–13)
TL:propAdM		0.01–0.01 (4–4)

DOI: 10.7717/peerj.9254/table-6

Final GLMs

The GLMs fitted to all data from 1910 to 2015 selected similar explanatory variables to the GLMs fitted at each timestep (Table 7). The interaction term ChaosAlt:TL was selected first for all models, and was the only term selected for the ALL model (all-species, with fishing and non-fishing runs included). The BP (biomass-pool) only species model also selected the interaction term NumL1cons:B₀. The AS (age-structured) species only model selected interaction ChaosAlt:NumL1cons whether fishing was included or not, and a third term, interaction NumL1cons:Informance was selected for the unfished AS model.

Table 7:

Explanatory variables selected and corresponding r² values for GLMs fitted to ALL (all species groups) model CVs, BP (biomass-pool species groups) only model CVs, and AS (age-structured species groups) only model CVs, using model outputs from 1910–2015, with fished and unfished versions for AS.

Model	ChaosAlt:TL	PrimCons: B0	ChaosAlt:PrimCons	PrimCons:Inf	Total r2
ALL	0.47				0.47
BP	0.38	0.53			0.53
AS	0.32		0.44		0.44
AS (fished)	0.33		0.45		0.45
AS (unfished)	0.31		0.44	0.54	0.54

DOI: 10.7717/peerj.9254/table-7

Note:

ChaosAlt, the set of runs, grouped by method for perturbing initial conditions and whether fishing was included or not; TL, trophic level; PrimCons, number of primary trophic connections; B0, virgin biomass; Inf, informance.

The Pearson’s residuals generally showed no concerning patterns against fitted values or explanatory variables for the final GLMs (Figs. 11–15). One exception was the residuals with respect to B₀ for the BP model, which suggested decreasing errors with increasing B₀, and a possible outlier (Fig. 12).

Higher trophic level was found to be associated with lower biomass CVs for all models and ChaosAlts (Fig. 16). CVs were generally lower for ChaosAlt ‘A’ and ‘B’, which were the model runs with all initial conditions shifted up or down and by the same scalar within each run. ChaosAlt ‘C’ and ‘D’, with initial conditions perturbed based on species group uncertainty, generally had slightly higher CVs across trophic levels (Fig. 16). This effect was also apparent in the interaction with primary connections in the AS model (Fig. 17). Biomass CVs were found to decrease with increased B₀ and with increasing number of primary connections for biomass pool species group (Fig. 17). The number of primary connections had the opposite effect for age-structured species groups, with more primary connections correlated with larger biomass CVs, although these CV effects were smaller (max. 11%) than for biomass pool species groups (max. 18%) (Fig. 17).

Discussion

Analysing sensitivities to initial conditions is an important part of developing complex models (Rabier et al., 1996; Rosati, Miyakoda & Gudgel, 1997; Payne et al., 2015; Cheung et al., 2016). If small perturbations to the initial conditions produce vastly different results, this may make interpreting results from the model challenging. Accounting for model uncertainty provides an envelope of model results, which tells us about the range of plausible outcomes rather than one possible instance. It is when the envelope is so wide that no result can be ascertained that it can be frustratingly un-useful, and it is important we are aware when this is the case. For example if scenarios exploring reduced fishing effort improved the general state of the ecosystem in some model runs, and deteriorated it in others, with all runs equally plausible, then we would be left none the wiser. It would be misleading to present results of only a subset or even a singular model run that does not adequately reflect the range of plausible outcomes.

We found the Chatham Rise Atlantis model was robust to initialisation uncertainty at the system level, in that we could perturb the initial conditions by small, and even quite large (up to 50%) changes, and the model produced very similar results with respect to ecosystem indicators. While the values of ecosystem indicators did retain some variability between model runs, the response to fishing was consistent, suggesting overall system dynamics were consistent under perturbed initial conditions. This puts us in a position to simulate scenarios using the Chatham Rise Atlantis model, including uncertainty of the initial conditions, and obtain an envelope of results with which to analyse and understand the likely responses of the Chatham Rise ecosystem.

While the system as a whole generally agreed within the range of results produced, the biomasses of some species groups varied between model runs more than others. The dynamics of some species groups appeared hyperstable as they promptly converged, while others retained variability between the runs, and for some the variability increased. We found the species groups that were more likely to have high biomass CVs were those of lower trophic levels. In nature, we expect to see more variability in the abundances of lower trophic level species, but most relevant field experts would likely suggest those patterns derive from variability within the environment (Dippner, Kornilovs & Sidrevics, 2000; Dippner et al., 2001; Molinero et al., 2008), which we are not applying in this study. If we combined varying the initial conditions with bootstrapping of the oceanographic variables, as carried out in McGregor et al. (2019b), we would likely see even greater variability in the lower trophic levels.

Another aspect of the trophic level effect on variability is the way in which we have modelled the species groups in the Chatham Rise Atlantis model. First, we have the difference between species modelled as biomass-pools and those modelled with age-structure. Biomass-pool representations are more dynamic as there is little/no delay structure built in—growth is pooled across its many forms (reproductive, somatic and otherwise), so can effectively occur instantaneously, unlike in age-structured groups where maturity may take years and specific events like spawning are restrained. Given biomass-pool groups are also generally lower trophic level (with naturally higher levels of productivity and turnover), the GLM fitted to CVs of all species groups could pick up trophic level as an explanatory variable that also accounts for this group structure. Within the age-structured species groups, trophic level could also be confounded with the proportion of additional mortality. The additional forced mortality would likely be a stabilising attribute, and the proportions applied were greater for the higher trophic level species, as these were the ones with less predation mortality in the system. That the stabilising aspect filters down through the trophic levels, with the lower trophic levels retaining variability, could suggest the extent to which this is a top-down controlled system.

The method used to perturb the initial conditions was found to be important in explaining the between-model run variability. The runs based on keystoneness did not result in the highest CVs, even though these runs perturbed the initial conditions of the species groups expected to have the greatest impact on the rest of the system. While Paine (1969) suggested keystone species have a stabilising effect on a system, it was more recently suggested to be more complicated than that (Mills, Soulé & Doak, 1993). In this study, the possible stabilising effect of keystoneness could be due to additional mortality applied to some of the high keystone species, and hence exerting a stabilising effect on the system. The runs perturbed based on uncertainty produced the greatest CVs. The effects of other explanatory variables, such as higher CVs for lower trophic levels, were consistent regardless of the method used to perturb the initial conditions. Hence, the method was not influential in how the system responded, only in how strongly it responded. It is possible the latter difference would diminish with a greater number of runs simulated for each set. In future simulations, perturbing the initial conditions based on uncertainty would seem appropriate, and should encompass the variability we would expect to see from other methods of perturbation.

One of the age-structured species groups that was most sensitive to initial conditions was the invert comm scav group (primarily scampi). When we account for uncertainty from initial conditions, the response of this group to heavy fishing is inconclusive. The heavy fishing on the system from the mid-1970s (Ministry for Primary Industries, 2017), some of which was targeted on scampi (Tuck, 2016), could easily be positive or negative for scampi based on these model results, and CVs for this species group remained high at just over 20%. In the base Chatham Rise Atlantis model (McGregor et al., 2019b), scampi were shown to respond to fishing in a very similar way to the fisheries stock assessment estimated biomass. The results here illustrate that the base model result for this species group, while convincing as it matched the fisheries models so well, was not representative of the many plausible results using this ecosystem model.

The effects of uncertainty from the oceanographic variables explored in McGregor et al. (2019b) had a greater range than those from perturbing the initial conditions. Similarly to the initial conditions uncertainty, the species most affected were also lower trophic levels. Diatoms had the highest CV at 79%, followed by carnivorous zooplankton CV at 46%. The effects of uncertainty from specification of the spawning stock recruitment relationship as it was applied to the small pelagic fish species group were explored in McGregor, Fulton & Dunn (2019a). These effects were seen right through to the ecosystem indicators, and would likely be more evident if more than one species group were directly affected in the study. Further work exploring the effects of parameter uncertainty ought to be carried out with the Chatham Rise Atlantis model, and these effects compared to those from the initial conditions, oceanographic variables and spawning stock recruitment already explored. One key area of parameters to explore is the predator/prey interactions, including the feeding functional response form and parameters. The effects of these will likely be noticed more as the model is used and hence taken away from its calibrated balanced state, where we might expect to see more variability in prey abundance.

An aspect not explored in this study is the effect of initialisation uncertainty on the model realised diets. There is plenty of scope for Atlantis diets to vary as they are the result of spatial and temporal overlap, gape-size limits, growth rates, feeding functional response, availability of other prey, predation from other predators, habitat refuge and prey preferences. As ecosystem models are generally developed to help understand flow-on effects within a system, understanding the effect of uncertainties on the species interactions could be important, and we recommend future work considers this aspect.

In the quest to provide meaningful and realistic results to simulations explored using complex ecosystem models, with high levels of uncertainty, we need to produce result envelopes, not single trajectories. It is important we move in the direction of simulating many instances of the model that account for its uncertainties, to understand how likely a given response is, and avoid presenting what may be errant or non-representative results. We know there is uncertainty in defining initial conditions of ecosystem models, so varying the initial conditions to reflect this uncertainty in model results is crucial. It is not the only area of uncertainty; there are many. Given the complexity of these models, exploring all possible uncertainties explicitly is unlikely to be tractable. It may be possible, however, to address subsets of uncertainty that encompass the broader range of the uncertainty of the model by targeting its key dynamics. The key dynamics of an ecosystem model generally consist of growth, recruitment, mortality, trophic connections, environmental effects, and initial state. Three of these (growth, mortality and trophic connections) relate directly to predation and consumption, and we could vary the feeding response function to explore the effects of uncertainties in these dynamics. Initial conditions were the topic of this study, and uncertainty from environmental effects were explored through bootstrapping the oceanographic variables in McGregor et al. (2019b). This leaves recruitment/productivity, for which we could vary the spawning stock recruitment parameters. The specifics of varying these will vary between models and systems, but accounting for uncertainty with respect to four main categories: (1) initial conditions; (2) environmental; (3) feeding functional response; (4) productivity/recruitment, is likely to cover the broad range for most systems and models.

Conclusions

The analyses presented here provided methods for testing and understanding the effects of initialisation uncertainty of an end-to-end ecosystem model. We present results of applying these methods to the Chatham Rise Atlantis model; an end-to-end ecosystem model of a deep-sea marine ecosystem of commercial importance since the mid-1970s. We found the lower trophic species groups to be more susceptible to initialisation uncertainty. We recommend results presented from ecosystem models are presented as envelopes of plausible or likely results that reflect uncertainty such as that from the initial conditions. This will help reduce the reporting of potentially errant results that can arise from single-run simulations that do not account for any uncertainty.

Appendix

Biomass trajectories (Appendix A).

Supplemental Information