Alaska Ecosystem Models

Multispecies Fish Stock Assessment Models that Include Predation

We have developed a multispecies virtual population analysis and forecasting (MSVPA and MSFOR) for the eastern Bering Sea (Livingston and Jurado-Molina, 2000) and incorporated uncertainty estimation through the creation of a new Bayesian multispecies model (MSM; Jurado-Molina et al. 2005). These models include predation interactions among several commercially important groundfish stocks and also predation by external predators (arrowtooth flounder and northern fur seal) on these stocks. Available results include estimation of predation mortalities produced by predators on prey species, and the annual consumption of prey by predators.

Multispecies Virtual Population Analysis Models

These models include predation interactions (Fig. 1)among several commercially important groundfish stocks (walleye pollock, Theragra chalcogramma; Pacific cod, Gadus macrocephalus; Greenland turbot, Reinhardtius hippoglossoides; yellowfin sole, Pleuronectes asper; rock sole, Lepidopsetta bilineata) and also predation by external predators (arrowtooth flounder, Atheresthes stomias; northern fur seal, Callorhinus ursinus) on these stocks. The MSVPA is a recursive algorithm that calculates the fishing mortality at age Fa,t, recruitment, stock size, suitability coefficients and predation mortality based on catch-at-age data, predator ration and predator diet information. MSVPA allows the estimation of vital population rates used in the management of fishing resources. One of the advantages of this model is the estimation of the predation mortalities produced by predators on preys species and the annual consumption of prey by predators. MSVPA is an extension of the VPA model that incorporates the predator stomach content data into the virtual population. The predation mortality estimation is based on the following equations:

(1),

where M1 is the residual mortality and M2 is the predation mortality.

(2)

Where represents the average stock size of the predator i of age j. Ri , j is the annual ration of the predator and Sp, a, i, j is the suitability coefficient for each combination of predator-prey-age. This parameter reflects the diet composition of the predator relative to the available food (Sparre 1991). The denominator of equation (2) represents the total suitable biomass available to the predator. In the denominator, represents the average stock size of the prey p of age a and Wp, a ;represents the weight at age of the prey in the stomach of the predator. The third equation estimates the suitability coefficients:

(3)

where Up, a, i, j represents the fraction of prey p of age a found in the stomach of the predator i of age j.

The MSVPA input data includes the catch-at-age data (1979-2002), percent of maturity-at-age, weight-at-age, terminal fishing mortalities, predator stomach content data, over 40,000 stomach samples (Table 1)and residual mortalities.

The terminal fishing mortalities were estimated with cohort analysis and these parameters were chosen to minimize the differences between our population estimates and the ones derived from the stock assessment models used by the Alaska Fisheries Science Center and the Alaska Department of Fish and Game scientists during 1979-2002. The estimation of the residual mortalities required an initial run of the MSVPA with M1=M(M from single-species). We estimated the average M2 of the adult age classes over the model time series. This average was subtracted from the single-species M to get the residual mortality for the multispecies models.

