# 48 Spatial implementation of multi-stanza and IBM

## Ecospace biomass dynamics

For functional groups not represented by multi-stanza population dynamics accounting, Ecospace represents biomass (*B*) dynamics over a set of spatial cells (*k*) with the spatially discretized rate formulation

[latex]dB_{ik}/dt=g_iQ_{ik}-Z_{ik}B_{ik}-(\sum\limits_km_{ikk'})B_{ik}+\sum\limits_{k'}m_{ik'k}B_{ik'} \tag{1}[/latex]

where *B**ik* is the biomass of functional group *i* in spatial cell *k;* *g*_{i} is conversion efficiency of food intake by group *i *into net production; *Q*_{ik} is total food consumption rate by group *i* in spatial cell *k;* *Z*_{ik} is total mortality rate of group *i* biomass due to predation, fishing, etc.; *m*_{ikk΄} is instantaneous movement rate of group *i* biomass from cell *k* to cell *k*΄; and *m*_{ik΄k} is movement rate of group *i* biomass from cell *k*΄ to cell *k*.

All of the terms on the right hand side of Eq. 1, except *g*_{i}, are treated as dynamically variable over time so as to reflect changes in food availability (*Q*_{ik}), fishing effort and predation risk (*Z*_{ik}), and seasonal changes in movement patterns (*m*_{ikk΄}). Food consumption rates *Q*_{ik} are calculated as sums over prey types *j* (i.e., *Q*_{ik} = Σ*jQ*_{jik}). Likewise, total mortality rates are calculated as sums over predator types and fishing fleets *f:* *Z*_{ik} = *M*_{oi} + Σ_{f} F_{ifk} + Σ_{j}Q_{ijk}/*B*_{ik }, where *M*_{oi} is unexplained mortality rate, the fishing rate components *F*_{ifk} by fleets *f* are predicted from spatial distributions of fishing effort for each “fleet” *f* over the grid cells *k*, and the *Q*_{ijk}/*B*_{ik} ratios represent predation rate components of *M* (i.e., *M*_{ijk} = *Q*_{ijk}/*B*_{ik}) calculated from predator *j* consumption rates *Q*_{ijk}.

The Ecospace grid cells are arranged as a rectangular grid with rows *r* and columns *c*, so that each cell *k* exchanges biomass directly only with those cells *k*΄ that are in adjacent rows and columns. If cell *k* represents row *r*, column *c*, then *k*΄ is restricted to cells (*r* – 1,*c*), (*r* + 1,*c*), (*c* – 1,*r*), and (*c* + 1,*r*). Exchanges at the map perimeter are set to zero, except for groups that are assumed to be advected across the map, in which case biomasses at the map boundary are set to constant (Ecopath base estimate) values.

**From the Ecosim Age-structured dynamics chapter**:

The basic accounting relationships for multi-stanza groups are

[latex]N_{a+1,t+1}=N_{a,t} \exp(-Z_{s,t}/12) \tag{2}[/latex]

[latex]W_{a+1,t+1}=\alpha_a q_{a,t}+ \rho W_{a,t} \tag{3}[/latex]

[latex]B_{s,t}=\sum\limits_{a=a1(s)}^{a2(s)} N_{a,t}W_{a,t} \tag{4}[/latex]

Here, *N*_{a,t} is number of age *a* (in months) animals in calendar month *t*, *W*_{a,t} is mean body weight of age *a* animals in month *t*, and *B*_{s,t} is the biomass of stanza *s*, defined as the mass (numbers × weight) of animals aged *a*1(*s*) through *a*2(*s*) months. *Z*_{s,t} is the total mortality rate of stanza *s *animals, defined the same way on the basis of fishing and consumption as for other model biomass groups *i* as *Z*_{s,t} = *M*_{os} + Σ_{f}F_{sf} + Σ_{j}Q_{sj}/*B*_{s}. All animals in stanza s are treated as having the same predation risk and vulnerability to fishing. The aggregated bioenergetics parameters *a*_{a} and *r* are calculated to make body growth follow a von Bertalanffy growth curve (with length-weight power 3.0) with user-defined metabolic parameter *K*. Exact von Bertalanffy growth occurs when predicted per-capita food intake *q*_{a,t} is equal to a base food intake rate that is calculated from the consumption per biomass parameter (*Q*_{s}/*B*_{s}) provided by the user for each stanza. The metabolic parameter *r*, which equals exp(–3*K*/12), is based on the assumption that metabolism is proportional to body weight^{[1]}.

