Preliminary Assessment of a Coupled Dynamic-Energy Budget and Agent-based Model (DEB-ABM) for Predicting Individual and Population-Level Dynamics: A Case Study on Anchovy, Engraulis japonicus

Baochao Liao1,2, Xiujuan Shan2,3* and Yunlong Chen2,3 1Department of Mathematics and Statistics, Shandong University, Weihai, Shandong 264209, China 2Function Laboratory for Marine Fisheries Science and Food Production Processes, Qingdao National Laboratory for Marine Science and Technology, Qingdao, Shandong 266237, China 3Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Qingdao, Shandong 266071, China Article Information Received 22 February 2020 Revised 04 April 2020 Accepted 23 April 2020 Available online 04 September 2020


INTRODUCTION
D ynamic energy budget (DEB; also known as the Kooijman-Metz DEB) theory links the physiological processes of individual organisms into a single framework (van der Meer, 2006;Kooijman, 2010). The theory describes the energy and matter fluxes within an individual, and between it and its environment (Kooijman, 2010;Marn, 2016;Agüera et al., 2017). Traditional DEB bioenergetic models are applied at the level of the individual, and can be used to deduce the energy input from individual-level growth and reproduction O n l i n e

F i r s t A r t i c l e
of a population of individuals, with each individual following an energy budget model (Railsback and Grimm, 2011;Sibly et al., 2013). This combined model can be used to predict how the energy dynamics at higher levels of biological organisation emerge from individual-level processes (Sibly et al., 2013). DEB models have been coupled with ABMs to study fish population dynamics in a large body of work Sibly et al., 2013;Grossowicz et al., 2017;Desforges et al., 2017;Smallegange et al., 2017). The pelagic Japanese anchovy (Engraulis japonicus; Temminck and Schlegel, 1846) is an important prey species and plankton feeder; it is also the most abundant fish in Chinese waters, including the Bohai (BH), Yellow (YS), and East China (ECS) seas (Jin, 2004;Wang et al., 2003). However, increased fishing effort in Chinese waters in recent decades has led to dramatic changes in the population structure of the anchovy fishery. This small, rapidly reproducing fish has a life span of approximately 4 years (Iversen et al., 1993); it acquires a portion of its reproductive energy reserves during a productive period a few months before spawning season (Zhao et al., 2003;Zhao et al., 2008). Individuals grow rapidly during their first year, and mature after their first winter (Zhu and Iversen, 1990;Chiu and Chen, 2001). Their life-cycle includes five stages: egg, yolk sac larva, exogenous feeding larva, juvenile, and adult (Wan and Bian, 2012;Wang et al., 2003).
In this present study, the DEB-ABM model analyzed the contributions of different energy sources and calculated the energy budget required to maintain growth and reproduction of a small and rapidly reproducing species E. japonicus in Chinese seas. Our aims are to: i) find a parameterized method and a validation of model output; (ii) present a novel model for investigating the dynamic population-level energy budget of anchovy; and iii) determine the most relevant quantitative variable to explain energy budgets relevant to growth and reproduction of anchovies across study areas.

Individual-level
The DEB model predicts individual energy allocation, acquisition, and use for growth and reproduction (Kooijman, 2010;Nisbet et al., 2012). Food intake was controlled by a scaled (dimensionless and Holling type-II) functional feeding response (f) curve that ranged between 0 (i.e., no food intake) and 1 (i.e., maximum food intake); a fixed fraction of κ was allocated to maintenance and growth, and the remaining energy (1−κ) was utilized for development and reproduction ( Table I). The food uptake was proportional to the surface area of the larva and to food density. The ingestion rate was equal to ṗ X =f(X)*{ṗ X }*V 2/3 , where f(X)=X/(X+K). The four state variables (fluxes, units: J day −1 ), including the Arrhenius temperature (T A ); the volume-specific somatic maintenance ([ṗ M ]); energy conductance (ὺ); allocation coefficient (κ); fraction of energy allocated to fixed growth (κ G ); and maturity maintenance rate coefficient (κ J ). DEB dictated that energy is assimilated (ṗ A ) from food and transferred into reserve (E). The fixed fractions (κ and 1−κ) of the catabolized energy (ṗ C ) are allocated to either soma (E v , units: J) or maturation/reproduction (E R ). Increases in anchovy size (V) regulated transitions between developmental stages. All energy allocated to the reproductive buffer before puberty (V <V P ) was lost to maturation (E R = 0). DEB model was parameterized based on the biometric data of juveniles and adult E. japonicus (SL range: 5.8-17.2 cm) collected between 2001 and 2004 (Wan et al., 2004;Zhao, 2006). Anchovy water mass W ww was predicted by the following equation: W ww = exp (−3.19+0.38*l + 0.59* log (e)−0.035 *l * log (e)); this length-weight relationship is typically expressed as W = a×L b . Larval growth data were consisted of individual specimens (Total length: 1.8-5.4 cm; W ww : 0.018-1.32 g). Larval growth and proximate composition data for anchovies were used to estimate the values of the state variables as follows: E R0 =0 J; E GAM0 =0 J; E 0 =0.11 J; E V0 =0.025 J; and V 0 =0.00013 cm 3 . We set the length-and weight-at-first-feeding (day 4) values to 0.292 cm and 0.031 g, respectively (Wan et al., 2002; Table II. Condition index (K ful ) Model output (g cm 3 ) Gonad somatic index (GSI)

Population-level
Population-level dynamics describe those of a population of individuals, where each individual follows an energy budget model Sibly et al., 2013). The i-states in our model follow DEB characteristics in their size (V), reproductive output (E R , E GAM , and E batch ), and energy reserves (E) (Pecquerie et al., 2009Sibly et al., 2013Grossowicz et al., 2017), as follows: Where the variables V i , E i , R i , L i , and WT i were associated with cohort i and simulated with DEB equations; h i represented the harvest rate coefficient; m i was the natural mortality rate coefficient; N i was the number of individuals in cohort i; P was the cumulative harvest of cohorts; N was the total number of individuals; and S was the weight of all cohorts after summation to estimate the standing stock at each time step. The total recruitment (R) was calculated as a function of spawning stock biomass (SSB) (Zhao et al., 2003). The characterized state variables of the individual were included reserves (E), structure (V and E v ), reproduction (E R ), and gametes (E GAM ; Table I) (Regner, 1996;Boussouar et al., 2001). The sensitivities of the intrinsic model-specific parameters (deterministic, nondeterministic, positive, and negatively correlated) were investigated using a traditional one-parameter-ata-time analysis (OPAT), in which each model parameter was varied separately with ±10% white noise (CV). We calculated a relative estimate error (REE) to compare results among CVs (Chen et al., 2010). The REE was O n l i n e

