\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{A modified niche model for generating food webs with stage-structured
consumers: The stabilizing effects of life-history stages on complex
food webs}



\author[1]{Etsuko Nonaka}%
\author[1]{Anna Kuparinen}%
\affil[1]{University of Jyväskylä}%


\vspace{-1em}



  
  \date{March 14, 2025}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
1. Almost all organisms grow in size during their lifetime and switch
diets, trophic positions, and interacting partners as they grow. Such
ontogenetic development introduces life-history stages and flows of
biomass between the stages through growth and reproduction. However,
current research on complex food webs rarely considers life-history
stages. The few previously proposed methods do not take full advantage
of the existing food web structural models that can produce realistic
food web topologies. 2. We extended the niche model by Williams \&
Martinez (2000) to generate food webs that included trophic species with
a life-history stage structure. Our method aggregated trophic species
based on niche overlap to form a life-history structured population;
therefore, it largely preserved the topological structure of food webs
generated by the niche model. We applied the theory of allometric
predator-prey body mass ratio and parameterized an allometric
bioenergetic model augmented with biomass flow between stages via growth
and reproduction to study the effects of a stage structure on the
stability of food webs. 3. When life-history stages were linked via
growth and reproduction, fewer food webs persisted while persisting food
webs tended to retain more trophic species. Topological differences
between persisting linked and unlinked food webs were small to modest.
Temporal variability of biomass dynamics and slopes of biomass spectra
were lower in the linked food webs than the unlinked ones, suggesting
that a life-history stage structure enhanced stability of complex food
webs. 4. Our results suggest a positive relationship between the
complexity and stability of complex food webs. A life-history stage
structure in food webs may play important roles in dynamics of and
diversity in food webs.%
\end{abstract}%



\sloppy


\textbf{Introduction}

A large body of research in the last several decades has investigated
potential factors that can promote the structural and dynamical
stability of complex food webs and their constituent populations. These
factors include hierarchically ordered feeding ), characteristic
predator-prey body mass ratios (Brose \emph{et al.} 2006), allometric
degree distributions of feeding links (Otto \emph{et al.} 2007),
compartmentalization (Stouffer \& Bascompte 2011), weak interaction
links including weak omnivory and reduced predation pressure at low
densities