Actual or realized food intake *q*_{s,t} at each time step is calculated from the total predicted food-intake rate for the stanza (*Q*_{s,t}) as *q*_{s,t} = *Q*_{s},t_{w}_{a,t}^{2/3}/*P** _{s,t},* where

*P*

_{s,t}is the relative total area searched for food by stanza

*s*animals and is computed as

*P*

_{s,t}= Σ

_{a}N_{a,twa,t}

^{2/3}. For foraging-arena food-intake and predation-rate calculations involving stanza

*s*,

*P*

_{s,t}is used instead of

*B*

_{s}as the predictor of total area or volume searched for food per unit time. The assumption that area searched and food intake vary as the ⅔ power of weight (i.e., as the square of body length) is a basic assumption that also underlies the derivation of the von Bertalanffy growth function. For notational simplicity, Eqs. 2-4 above are presented without a species index.

## Original representation of multi-stanzas in Ecospace

When the multi-stanza option was originally developed for Ecosim, it was not incorporated directly into Ecospace. Instead, each stanza was treated as its own higher-order functional group for Ecospace biomass-dynamics calculations without accounting for age structure within the stanza. Rather, the age-structure of each stanza was assumed to be in equilibrium. Body weight was computed grid-wide (not cell-specifically) for each stanza. Feeding rates were assumed proportional to a relative search-area index *P _{s}* calculated from a prediction of the numerical abundance of the stanza

*N*as

_{s}*P*, where

_{s}= N_{s}̅P_{s}*̅P*is the initial (t = 0) per-capita mean of the relative area-searched index

_{s}*P*

_{s,t}, i.e.,

[latex]\bar P_s=\sum\limits_aN_{a,0}W_{a,0}/\sum\limits_aN_{a,0}\tag{5}[/latex]

Dynamics of the numbers in each stanza *N _{s}* were computed for each cell by the differential equation

[latex]dN_{sk}/dt=R_{sk}-Z_{sk}N_{sk}-(\sum\limits_{k'}m_{ik'k}N_{sk})\tag{6}[/latex]

where *R _{sk}* is an approximate difference between recruitment (incoming) rates and exit (to next stanza) rates for stanza

*s*in spatial cell

*k*. If the age structure within the stanza is assumed to remain near equilibrium, the

*R*term in Eq. 6 can be approximated as

_{sk}[latex]R_{1k} = E_{tk} (1-exp(Z_{sk}(a_2(1))/12)) \; \text{for} \; s = 1 \tag{7a}[/latex]

[latex]R_{sk} = N_{s-1,k} Z_{s-1,k} / (\exp(Z_{s-1,k}(a_2(s-1)-a_1(s-1))/12)-1) \; \text{for} \; s>1 \tag{7b}[/latex]

Eq. 7a represents egg production rate minus survival rate to the age at exit from stanza *s* = 1. Egg production is assumed to be approximately proportional to biomass *B _{sk}* of the oldest (adult) stanza

*s*in cell

*k*. Eq. 7b is derived from the equilibrium of the delay-differential equation for

*N*that results from assuming spatial gain and loss rates to be approximately balanced, so that the dominant effects on

_{s}*N*are gains from individuals progressing from the previous stanza and from losses of individuals as they progress to the next stanza and mortality within the stanza.

_{s}The equilibrium assumption needed for derivation of Eq. 7 can lead to inaccurate predictions because it can result in incorrect size distributions if incoming and outgoing numbers are not in balance, and size then affects the predation-rate parameters (areas searched, maximum prey-consumption rates). Eq. 7 is a relatively poor approximation for both egg production and net rates of numbers gained through graduation from younger stanzas and loss to older stanzas, so this early version of Ecospace tended to predict incorrect absolute values for cell-specific numbers *N _{sk}* relative to Ecospace-predicted cell-specific biomasses

*B*, but the predicted spatial distributions of abundances were at least qualitatively reasonable. In past applications, Ecospace generally predicted that

_{sk}*N*was relatively high in cells with high egg production, in cells with favourable habitat, and in cells near seasonally varying optimum migration positions for migratory stanzas.

_{sk}## Predicting multi-stanza spatial distribution from continuous mixing-rate models

In an effort to avoid the large computer-memory requirements and massive accounting calculations (for typical models, on the order of 10^{3} more calculations per time step) required for replicating the full age-structure accounting for multi-stanza groups for every grid cell of large Ecospace models, a simple approach based on combining the overall multi-stanza population accounting (as described in the Ecosim Age-structured dynamics chapter, eq. 2-4 above) with the relatively simple Eq. 7 diffusion model for predicting relative spatial abundances by stanza. This approach depends on two key assumptions:

- the Ecosim multi-stanza accounting system can be applied for each multi-stanza population as a whole (summed over all Ecospace grid cells), given reasonable estimates of mean food consumption rates
*q*and mortality rates_{s,t}*Z*averaged over the grid cells (a basic assumption that is made anyway in the non-spatial Ecosim representation of any large area) and_{s,t} - the diffusion model, Eqs. 6-7, gives reasonable predictions of the relative distribution of the biomass of each stanza over grid cells whether or not the absolute numbers
*N*are predicted correctly, hence preserving effects of complex spatial-overlap patterns among stanzas._{sk}

We then perform the Ecospace time solution on monthly time steps using the following four-step procedure:

- We use the results from integration of Eqs. 6-7 to apportion the spatial distribution of total stanza biomass
*B*_{st}_{ }over spatial cells*k*to give*B*relative cell biomasses using_{sk}*B*._{sk}= B_{st}N_{sk}/Σ_{k}N_{sk} - The spatial
*B*biomasses (and relative predator-search areas_{sk}*P*) are then used in the Ecosim foraging arena and fishing rate calculations for each cell_{sk}= P_{st}N_{sk}/Σ_{k}N_{sk}*k*to predict food-consumption rates*Q*and mortality rates_{sk}*Z*._{sk} - Biomass-weighted average food-consumption rates
*̅q*, and mortality rates_{s,t}*̅Z*, for the whole population are calculated as,[latex]\bar q_{s,t}=\sum\limits_k B_{sk}q_{sk}/B_{st} \ \ \text{and} \ \ \bar Z_{s,t}= \sum\limits_k B_{sk}Z{sk}/B_{st} \tag{4a,b}[/latex]_{s,t} - The system-scale multi-stanza accounting is done by means of Eqs. 2-4 with the biomass-weighted averages (
*̅q*_{s,t}*,**Z*) to give predicted total population age and size structure and total stanza biomasses_{s,t}*B*at the start of the next month._{s,t+1}

This procedure retains some information about predicted changes in spatial abundance patterns due to mixing processes and spatial variation in mortality rates *Z _{sk}* because

*Z*is included in the prediction of relative numbers

_{sk}*N*by cell from Eqs. 6-7, but it discards information about spatial variation in growth rates

_{sk}*q*so in favour of using a single system-scale prediction of body growth (Eq. 3 with consumption rate

*q*represented by

_{s,t}*̅q*).

_{s,t}Further, it fails to account for the cumulative divergence that can take place in both age and size structure for relatively sedentary species resident in spatial cells that are protected from fishing. That is, for “adult” stanzas containing many age classes, it fails to represent the potential accumulation of older, more fecund animals in protected areas, considered by some to be a key benefit of marine protected areas. One possible solution to allowing accumulation of large adults in specific areas is to split the oldest stanza group into a number of stanzas, but doing so is an approximate fix rather than a solution.

For resident species, the mixing-model approach also fails to account for regional variation in growth rates associated with spatial cells that have higher basic (primary and lower-trophic-level) productivity or reduced intraspecific competition due to limited recruitment. For these reasons, the mixing-model approach is best suited to analyses of pelagic systems, where relatively high mobility results in averaging of feeding and mortality rates over substantial areas. When used for systems with many resident or sedentary species, the approach is potentially misleading and should be used only to provide computationally “quick-and-dirty” policy screening for options such as size and spacing of MPAs, to be followed by more careful screening according to the more detailed individual-based approach described below.

## Individual-based approach for predicting spatial patterns in growth, survival, and distribution

Most regional populations exhibit at least some degree of localized or cell-scale variation in recruitment, body growth, and survival rates, and erosion of this local structure has serious implications for maintenance of both biodiversity and overall productivity. The original and mixing-model approaches described above cannot adequately capture such local structure, which can result from the cumulative effects of the development of a fishery or from MPAs. Ecospace therefore includes a much more detailed and realistic approach to the representation of localized trophic-interaction effects based on concepts of individual-based modeling (IBM).

In the IBM approach, we retain the spatial biomass-dynamics accounting for non-multi-stanza species represented by Eq. 1 and the multi-stanza population dynamics accounting of Eqs. 2–4, but rather than solving Eqs. 2–4 once for each stanza using spatially averaged (grid-wide) food-consumption and mortality rates, we divide the age-0 recruits for each multi-stanza population (*N _{0,t}*) into a large number

*np*of packets (cohorts). Each packet is assumed to represent some number of identical individuals of the population, and all packets from the monthly recruitments start out with the same individual biomass and numbers at recruitment (

*N*). Each packet is then followed independently as it moves among spatial cells on the grid. This approach is similar to that recommended by Rose et al.

_{p,0,t}= N_{0,t}/n_{p}^{[2]}and Scheffer et al.

^{[3]}

The growth-survival Eqs. 2-3 are then solved for each packet, yielding its predicted age and size dynamics (*N _{p,a,t}* and

*W*). Packets are discarded from the overall population when they reach a maximum age (denoted amax) beyond which

_{p,a,t}*N*is negligible. Each packet

_{p,a,t}*p*is assigned an initial spatial position

*X*, and movements of the packet over time are predicted from both random (diffusive) and oriented (migratory) changes in position. At each simulation time step, the ecological conditions (food intake rates, mortality rates) for the spatial cell in which each packet is located are used in Eqs. 2–3. The overall accounting for

_{p,0,t},Y_{p,0,t}*B*and

_{s,k}*P*needed for trophic-interaction predictions (impacts from and on biomasses of non-stanza species in each cell

_{sk}*k*) then involves simply summing

*B*and

_{p,k,t}*P*packet biomasses and predation search areas over those packets present in each cell

_{p,k,t}*k*, before foraging arena predictions of

*Q*are performed for that cell.

_{sk}, Z_{sk}The obvious advantage of the IBM approach is that it retains the cumulative history of each packet’s space-use pattern, in the form of the packet’s numerical (worth) and body size (weight) states. For sedentary species, local differentiation in growth and accumulation of older animals is represented by how packets in different local areas (cells) fare over time. Further, through use of restricted movement rules, collections of packets can easily be made to form distinctive local populations, presumably key units of local adaptation and biodiversity, including as a consequence of local environmental conditions.

A disadvantage of the approach is that it requires massive computation, both for the survival-growth calculations and for movement of a sufficient number of packets over the simulated grid to permit realistic spatial distributions and variation. This number must be determined by trial and error; the number of packets must be increased until results stop changing. Most of the computational effort (typically about 90%) ends up being in the simulation of movement as changes in the locations *X _{p,a,t}, Y_{p,a,t}*.

Monthly survival-bioenergetics updates for each packet are based on food intake and mortality rates predicted for the spatial cell where the packet is located at the start of the month. No attempt is made to integrate *q* or *Z* rates over times within the month spent in different cells; doing so would be prohibitively computationally intensive. This omission amounts to assuming either that cell sizes are set large enough that most movements over any month occur within a single cell or that spatial correlation in productivity and predation risk among nearby cells there is reasonably high, so movements over such cells would result in the same predicted food intake and mortality rates obtained from the initial cell. Effects of violating this assumption could be tested if the model were run with varying grid cell sizes.

The initial or spawning position for each packet (*X _{p,0,t}, Y_{p,0,t}*) is set to the centre of a cell

*k*, where the probability of recruiting to cell

*k*is set equal to

*E*and

_{kt}/Σ_{k}E_{kt}*E*is the predicted total egg production in cell

_{kt}*k*for month

*t*summed over all packets that are in cell

*k*at the start of the month. This procedure allows spawning to occur well away from locations of larval settlement or juvenile growth because larval dispersal and juvenile migration can be explicitly represented, through either different or similar movement-simulation rules as used for packets of older fish. In particular, the IBM approach “encourages” formation of local stock structure; recruitment tends to occur near centres of egg production. In the context of MPAs, lower mortality rates Z in designated cells can result in the accumulation of older, more fecund fish, and those cells can thus become local areas of high reproduction.

Monthly movements by each packet are simulated as a set of ns increments to the *X,Y* values that determine location on the grid. The user specifies an average annual movement distance, which implies an average monthly movement distance. The number of moves *ns* each month is then set so that the distance per increment cannot exceed the width of one cell. Each movement is made only in a cardinal direction (*N,S,E,W*), so that only *X* or *Y* (not both) changes for each move. The probability of choosing each of the four directions, *k΄*, is set to *m _{skk΄}/Σ_{k΄}m_{skk΄},* where

*m*is the instantaneous movement rate from cell

_{skk΄ }*k*to

*k΄*calculated for the continuous biomass model (see Eq. 1). As noted above, the

*m*can be set equal for all

_{skk΄}*k΄*, to represent purely diffusive movement, or biased to represent avoidance of cells with unsuitable habitat, movement toward preferred habitats, or seasonal migration patterns. This method for choosing movement directions makes it possible to use the same user interface for entering assumptions about movement distances and orientation for multi-stanza populations as for groups represented only by biomasses, and it ensures that the multi-stanza movement patterns are broadly comparable with predictions from the computationally faster continuous mixing-model version of Ecospace.

**Attribution **The chapter is adapted from Walters et al. 2010, *Bulletin of Marine Science*^{[4]}*, *which permits authors to use figures, tables, and brief excerpts in scientific and educational works provided that the source is acknowledged and the use is non-commercial.

- Essington, T. E., J. F. Kitchell, and C. J. Walters. 2001. The von Bertalanffy growth function, bioenergetics, and the consumption rates of fish. Can. J. Fish. Aquat. Sci. 28: 2129–2138. https://doi.org/10.1139/f01-151 ↵
- Rose, K. A., S. W. Christensen, and D. L. DeAngelis. 1993. Individual-based modeling of populations with high mortality: a new method based on following a fixed number of model individuals. Ecol. Model. 68: 273–292. https://doi.org/10.1016/0304-3800(93)90022-K ↵
- Scheffer, M., J. M. Baveco, D. L. DeAngelis, K. A. Rose, and E. H. Van Nes. 1995. Super-individuals: a simple solution for modeling large populations on an individual basis. Ecol. Model. 80: 161–170. https://doi.org/10.1016/0304-3800(94)00055-M ↵
- Walters, C., Christensen V, Walters W, Rose K. 2010. Representation of multi-stanza life histories in Ecospace models for spatial organization of ecosystem trophic interaction patterns. Bull. Mar. Sci. 86(2):439-459 ↵