F i r s t A r t i c l e
A Coupled DEB-ABM Model to Predict Population-level Dynamics 5 evaluated by calculating a simple sensitivity index (SSI) as follows: Where W t 0 is the total estimated wet body weight in simulation (predicted with the new parameter value at time t); W t 1 is the total actual wet body weight in simulation (predicted with the standard simulation at time t); k is the number of simulated days; L t 0 is the total body length (predicted with the new parameter value at time t); and L t 1 is the total body length (predicted with the standard simulation at time t). The effects of external factors, including sea surface temperature (SST), natural mortality (N 0 ), and food availability (X), on the growth rate of individual anchovies (or anchovy populations) was evaluated using a number of scenario-based simulations.

Individual traits
The Arrhenius temperature (T A ) is an expression of the effects of temperature on biological reaction rates. Here, we obtained a mean T A value of 9800 ± 835 K, which was determined by plotting ln(1/D) against 1/T, where D and T are egg development time and temperature, respectively. We estimated two different Arrhenius relationships for the moving average (MA) of ingestion rate and respiration rate (Fig. 1). The average length-and weight-at-age of an individual anchovy collected in May predicted by our simulation were consistent with the actual length-and weight-at-age recorded during the same period by our surveys (2000)(2001)(2002)(2003)(2004). We calculated energy density as the ratio of the total energy reserves to the total anchovy wet weight. By simulating differences in hatching dates in a seasonal environment, over several years, our model produced variability in anchovy body length and weight (Fig. 2). To convert state variables into units of weight (gW dw ), we used separate conversion factors for structures (E V ; μ V = 19.9 KJ g −1 ), reserves (E and E R ; μ E = 35.2 KJ g −1 ), and gametes (E GAM ; μ G = 33.2 KJ g −1 ). As anchovies have different shapes as juveniles and adults, separate shape coefficients were estimated by fitting the weightlength relationship as W ww 1/3 /L. The coefficients δ Adult and δ larvae were 0.169 and 0.154, respectively. The threshold structural volumes at first feeding (V ab ), at metamorphosis (V morph ), and at puberty (V p ) were estimated to be 0.000164, 0.53, and 5.43 cm 3 , respectively. The simulation was designed to validate the predictions of our model under conditions of prolonged starvation. The shape coefficient δ m was calculated as 0.172, and the scaled reserve to mass converter {ṗ AM }/ρ E was calculated as {ṗ AM }/ρ E = 0.00275 g cm −2 day −1 . The life-stage specific variables for the anchovy DEB-ABM are shown in Table III.

Food (mg m −3 ) f X mean a ω K Prey
Yolk sac larvae 0 0 -0 -0 Sinusoidal functions were applied using the equation T(t) = T mean + a sin (2π(t + ω)/365), X(t) = X mean + a sin (2π(t+ω)/365), where X mean or T mean is the mean value, a is the amplitude, and ω is the phase shift of the sinusoid curve. For larvae, prey consists of phytoplankton sized <0.2 mm, juveniles feed on small zooplankton sized 0.2-0.5 mm, and adults feed on large mesoplankton consisting of zooplankton sized >0.5 mm. A max is the given maximum age for five stages; M = 1 at the end of age-4 (1825 days).