(, pairwise negative correlation between interaction strengths , and
self-regulation (e.g., cannibalism, intraspecific interference, ), among
others (reviewed by ). More recently, studies have started incorporating
different types of interactions in complex food webs (multiplex or
multi-layer networks; ) to account for multiple types of ecological
interactions, such as mutualisms and parasitism, in which organisms
simultaneously engage in natural communities.

A ubiquitous feature of natural systems is that almost all organisms
grow in size during their lifetime and switch diets, trophic positions,
species interacting with, and habitats as they grow . Such ontogenetic
development introduces life-history stages and flows of biomass between
the stages through growth and reproduction to food webs, collectively
forming complex multi-layer ecological networks. Studies have shown that
ontogenetic diet shifts have far-reaching effects on competitive and
predator interactions, population dynamics, and community structure in
small food web modules . The persistence of consumers can be enhanced in
life-history structured communities through biomass overcompensation in
consequence of ecological asymmetry between different stages (e.g.,
juveniles are better competitors than adults; ). Such asymmetry,
however, can also be expected to destabilize populations by inducing
cohort cycles or alternative stable states without a predator . Research
on how these effects in small food web modules may scale up to an entire
complex food web is still in its infancy, and so are the tools to
generate life-history structured complex food webs in a biologically
justifiable manner.

Studies have reported the mixed effects of including a stage structure
on the stability of complex food webs . Rudolf \& Lafferty (2010) found
that, using static topological models of food webs, structural
robustness to species removal was lower with a stage structure than
without. They pointed out that species might be more sensitive to
resource loss when ontogenetic stages were sequential resource
specialists. Bland et al. (2019) used population dynamical models of
complex food webs and showed that non-stage-structured food webs lost
twice as many consumer taxa as stage-structured webs, while the
variability of biomass dynamics did not differ. Mougi (2017) also used
similar population dynamical models and concluded that species
persistence (the fraction of species persisting in a food web) increased
as the proportion of stage-structured species increased in the food webs
and that the effect was more pronounced in food webs with a greater
number of species and interactions. More studies are needed to elucidate
the role of a stage structure on persistence and stability and how it
may come about in complex food webs.

Rudolf \& Lafferty (2010) and Bland et al. (2019) used the niche model
(Williams \& Martinez 2000) to generate network topologies and split a
node into stages to create a stage-structured taxon (nodes represent
taxa, and interacting taxa are connected by links in ecological
networks). The niche model has a demonstrated ability to produce many
observed structural properties of empirical food webs despite its
simplicity and has been the most widely used food web structural model.
Splitting a node, as in Rudolf \& Lafferty (2010) and Bland et al.
(2019), can nontrivially modify the food web topology generated by the
niche model, likely compromising the desirable properties of the food
webs. Therefore, it is unclear how realistic the modified food webs in
these studies would still have remained after new nodes and links were
added to incorporate a life-history structure. Firstly, minimizing the
alteration of the network topology generated by the niche model is
desirable because it is capable of producing realistic food web topology
and because food web data resolved to life-history stages to verify the
topology of food webs with stage-structured taxa are currently very
scarce. Secondly, the niche model generates a ``trophic species,'' which
is a functional group defined to consist of one or more taxa (e.g.,
species, genus, ontogenetic stages) that share the same sets of
predators and prey . Life-history stages of a species are distinguished
for their distinct ecological roles, at least partly by their
characteristics related to feeding, so that a life-history stage can be
considered as a whole trophic species (not a fraction of it). Based on
these interpretations and the observation that ontogenetic diet shifts
are widespread in nature (Werner \& Gilliam 1984), a plausible
alternative approach is to instead group nodes generated by the niche
model to assemble a stage-structured taxon. This approach would allow
preserving mostly the original topologies of the food webs from the
niche model. No study has investigated this approach before.

We took the grouping approach to introduce a stage structure into
complex food webs. Following Bland et al. (2019), we applied the
allometric trophic network (ATN) model of biomass dynamics to the
stage-structured food webs on which we linked stages by biomass flow via
growth and reproduction. We motivated the food webs studied here from
aquatic communities at temperate and northern latitudes. It is well
established that consumer-resource interactions are hierarchically
structured largely by body size in aquatic communities because of the
indeterminate growth of fishes and gape-limited predation . We found
that food webs with stage-structured consumers tended to be less
dynamically variable and supported a greater number of species than food
webs with non-stage-structured consumers.

\textbf{Materials and Methods}

\emph{Introducing life-history stage-structured consumers to the food
webs generated by the niche model}

We built on the original niche model by Williams \& Martinez (2000; Box
1) and incorporated an additional algorithm to construct life-history
structured consumers by grouping trophic species based on the extent of
overlap between feeding ranges. As organisms grow in size during their
ontogeny, they experience changing competition, predation, and energetic
demands and may shift diets to maintain positive growth and minimize
mortality (Werner \& Gilliam 1984). Ontogenetic diet shifts among
life-history stages within a species are widely observed in nature
(Werner \& Gilliam 1984, Werner 1986), with diet ranges overlapping at
various degrees . Diet overlap is negligible in the case of diet shifts
associated with habitat shifts (e.g., riverine vs. marine) or
metamorphosis (e.g., aquatic vs. terrestrial), nested when organisms add
larger prey to the diet as they grow in size, or partially nested
because smaller prey are successively dropped from diet for energetic or
mechanical reasons . We applied this concept to construct food webs with
consumers with a life-history stage structure.

After obtaining food webs from the niche model (Box 1), we assigned two
measures of trophic position, the short-weighted trophic level
(\(T\)) and the prey-averaged trophic level
(\(T2\)), to each taxon . The short-weighted trophic level
is the average of shortest trophic level (\(T1\)) and the
prey-averaged trophic level (Williams \& Martinez 2000, 2004). The
shortest trophic level is equal to 1 + the shortest chain length from a
basal species to the consumer, and the prey-average trophic level is
equal to 1 + the mean trophic level of all the consumer's resources
(Williams \& Martinez 2004). An integer (\(N_{\text{fishes}}\); Table 1)
was uniformly randomly drawn from an interval between the desired
minimum and maximum numbers of stage-structured taxa. We assumed that
they were fishes (assume no stage structure in autotrophs and
invertebrates) and that fishes fed on at least one non-basal taxon
(i.e., \(T2>2\)). To create a stage-structured fish taxon, we
first selected a species with the highest \(T2\) (a ``focal
taxon'') that was greater than 2 (\(\text{Th}_{\text{fish}}\)), indicating that
this taxon ate at least one non-basal taxon (Fig. 1 and Fig. A1 for an
extended graphical example). We then identified taxa whose feeding range
maxima fell within the range of the focal taxon with the overlap of the
two feeding ranges greater than a specified minimum overlap
(\(\text{OL}_{\min}\)) of the union of the two and whose niche value was
smaller than and closest to the focal taxon's niche value, to become the
next lower stage. This stage became the next focal taxon, and we
repeated the steps until either the specified maximum number of stages
(\(\text{Nstage}_{\max}\)) had been assigned or taxa whose range maxima fell
in the range of the focal taxon with sufficient overlap ran out. When a
focal taxon did not have any overlapping taxa to choose from to form the
minimum number of stages (\(\text{Nstage}_{\min}\)), this taxon was
disqualified and another taxon was chosen in the same way as the current
focal taxon if other choices for the previous focal taxon were
available. If it was impossible to find the minimum number of stages for
the first focal taxon, it (but not the other ones that had been
subsequently considered) was reclassified as a taxon without a stage
structure and classified as an invertebrate. This occurred when a focal
taxon happened to have a small feeding range so that no range optima
fell in the range. The multiple taxa (nodes) selected in this procedure
collectively made up one stage-structured fish taxon and were removed
from the pool of available taxa. We then chose another focal taxon with
the next highest\(T2\) from the remaining taxa and repeated
the steps. We repeated these procedures until the chosen number of
fishes had been created or taxa with suitable range overlaps had run
out. If the minimum number of stage-structured taxa could not be
created, the food web was discarded. Predation of lower stages by higher
stages within a stage-structured taxon was interpreted as cannibalism.
Cannibalism within a stage and cannibalism of higher stages by lower
stages (very rare) were removed (cannibalism in non-structured taxa was
kept). Taxa that were not identified as autotrophs nor fishes were
identified as invertebrates (non-structured consumers). Therefore, nodes
in the network represented autotrophic trophic species, invertebrate
trophic species, or life-history stages of a fish. Hereafter, a node or
a taxon refers to a non-stage-structured species (invertebrates and
autotrophs) or a fish stage (Fig. 2). A species refers to an autotrophic
species, an invertebrate species, or a fish species that consists of
three or more life-history stages, while a stage refers to a fish
life-history stage.

\emph{Coupling stage-structured food webs and biomass dynamics in the
allometric trophic network (ATN) framework}

To simulate deterministic population dynamics of the species, we
employed a bioenergetic model in the allometric trophic network (ATN)
framework developed by Brose et al. (2006) and expanded by Bland et al.
(2019) to food webs with stage-structured fishes (see Brose et al. 2006
for a complete description). Consequently, we used many parameter values
and sub-models used in their work.

\emph{Body mass:} In this framework, body mass plays integral part in
determining bioenergetic parameter values . More specifically, the rates
of metabolism and maximum consumption are approximated by means of
body-mass scaling relationships (Yodzis \& Innes 1992). We calculated
relative masses of the taxa based on the short-weighted trophic position
\((T)\) in accordance with the theory of allometric
predator-prey body mass ratio (Brose et al. 2006). We set the body mass
ratio (\(Z\)) of fish predators and their prey to
10\textsuperscript{2.6} and of invertebrate predators and their prey to
10\textsuperscript{1.15} (; Table 1b). The function, body
mass,\(\ M=\ Z^{T-1}\), was used to define the body masses of
invertebrates and the terminal stages of stage-structured fishes. Hence,
the body masses were relative to those of autotrophs whose body masses
were defined to be equal to 1 (Brose et al. 2006; Bland et al. 2019). As
in Bland et al. (2019), to model the well-known pattern of fish growth
with time, we used a von Bertalanffy isometric growth curve to define
the body masses of lower stages (Table 2). We assumed that the
individuals of terminal stages reach 90\% of their asymptotic weight
(Bland et al. 2019). Although body masses in lower stages no longer
strictly conformed to the allometric body mass ratios, the median ratios
from our model fell near the modes of the empirical distributions (Fig.
3 in Brose et al. 2006; Figure A2).

\emph{Dynamical model:} The population dynamics within the food webs
were formulated as a multispecies consumer-resource model (Yodzis \&
Innes 1992; Williams \& Martinez 2004; Brose \emph{et al.} 2006; Bland
\emph{et al.} 2019). They were described by a set of ordinary
differential equations (ODE)

\begin{equation}
\frac{dB_{i}}{\text{dt}}=\overset{\text{logistic\ growth\ of\ autotrophs}}{\overbrace{g_{i}\left(1-\sum_{j\in autotrophs}\frac{B_{j}}{K}\right)B_{i}}}-\overset{\text{loss\ to\ grazing}}{\overbrace{\sum_{j\in consumers}{x_{j}y_{\text{ji}}B_{j}\frac{F_{\text{ji}}}{e_{\text{ji}}}}}}\nonumber \\
 \end{equation}\begin{equation}
\frac{dB_{i}}{\text{dt}}\underset{\text{metabolic\ loss}}{}+\underset{\text{dietary\ intake}}{}-\underset{\text{loss\ to\ predation}}{}\nonumber \\
 \end{equation}

where \(g_{i}\) was the intrinsic growth rate of autotroph
\(i\), \(K\)was the carrying capacity,
\(x_{i}\) was the metabolic rate of consumer\(i\),
\(y_{\text{ij}}\) was the maximum consumption rate relative to
metabolic rate, \(e_{\text{ij}}\) was the assimilation efficiency of
predator \(i\) eating prey \(j\),
\(f_{m}\) was the fraction of assimilated carbon lost for
maintenance, and \(f_{a}\) was the fraction of assimilated
carbon that contributes to biomass growth (see Table 1b for parameter
values). The model deterministically simulated the biomass dynamics
during growing seasons. \(F_{\text{ij}}\) was the functional response
of consumer \(i\) when dealing with prey
\(j\)

\begin{equation}
F_{\text{ij}}=\frac{\frac{\omega_{\text{ij}}}{\sum_{l\in resources}\omega_{\text{il}}}B_{j}^{q}}{B_{0_{\text{ij}}}^{q}+\sum_{k\in consumer}{\left(c_{\text{kj}}p_{\text{ik}}B_{k}B_{0_{\text{ij}}}^{h}\right)+\sum_{l\in resources}\left(\frac{\omega_{\text{ij}}}{\sum_{l\in resources}\omega_{\text{il}}}B_{l}^{q}\right)}}\nonumber \\
 \end{equation}

where \(\omega_{\text{ij}}\ \) was the preference of
consumeri\selectlanguage{ngerman} toward prey \(j\),
\(B_{0_{\text{ij}}}\) was the half saturation density for consumer
\(i\) eating prey \(j\),\(c_{\text{kj}}\)
was the predator interference competition coefficient of
\(k\) eating \(j\), and \(p_{\text{ik}}\)
was the fraction of resources of consumer \(i\) shared
with consumer \(k\). The values of\(\ B_{0_{\text{ij}}}\) and
\(c_{\text{kj}}\) varied among taxa and were taken from and Bland et
al. (2019, their Figure 1) with modifications (Table 1; also see ). The
parameters for interspecific or between-stage interference competition
were set to zero (i.e.,\(c_{\text{kj}}=0\) for \(k\neq i\)) for
simplicity (sensitivity to these assumptions were checked in the
sensitivity analysis). Previous studies that used the ATN framework for
aquatic systems (Brose et al. 2006, Boit et al. 2012, Bland et al. 2019)
differentiated the assimilation rates of consumers between non-basal and
basal species only. We added a rate for fish prey because fish is highly
effective food for fish growth (Table 1, and lowered the assimilation
rate for non-basal species (i.e., invertebrates) to have the average of
the two rates remain the same.

We added an ecologically plausible assumption that fishes preferred to
feed on fish over invertebrates and on invertebrates over autotrophs, if
they were included in their diets, to quickly grow beyond a size
vulnerable to predation and for higher fecundity. To achieve these
preferences in the absence of such empirical data, we set the
parameter\(\omega_{\text{ij}}\ \)such that fishes whose diets included both
autotrophs and animals fed almost exclusively on fish, to a lesser
extent on invertebrates, but not much on autotrophs (Table 1b).
Similarly, we assumed that invertebrates preferred invertebrates the
most, followed by fish, to autotrophs. Growth of fish depends on the
quantity and quality of food they eat, and shifting to piscivory
invariably increases fish growth rate . As fish grow, piscivory could be
necessary to meet energetic demands (Juanes et al. 2002). Also, because
optimal morphologies for different diets (e.g., planktivory, benthivory,
piscivory) are quite different, tradeoffs often arise and a diet
specializing on the most profitable is likely preferred (Persson 2002).
Herbivory by fish occurs mostly in tropics and is much less common above
55° latitude because the enzyme to digest plant material is not active
at low temperatures . If we assumed no preference of fish for prey items
(consumption proportional to relative availability), the majority of
fish would consume high proportions of autotrophs due to their high
abundance, an unlikely scenario in temperate and northern systems. If
prey taxa went extinct (\textless{} 10\textsuperscript{-6}), they were
removed from preference calculation.

The Hill exponent \(q\) of the functional response was set
to 1.8, higher than the value commonly used in previous ATN models
(1.2--1.5), to ensure sufficient dynamical stability in large food webs
(see Fig. A5 for sensitivity analysis; Williams \& Martinez 2004). The
high value of the exponent was desired especially because food
preferences of consumers increased energy flow higher up in the food web
and reduced stability of the food webs in the model

(Martinez \emph{et al.} 2006). Higher values of \(q\)
effectively converted the functional response closer to Holling type
III\(\left(q=2\right)\), which implicitly incorporates prey refugia,
other evasive behavior, or adaptive foraging .

\emph{Growth and reproduction:} Growth and reproduction from surplus
energy (dietary intake -- metabolic loss -- loss to predation; Eqn. 1)
were accounted for at the end of the growing season when the ODE model
was paused, which implicitly assumed that fishes all reproduced at the
beginning of each growing season . The fraction of mature fish in each
stage was determined by using a logistic function (Table 2). We assumed
that 50\% of individuals were mature halfway through to the terminal
stage. For example, if the taxon had five stages, about 50\% of
individuals were mature in Stage 3. We further assumed that fish in
immature stages invested all their surplus energy in somatic growth,
while mature fish allocated surplus energy to both growth and
reproduction (Kuparinen et al. 2016). The allocation to reproduction
(\(I\)) linearly increased with stage, and the terminal
stage allocated 20\% of surplus energy to reproduction (Table 2).
Therefore, the biomass of the first stage class produced through
reproduction was surplus energy multiplied by the probability of being
mature and reproductive investment. We used the Leslie matrix to shift
somatic biomass to the stage above via growth and to convert it to new
recruitment (Table 2). The model allowed phenotypic variability within a
stage such that some individuals did not grow enough during the
preceding growing season to be recruited to the higher stage. We assumed
that fish in the terminal stage reproduced without having energy surplus
in exchange for somatic mass (Wootton 1998). Each column added up to 1
in this formulation; therefore, there was no loss of biomass between
consecutive growing seasons (i.e., fish did not gain or lose mass or die
during winter).

\emph{Simulation design}

We generated 30,000 networks of 60 nodes and connectance equal to 0.15
containing between 2 and 6 stage-structured fishes with at least 3 and
up to 5 stages (Table 1a). We generalized the fixed numbers used by
Bland et al. (2019; three fish species with three stages). Food web
studies on northern temperate systems typically include two to three
species of fish, each of which has two to four stages (e.g., The minimum
overlap between feeding ranges to be qualified as consecutive stages of
the same stage-structured fish species was 20\% . We then ran the ATN on
each network from random initial biomasses uniformly distributed between
0.1 to 100 \(\mu gC/L\). We ran two sets of simulations on the
same networks. In one set, stages were linked via growth and
reproduction as described above (``linked''), while in the other set we
removed the growth and reproduction links between stages (``unlinked'')
to examine the effects of the additional biomass flow on the resulting
food web dynamics and persistence. The simulations were otherwise
identical (including the seed for the random number generator). Taxa
were considered extinct when the biomasses became \textless{}
10\textsuperscript{-6}\(\ \mu gC/L\), and if any of the taxon
reached the biomass of 10\textsuperscript{12} \(\mu gC/L\) the
simulation was terminated (Table 1b). These values were many orders of
magnitudes smaller and larger, respectively, to the mean total fish
biomass (10\textsuperscript{1.1} \(\mu gC/L\) ). At the end of
each generation, fishes that retained only older stages but not younger
ones for more than 10 generations were removed as extinct. This happened
in some simulations because the biomasses of older stages without
younger ones lingered although they were to decay over time. Each year
consisted of 90 time-steps, representing one growing season, followed by
a non-growing season where reproduction and growth were accounted for.
One time-step corresponded to the generation time of the autotrophic
species the growth rates of other taxa were normalized to the time scale
of the basal species, as typically done in the ATN models; Brose et al.
2006). Food webs were regarded as persistent if simulations lasted for
1000 growing seasons, which was sufficient for transient dynamics to die
out, with all the nodes of persisting species connected in one network,
no invertebrates without predators, no species without food (i.e.,
transients did not completely die out in a handful of webs), and at
least one fish with 3 or more stages for the linked case or at least 3
fish taxa for the unlinked case. We note that we regarded food webs as
persistent if they retained fishes and met other criteria we just stated
(Fig. 3) to the end of the simulations, not by the proportions of taxa
(nodes) persisted as typically reported by food web studies.

\emph{Analyzing model outputs}

We first assessed various food web properties of the 30,000 networks
(``original webs'') to check whether the generated food webs were
reasonable. Then, we identified the linked and unlinked food webs that
persisted for 1000 growing seasons (i.e., years) and passed the
abovementioned criteria (Fig. 3), and retained only persisting taxa in
these webs (``ATN-filtered webs''). We then examined food web properties
of the ATN-filtered food webs using the last 100 years of the
simulations. We calculated the mean of the CVs of biomasses of fish
stages to characterize dynamical stability of food webs and the number
of nodes and link density to characterize food web complexity. We also
compared the unlinked and the linked webs for total fish biomass, the
mean body masses of fishes weighted by relative abundance, the number of
fish stages, and the maximum trophic levels and median predator-prey
body mass (PPMR) ratios for fish stages. We computed a surrogate of
total interaction strengths by summing up all energy flows into fish
stages from their resources, CVs and skewness of the energy flows, which
were individually normalized by the total gain of the recipient fish ,
and the slopes of biomass spectra by fitting a linear model to estimate
biomass at trophic levels as a function of body mass . The presence of
weak interactions is theorized to stabilize food webs . Less negative
slopes of biomass spectra indicate less bottom-heavy food webs (Treblico
et al. 2013), which tend to be less stable than more bottom-heavy food
webs . Some measurements were log\textsubscript{10} transformed for
interpretability. We computed 10000 bootstrap estimates of the mean
effect sizes of a stage structure (linked -- unlinked) for these metrics
by taking the differences between the means of resampled values from
persisting linked and unliked food webs.

We checked the sensitivity of the model outputs to major input
parameters (Table 1). We ran 10000 simulations for each parameter
variation (16 sets in total) using the same food web topologies as the
baseline simulations. We computed 1000 bootstrap estimates of the mean
effect sizes of a stage structure (linked -- unlinked) for these metrics
by taking the differences between the means of resampled values from
persisting linked and unliked food webs and compared the mean effect
sizes with those from the baseline simulations.

\textbf{Results}

\emph{Structural properties of the stage-structured food webs}

The structural properties and characteristics of the 30000 original
stage-structured food webs generated by the modified niche model were
quantified on the unweighted networks (Fig. A3). The majority of the 60
nodes were invertebrates, followed by 8 to 30 life-history stages of
fishes (2 to 6 species as specified) and 3-12 autotrophs. Between 6 and
30 fish stages were piscivorous (their diets included fish), while 3 to
28 of them had autotrophs in their diets. Fishes were almost all
cannibalistic at some stage. Omnivores (feeding at 2 or more trophic
levels) were abundant (45 to 50 species). Each fish stage had about 16
prey species and 7 predators on average. The mean maximum trophic level
was near 5. Almost all interactions involved intermediate taxa (taxa
with both prey and predators). About a half of the webs did not have a
top predator (taxa without a predator), while the rest had one top
predator. The sensitivity analysis showed that these patterns remained
almost identical when \(\text{OL}_{\min}\) was reduced to 10 \%, while the
number of fish stages decreased along with other measures as a direct
consequence of having fewer fish stages when\(\text{NStage}_{\min}\), and
\(\text{NStage}_{\max}\) were decreased (Fig. A3).

\emph{The impact of a stage structure on the food webs}

Out of 30,000 food webs, 5865 linked (19.6\%) and 9099 unlinked webs
(30.3\%) persisted for 1000 growing seasons and passed the criteria, of
which 2628 webs were shared. Assessed by a set of network structural
metrics and characteristics quantified on the unweighted networks after
removing extinct taxa (``ATN-filtered''; note that taxa that went
extinct could be different between the linked and unliked webs), the
persisting food webs in both scenarios largely overlapped in terms of
the majority of the properties measured (Fig. A4). The number of
persisting nodes ranged from 15 to 60, of which invertebrates comprised
the majority and 3 to 30 nodes were life-history stages of fishes. The
majority of the food webs did not have a top predator.

The distributions of measured food web properties between the linked and
unlinked food webs overlapped substantially (Fig. 4a). Nonetheless, the
bootstrap differences in means between the linked and unliked food webs
(linked -- unlinked webs) highlighted that some of the mean effect sizes
were noticeable (Fig. 4b); the mean CVs of fish biomasses was about
3.3-orders of magnitude smaller, the number of nodes was greater by
about 6, the slope of biomass spectra was 0.17 smaller, and the CV and
skewness of normalized fish energy gain were 0.17 and 2.35 larger,
respectively, in the linked webs than in the unlinked webs on average.
These suggested that, when stages were linked relative to when unlinked,
on average fish biomass dynamics were less variable, the biomass
pyramids were more bottom-heavy, food webs sustained several more
species or fish stages, and fish energy gain was more variable among
links and dominated by small gains (indicated by the positive skews).
Taken together, the results suggested that the life-history stage
structure of fishes on average moderately promoted the stability of and
diversity in these complex food webs, once they were able to persist.
The sensitivity analysis showed that the patterns in the bootstrap
differences in mean effect sizes remained in the same direction
(negative or positive) as the baseline simulations across the
sensitivity analysis simulations except several cases (Fig. A5). The
bootstrap differences in mean effect sizes were qualitatively deviated
from the baseline runs for a few measures when the maximum and minimum
numbers of fish stages\(\left(\text{NStage}_{\max},\ \text{NStage}_{\min}\right)\) were smaller (\#1), full
interference competition was included (\#13), or when \(q\)
was reduced to 1.2 (\#15). The number of persisting webs was markedly
smaller (\textasciitilde{}25\%) when \(q=1.2\) (data not
shown).

\textbf{Discussion}

\emph{The modified niche model}

To introduce life-history stages into food webs, we developed additional
algorithms that modified food webs from the niche model by Williams \&
Martinez (2000). Food web topologies generated by other food web
structural models can also be used as long as they assign feeding
hierarchies and ranges to every taxon, such as the variants of the niche
model (generalized niche model ; the relaxed niche model ; the minimum
potential niche model ). Our method utilized the concepts of ontogenetic
diet shifts and niche overlap among ontogenetic stages (Werner \&
Gilliam 1984; Werner 1986) to identify life-history stages and
heuristically assembled the specified number of taxa with a stage
structure. It can be adapted to other situations where, for example,
stage-structured taxa feed higher in the trophic level (e.g., set the
minimum trophic level \textgreater{} 3 to become a fish or exclude
consumers feeding on autotrophs from the pool of fish candidates) or
feeding range overlaps are smaller or larger. The outputs can be fed
into the ATN framework or other dynamical models that can accommodate
biomass flow via growth and reproduction.

We grouped trophic species created by the niche model to assemble a
stage-structured fish taxon, unlike the previous models where a trophic
species was split into stages (Rudolf \& Lafferty 2010; Bland et al.
2019). As a result, our method generated food webs that largely
preserved topologies (i.e., besides removing rare within-stage and
reverse cannibalism) produced by the niche model, which has been shown
to reproduce empirically observed food web properties (Williams \&
Martinez 2000, 2008). The previous methods introduced new nodes and
links, likely compromising the merit of using the niche model. Our
approach also agrees with the method employed by Williams \& Martinez
(2000, 2008) to evaluate the niche model's performance, where some of
the empirical data they used distinguished different stages of the same
species (e.g., larval/young-of-year and adult fish in Little Rock Lake,
Ythan estuary, and Chesapeake data). Our approach hence followed from
the definition of trophic species, a group of taxa sharing predators and
prey, from the common phenomena of the ontogenetic diet shift, and from
the fact that the niche model creates trophic species. What constitutes
a trophic species should depend on the level of aggregation appropriate
for a given study. Because we were interested in trophically distinct
roles of ontogenetic stages on food web dynamics (Werner 1986), it was
both convenient and reasonable to interpret trophic species as
ontogenetic stages and group multiple trophic species into a stage
structured species. It appears to be a great advantage to minimize
alteration of food webs obtained from the niche model. As a by-product,
we also eliminated the convoluted steps to assign niche values to newly
created nodes in the method by Bland et al. (2019). We think that our
approach improves and simplifies their method, making it more
conceptually accessible to food web researchers.

In our results, a far greater number of food webs with unlinked stages
persisted than those with linked stages. Linked stage-structured food
webs were qualified with more stringent criteria (namely, higher stages
cannot persist without lower stages for more than 10 generations vs.
independent stages; at least one fish with 3 or more stages persisting
vs. at least any 3 fish nodes persisting). Once the food webs were
persisting, linked life-history stages stabilized food webs relative to
when they were unlinked, as indicated by the lower variability of
biomass dynamics (CV) and higher numbers of taxa (nodes) persisting in
the food webs (Fig. 4). The relative frequency of food webs with
oscillating biomass dynamics (i.e., higher CV) was higher when stages
were unlinked. In contrast, the linked and unlinked food webs in Bland
et al. (2019) behaved similarly in terms of these measures. Furthermore,
the mean of the CVs of energy flow into fishes was modestly higher in
the linked webs on average, mirrored by greater positive skewness, which
indicated that the linked webs contained more weak links than the
unlinked webs. A similar difference in CVs of energy influx resulted in
approximately a 5\% change (an increase in webs of 20 species, while a
decrease in larger webs) in the proportion of stable webs in Gross et
al. (2009). Having weak interactions is one of the key properties that
can increase stability of food webs . Also, we observed that the linked
webs had lower slopes of biomass spectra and hence exhibited more
bottom-heavy biomass pyramids than did the unlinked webs. Bottom-heavy
biomass pyramids tend to relate to dynamically stable consumer-resource
dynamics, while top-heavy biomass pyramids tend to suggest unstable
dynamics . Therefore, the stabilizing effects of life-history stages
that we saw in our simulations appear in agreement with what current
food web theories predict.

Our method and the method by Bland et al. (2019) also differed somewhat
in modeling demographic shifts via growth and reproduction at the end of
growing seasons. The differences were in how surplus energy were dealt
with and in the proportion of the biomass of the terminal fish stage to
be transferred to the first stage. Thus, the differences between our
results and those of Bland et al. (2019) may not be attributable only to
how life-history stages were constructed (grouping nodes vs. splitting a
node). Further research should systematically explore how a life-stage
structure can affect food web stability. Our method can serve as a tool
to generate biologically justifiable stage-structured food web
topologies to facilitate such explorations in future studies.

We noticed that the original niche model by Williams \& Martinez (2000)
produces many consumers that include autotrophs in their diets. In
temperate and northern regions, fishes feeding on autotrophs are
uncommon because of low activity levels of digesting enzymes . We added
this realism through prey preferences in the dynamic model. In
simulations, fishes consumed little autotrophic biomass as a
consequence, despite including autotrophs in their diets. We also
noticed that persisting webs often did not have a top predator based on
their topologies (Fig. A4), which implied that food webs with top
predators tended to be dynamically unstable. In natural systems, there
may be no true top predators as even the adult stages of piscivorous
fishes could be eaten by vertebrate and invertebrate predators (e.g.,
birds, octopuses) or parasitized (e.g., scale eaters, leeches,
pathogens). Alternatively, we could have retained in the ATN filtered
webs only interactions that were substantial enough to be observed in
the fields. If we removed interactions that contributed a very small
fraction to the consumers' total energy intake (\(<10^{-4}\%\)) with
the specified diet preference, 70\% of persisting webs had at least one
top predator (data not shown).

\emph{Life-history structure and food web stability}

How a life-history structure may affect the stability and persistence of
complex food webs has not been much studied. It is not immediately clear
whether it is stabilizing or destabilizing based on existing theories.
We can expect multiple aspects of stage-structured populations to
contribute to instability. As discussed by Rudolf \& Lafferty (2010),
when stage classes have smaller subsets of the feeding range than the
species as a whole, resources essentially become less substitutable,
especially when the overlaps between feeding ranges are small. Thus, if
resources for one stage become scarce, the persistence of the entire
species is greatly endangered unless growth and reproduction can
constantly replenish the dwindling stage. Similarly, as stages become
more specialized, consumer-resource interactions may become less diffuse
and some of the remaining interactions may strengthen. Because weak
interactions tend to stabilize trophic interactions (McCann 2012;
Gellner \& McCann 2016), specialized stages likely reduce stability of
food webs. Also, a stage structure introduces delays and asymmetry
between stages into population models, both of which are known to often
cause population instability in the forms of cohort cycles and
alternative stable states . Therefore, the odds seem to be against
increased food web stability by introducing life-history stages.

Stages in structured populations can subsidize dwindling stages through
growth and reproduction, which is probably one of the main reasons why a
stage structure in food webs could enhance the persistence of
stage-structured populations and other dependent populations.
Furthermore, biomass flow via growth or reproduction between competing
stages with overlapping diet might moderate the destabilization effects
of exploitative competition. showed that the exploitative competition
module reduced food web persistence as it increased in frequency in
dynamical models of complex food webs. This effect of diet overlap may
appear contradictory to the result from Rudolf \& Lafferty (2010), which
showed that, when feeding niches were overlapping by more than about
30\%, the inclusion of stages increased the robustness of food webs.
Because they studied the robustness of static food webs (only topology,
no dynamics), diet overlap reduced reliance of a stage-structured
population on any particular resource, as the authors explained. In
dynamic models, exploitative competition can ensue and drive one of the
competitors and possibly other populations to extinction, but if the
competitors are ontogenetic stages of the same species, biomass flow
between the stages could alleviate competitive exclusion. In the study
by Stouffer \& Bascompte (2010), the frequency of the tritrophic food
chain module had positive effects on persistence in large food webs. In
a sense, life-history structured populations contain a biomass flow
chain inside. We conjecture that this might also contribute to food web
stability. In addition to the possible adverse effects on stability we
discussed above, ontogenetic asymmetry may also help some populations
persist in food webs. showed in stage-structured food web modules that
persistence of consumers could be promoted in communities with
stage-structured prey through emergent facilitation due to biomass
compensation in the prey population. It seems reasonable to state that
the effects of life-history stages on the stability of complex food webs
are complex and contingent on the balance of the effects of different
processes.

\textbf{Conclusions}

Life-history stages constitute part of biological diversity and increase
complexity of food webs. As a large majority of organisms grow in size,
often over orders of magnitude, during their lifetime and experience
various degrees of ontogenetic diet shifts (Werner \& Gilliam 1984),
life-history structures are important to be considered in studies on the
stability of complex food webs. In this study, we demonstrated a
positive relationship between the complexity and stability of complex
food webs; food webs with stage-structured taxa exhibited lower biomass
variability and supported more taxa than did those with unlinked stages.
These results are qualitatively in agreement with the findings by Mougi
(2017). For aquatic systems and fishes in particular, ontogenetic stages
are well recognized and studied so that including life-history stages
explicitly in models can facilitate linking theory and data.
Practically, including separate stages makes it more mechanistic and
straightforward to implement allometrically scaled functions or
parameters and differences in behaviors among stages. For example,
simulating size-selective fishing and the evolutionary impacts of such
fishing on the population dynamics of exploited species in food webs
becomes more straightforward once a life-history structure is explicitly
incorporated (e.g., Kuparinen et al. 2016). Moreover, our work
contributes a way of incorporating another aspect of interaction
diversity via life histories to the growing research on multi-layered
networks . Biomass flow via growth and reproduction forms networks of
energy transfer analogously to consumer-resource interactions. Research
on multi-layered networks has so far revealed that non-trophic
interactions (thus interaction diversity) can ameliorate or degrade the
stability of trophic interactions and the persistence of species
(reviewed by Kéfi \emph{et al.} 2018). Interestingly, showed that
network structures known to stabilize mutualistic interactions became
less effective when combined with trophic interactions in a
multi-layered network. Tritrophic food chain and omnivory modules have
been shown to stabilize complex food webs (Stouffer \& Bascompte 2010),
and it will be instructive to examine if they still do so when embedded
in complex food webs including ontogenetic biomass flow.

\textbf{Literature cited}

\emph{Box 1. The niche model by Williams \& Martinez (2000)}

The niche model requires the number of trophic species (the number of
nodes, \(S_{0}\)) and connectance (the fraction of realized
feeding interactions out of all potentially possible,
\(C\)) as input parameters (Table 1). It hierarchically
ranks species according to the ``niche value,'' \(n_{i}\),
randomly drawn from a uniform distribution and assigns a feeding range
to each species as follows. The range size,\(r_{i}\)
(\(i\) is the index for taxa), is determined by first
drawing a random variable, \(\eta_{i}\), from a beta distribution
calibrated to obtain the desired connectance and then multiplying
\(\eta_{i}\) by the niche value for \(i\ (r_{i}=\eta_{i}n_{i})\). The center
of the feeding range, \(c_{i}\), is randomly chosen from a
uniform distribution in \(\left[\frac{r_{i}}{2},n_{i}\right]\), and the range is then
determined as\(\left[c_{i}-\frac{r_{i}}{2},c_{i}+\frac{r_{i}}{2}\right]\). Therefore, species with larger niche
values tend to have larger feeding ranges. The ranges are set such that
cannibalism is allowed (i.e.,\(n_{i}\) can fall in the range
of \(i\)). All the taxa whose niche values fall in the
feeding range of another are regarded as the prey of the latter. The
taxa with no prey are identified as basal taxa (i.e., autotrophs). The
taxon with the lowest niche values is designated as an autotroph. We
discard disconnected webs, webs with connectance beyond a given
tolerance level (\(C_{\text{error}}\)), and webs with taxa not connected
to a basal taxon.

\par\null

\textbf{Data accessibility statement}

The Matlab code for the modified niche model and the ATN model used in
this study is available in the Dryad Digital Depository https:/ XXX

\textbf{Competing interests statement}

The authors declare no competing interests.

\textbf{Author contributions}

A.K. and E.N. conceived the ideas. E.N. led the development of the model
to generate life-history structured food webs, wrote the code,
implemented and analyzed the models, interpreted model outputs, and
wrote the manuscript. A.K. provided advice on the interpretation of
outputs and critically commented on the manuscript. Both authors gave
final approval of the version to be published.

\textbf{Acknowledgements}

The authors thank Dr. Silva Uusi-Heikkilä for critical and insightful
comments on the manuscript. This study was funded by the Academy of
Finland (project grant 317495 to A.K.), Natural Sciences and Engineering
Research Council of Canada (NSERC; Discovery Grant to A.K.) and the
European Research Council (COMPLEX-FISH 770884 to A.K.). The present
study reflects only the authors' views, and the European Research
Council is not responsible for any use made of the information it
contains.

\selectlanguage{english}
\FloatBarrier
\end{document}