Results from the MSVPA are deterministic and include adult stock size (Figs 2-7),average class size, 1979-2002 (Figs 8-13,age-0 recruitment (Table 2), fishing mortalities, average predation mortality (Table 3), prey biomass consumed by predators (Figs 14-19), comparison of biomass consumed by predators and yield (Figs 20-25), and comparison of total mortality z (Figs 26-31).

Multispecies Statistical Model

Predation interactions are incorporated into a statistical catch-at-age model.

Single-species models play an important role in fisheries management allowing the establishment of target reference points and setting of allowable harvest limits. Scientists involved in stock assessment have acknowledged the possibility that predator-prey interactions can influence population dynamics estimation and reference points, but the general belief is that the effects of predation and competition are subordinate to the direct effect of fishing (Sissenwine and Daan 1991). Fishery managers are increasingly asked to consider multispecies interactions in their harvesting decisions because there is an increasing tendency to recognize that fish populations are not isolated entities. Scientists acknowledge the fact that harvest limits of some prey species may depend on harvest limits of predators (Collie and Gislason 2001), thus there is a real need to link assessment models. Therefore, several attempts have been made to include biological interactions in models. The multispecies virtual population (MSVPA) and the multispecies forecasting model (MSFOR) have the potential to provide additional information to fisheries managers to improve management policies. However, their lack of statistical assumptions impedes the inclusion of uncertainty into multispecies model parameter estimation. Thus, it is important to include assumptions on process and/or observation errors in multispecies models. A new approach requires the inclusion of the separable fishing mortality assumption (Doubleday 1976; Pope 1977), which allows for common statistical estimation procedures in a multispecies context, a task that is not possible with the current available MSVPA-MSFOR technology. In the Alaska Fisheries Science Center we have developed a multispecies statistical model MSM.

In addition to the predation equations (see MSVPA), MSM incorporates the following equations:

(4)

(5)

(6)

(7)

where Na, t, q represents the number of individuals of age a in year t in the quarter q Fa, t the fishing mortality at age, Ma the natural mortality at age, Ca, t the catch-at-age, sa the age-dependent gear selectivity (u and v are its parameters),Ft the full fishing mortality.

In this model we assumed observation error in the catch-at-age and the relative indices of abundance. The log-likelihood function for each index was defined as:

(8)

where CV is the coefficient of variation and I t the observed relative index of abundance. The predicted relative index of abundance was estimated as:
(9)

The catchability coefficient was estimated using the following equation:

(10)

Similarly, the log-likelihood for the catch-at-age was assumed to be log-normal distributed:

(11)

The sum of the log-likelihood components associated with each species was used as the objective function. In the estimation process for both species a CV of 0.2 was used.

The MSM (quarterly form) for the Bering Sea include only walleye pollock (Theragra chalcogramma) and Pacific cod (Gadus macrocephalus) as predator-prey species. The MSM estimates stock size and predation mortality based on catch-at-age data 1979-2002), relative indexes of abundance (Alaska Fisheries Science Center's bottom trawl BTS and the echo integration trawl EIT surveys), predator annual ration and predator stomach contents (data from the center's food habits data base assumed to be measured without error) using estimation procedures for the statistical part and the predation mortality. The MSM statistically estimates parameters (the initial age structure in 1979 Na,0, 12 parameters; yearly age-0 recruitment Rt, 24 parameters; yearly full fishing mortality Ft, 24 parameters; and selectivity coefficients, 2 parameters) using either an optimization algorithm (e.g. Newton-Raphson) or Bayesian methods. In both cases, for a given set of parameter values, MSM projects population trajectories (for each species) over the specified time frame. These trajectories are computed based on catch-at-age (assumed known) and on predation mortality. The predation mortality requires an iterative solution since its value in each year is confounded with the abundance of other species in that year. This iterative process is referred to as the "predation algorithm" (Sparre 1991). Given fixed values for the parameters (Na,0, Rt, F>t and selectivity coefficients) and an initial guess for M2, the population trajectories are computed. These population estimates together with given values of the suitability coefficients allows the estimation of the predation mortality. The predation algorithm then adjusts the M2 values and updates population estimates until two consecutive iterations converge to marginally different M2 and suitability coefficients values according to established criteria (Sparre 1991). Once the criteria are reached and the estimates of predation mortality, suitability coefficients and population have converged, the likelihood (Equations 9-10) is used in the criteria of the main data fitting routine. This procedure (the main part and the internal estimation for M2) is repeated until the negative log-likelihood is minimized (by adjusting main parameters) or in the Bayesian methods until the posterior distribution is adequately represented. We also updated the MSVPA model for comparison purposes. Results from these models were also compared with results with the single-species statistical stock assessments carried out in the Alaska Fisheries Science Center.

Results from MSM fitted the walleye pollock and Pacific cod abundance indices well. The estimates of the suitability coefficients from both models (MSVPA and MSM) for both species were similar. The regression between both types of suitability estimates (suitability coefficients from MSVPA as the independent variable) was significant (p-value ~ 0) with a slope of 0.99 ± 0.01 and explained 99% of the variability observed.

The MSVPA and MSM (MLE) estimates of average predation mortality of walleye pollock (ages 0 and 1) followed the same trend. The estimates of Pacific cod predation mortality were not as similar as the pollock case.

Estimates of N3+ walleye pollock from the MSM, the singles-species stock assessment from the center, MSVPA (two species) and the MSVPA base-run (all predator species) were also similar.

The MSVPA model is an important tool, which can provide advice to fisheries managers on potential indirect effects of fishing in a multispecies context. However, its lack of statistical assumptions does not allow the use of tools regularly used in age-structured single-species stock assessment such as Bayesian analysis, risk analysis and decision analysis. It is likely that the predation equations can be incorporated into most of the statistical-catch-at-age models currently used in fisheries management that have the separable fishing mortality assumption (Doubleday 1976; Pope 1977).

Several improvements could be made in future versions of the MSM including the addition of more predator species. In particular, the addition of Pacific herring and rock sole (prey), which seem to be the most sensitive species to indirect effects of fishing mortality (Jurado-Molina and Livingston 2002a; Jurado-Molina and Livingston 2002b). It is also necessary to include the complete set of predator stomach content data and new data in future versions of the model. It is also important to incorporate a statistical model for the diet composition data in the future. A formal evaluation of the model through simulation of the system (adding sampling error) given a known dynamic will be necessary to see if the model is able to generate good estimates of the parameters. The final goal would be to incorporate all the data sources used in the single-species stock assessments into the multispecies statistical model. Forward simulations to test management scenarios also need to be implemented. In those simulations, future recruitment can be log-normal distributed or modeled as either a Ricker or Beverton and Holt spawner-recruit relationships to account for compensatory effects.