Population traits
In a 5-year simulation, we yielded a mean population (log-transformed biomass) growth rate (r p ) estimate of 3.4 year −1 . On the basis of a comparison of individual and population-level processes, the present study found that the sensitivities of age-specific annual fecundity and population rates to short-term environmental scenarios were similar. The simulation was designed to validate the predictions of the ABM (details of the anchovy ABM variables are shown in Table IV. By simulating multiple cohorts in an ABM, inter-individual variability in growth trajectories and reproductive potential became apparent. Age class 1 anchovies accounted for 30%−60% of the spawning stock; age class 2 accounted for 20%−50%, age class 3 accounted for 10%−30%, and age class 4 accounted for less than 10%. Our model produced a pattern of biphasic growth when we varied food availability during the first stage of the anchovy life cycle. Longterm simulations were performed under different annual environmental conditions and hatching dates for length-atage and weight-at-age.

DISCUSSION
The DEB model describes vital rates at which organisms acquire and use energy for various activities, but it does not model energy allocation or growth at the population level (Sibly et al., 2013;Grossowicz et al., 2017). A combined DEB-ABM can, however, describe the dynamics of a population of individuals, where each individual follows an energy budget model (van der Meer, 2006;Kooijman, 2010). Therefore, an understanding of individual-level strategies is critical to determining the relationship between growth and time at a population level (Serpa et al., 2013;Gelman et al., 2014). The coupling of a DEB and ABM represents a promising method for prediction of population-level dynamics of species using individual data (Sibly et al., 2013).

O n l i n e F i r s t A r t i c l e
The g(P) is a function expressed by g(P) = (ΣP i v i /K)/(1 + (ΣP i v i /K)); i denotes the prey type I = 1,2,3; K is the half-saturation constant; h(T) = g Tc1 ×g Tc2 , g Tc1 = (x 1C * t * )/(1+ x 1C (t * -1)); and g Tc1 = x 4C * t # )/(1+x 4C (t # -1)), where x 1C , x 2C , x 3C , and x 4C are the values of the Cmax function corresponding to four temperatures, T 1C , T 2C , T 3C , and T 4C , respectively; x 1f , x 2f , x 3f , and x 4f are the values of the fitness function corresponding to four temperatures, respectively; . Special dynamic action (S) represents the energy allocated to the food digestive process: S = k S (C-E g ); and egestion (E G ) is a constant proportion of consumption; E g = k Eg C, where k S is the proportion of assimilated energy lost to special dynamic action, and k Eg is the scaling factors for egestion (Jin, 2004;Wang et al., 2003).
For the population-level ABM, relatively simple responses emerged from complex individual processes (Sibly et al., 2013). Population-level biomass and daily production varied with season, coinciding with numerical abundance. The length size of an individual often affects its survival, chance of reproduction, growth, and number and size of its offspring (Zuidema et al., 2010). In a 5-year simulation under standard environmental conditions, our model estimated a mean anchovy population (logtransformed biomass) growth rate (r p ) of 3.4 year −1 .
Many species have complex life cycles, with individuals at different stages classified by different attributes-e.g., plant populations have long-living seed banks, and an additional discrete-state variable is thus required to keep track of seed numbers (Ellner and Rees, 2006;Zuidema et al., 2010). After simulating multiple cohorts, inter-individual variations in reproductive potential due to inter-annual variability in environmental conditions became apparent. Intra-specific variation in population size was high during the first year and thereafter decreased with age. Using the DEB-ABM, we can begin to observe the action of population size and structure at the level of the individual. Similar observations, including population-level empirical generaliations, may be obtained once energy budget models are more routinely used in ABMs.
The DEB-ABM model explained many anchovy life history parameters well, including age-at-maturity, spawning interval, batch fecundity, seasonal fecundity, and condition dynamics. It increases our understanding of Yellow Sea anchovy reproductive strategies. Predicted feeding rates as a proportion of body weight are consistent with estimates based on stomach content (2.5%-4%) reported by Sun et al. (2006). Predicted inter-spawning intervals at peak spawning season are within ranges typically given for anchovies (3-5 days). Individual growth of anchovy results from the integration of a series of processes, such as growth, reproduction, maturation, and activity (Celisse, 2008). Formultiple-batch spawners, the level of energy reserves available for reproduction determines the number of egg batches an individual will spawn (Jager and Klok, 2010).

CONCLUSION
This modelling approach, "the coupling of an ABM with DEB theory", captured the acquisition/allocation of energy throughout the anchovy lifecycle. For a 5-year simulation, we calculated the mean population growth rate (r p ) to be 3.4 year −1 . The primary parameters [ṗ M ], [E G ], and [E M ] were estimated at 48 J cm −3 day −1 , 4,000 J cm −3 , and 2,700 J cm −3 , respectively. The population-level simulation was initialised with one cohort, consisting of one yolk-sac larva. The stage durations of individual traits are dependent on two forcing variables: SST and food density. The present study serves as a basis for further