The main goal in the development of this model is to incorporate the predation equations into the actual stock assessment models used to provide groundfish management advice for the Bering Sea. The development of the multispecies statistical method is an important step in its adoption in providing advice to fisheries managers because it allows consideration of uncertainty in a multispecies predator-prey modeling context and it will help to establish useful scenarios for the management of groundfish resources in the Bering Sea.

Research

• Collie, J.S. and Gislason H. 2001. Biological reference points for fish stocks in a multispecies context. Can. J. Fish. Aquat. Sci.58: 2167-2176.
• Doubleday, W.G. 1976. A least squares approach to analyzing catch at age data. Res. Bull. Int. Comm. NW Atl. Fish.12:69-81.
• Jurado-Molina, J., and Livingston, P. A. 2002a. Multispecies perspectives on the Bering Sea groundfish fisheries management regime. N. Am. J. Fish. Manage.22:1164-1175.
• Jurado-Molina, J., and Livingston, P. A. 2002b. Climate forcing effects on trophically-linked groundfish populations: implications for fisheries management. Can. J. Aquat. Fish. Sci. 59: 1941-1951.
• Pope, J. G. 1977. Estimation of the fishing mortality, its precision and implications for the management of fisheries.InFisheries mathematics.Edited byJ.H. Steele. Academic Press, New York, pp. 63-76.
• Sissenwine, M. P., and Daan, N. 1991. An overview of multispecies models relevant to management of living resources. - ICES mar. Sci. Symp.193: 6-11.
• Sparre, P. 1991. Introduction to multispecies virtual population analysis. ICES mar. Sci. Symp.193:12-21.

For further information contact Kerim Aydin, Program Leader or Jesus Jurado-Molina, Author. A portion of this research was supported by the North Pacific Research Board.

Ecosystem Mass-Balance and Simulation Models

Ecosystem level models have been developed for the eastern and western Bering Sea and are under final development for the Gulf of Alaska and Aleutian Islands. Current models include separate accounting dynamics for over 120 species or species groups and include simulation of multiple gear and fishing methods and age-structured simulation of major commercial fish. Bayesian Synthesis techniques are being implemented to better estimate the relative impacts of predation, climate and fishing on marine populations, and methods of linking these simulation models to climate and lower-trophic level production are under development (Aydin et al. 2005).

Research

• Aydin, K., GA. McFarlane, JR. King, BA. Megrey, and KW. Myers. 2005. Linking oceanic food webs to coastal production and growth rates of Pacific salmon (Oncorhynchus spp.), using models on three scales. Deep-sea Res, II. 52: 757-780

Single-Species Stock Assessment Models that Include Predation

So far we have developed two of these models: one for Eastern Bering Sea pollock (Livingston and Methot 1999) and one for Gulf of Alaska pollock (Hollowed et al, 2000). We are planning on doing one for Aleutian Islands Atka mackerel in the near future. The idea of these models is that we can better examine the sources and time trends of natural mortality for these predators. We have learned that not only is natural mortality at younger ages much higher than adults but that it varies across time depending on time trends in predator stocks. This has given us better ideas of what influences predation has on fish recruitment over time and helps us to separate predation and climate related effects on recruitment. We can better show the demands of other predators such as marine mammals for a commercially fished stock and how it might influence the dynamics of that stock (although there is need for more progress in modeling the effects on marine mammals). These models will be updated periodically.

Results

Pollock cannibalism and predation on pollock by Pacific cod and northern fur seals were added to a Bering Sea pollock population model. Results indicated cannibalism was the dominant source of predation mortality, particularly of juvenile pollock and predation mortality varied on a year-to-year basis. Cannibalism appears to explain a large portion of the declining recruitment seen during periods of high adult pollock biomass (Figure 1) although the largest recruitment outliers from the estimated stock-recruitment relationship appear to be years with anomalous oceanographic conditions.

Current Model Work

The Gulf of Alaska pollock model is being updated in 2003-2004 to include two new years of groundfish food habits data and to address questions of arrowtooth flounder influence on pollock population dynamics.

Bioenergetic Models

Bioenergetic modeling of groundfish is important to quantify the relationship between fish growth, prey availability, and physical factors on the amount of food consumed by groundfish. A bioenergetic model of walleye pollock has been developed for the eastern Bering Sea (Buckley and Livingston, 1994; 1996).

Contact for More Information

• Kerim Aydin, Resource Ecology and Ecosystem Modeling Program