A nonlinear programming approach to the analysis of perturbed marine ecosystems under model parameter uncertainty

Ecological Modelling, 35 (1987) 1-28

1

Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

A NONLINEAR PROGRAMMING APPROACH TO THE ANALYSIS OF PERTURBED MARINE ECOSYSTEMS UNDER MODEL PARAMETER UNCERTAINTY C. ALLEN A T K I N S O N

System Science Applications, Inc., 121 Via Pasqual, Redondo Beach, CA 90277 (U.S.A.) (Accepted 18 November 1985)

ABSTRACT Atkinson, C.A., 1987. A nonlinear programming approach to the analysis of perturbed marine ecosystems under model parameter uncertainty. Ecol. Modelling, 35: 1-28. A general class of applied marine ecosystem problems is modeled using nonlinear programming (NLP) techniques. These problems involve multiple interacting population groups, one or more of which has been perturbed from some approximately known state via man-induced or natural phenomena. Critical ecosystem response characteristics and optimum management strategies are to be determined. The NLP modeling approach is shown to have substantial advantages for treating complex ecosystems and, in particular, for dealing with the practical issue of parameter uncertainty in the dynamic modeling. The NLP modeling approach is demonstrated by analyzing the collapse of the sardine fishery off the Pacific Coast of North America during the 1930's and 1940's. A priori predictions are shown to bound the measured sardine ecosystem response and provide evidence of overfishing as the reason for the collapse.

INTRODUCTION

In reviewing world fisheries problems, the Food and Agriculture Organization of the United Nations (1978) stressed the need for multispecies modeling approaches. Traditional single-species fish models (e.g., Beverton and Holt, 1957; Ricker, 1958) are usually inadequate since the exploited species interact strongly with predators, prey, and competitors within the ecosystem. Examples of multispecies models are described in papers by Riffenburgh (1960), Saila and Parrish (1972), Steele (1979), and May et al. (1979). The ability to estimate model parameters, which is essential to the practical application of all ecosystem models, is particularly critical in multispecies models. Data is not generally available on the temporal and spatial scale needed to estimate population interaction terms with any degree of accuracy (Goodall, 1972). Unless parameter estimates can be 0304-3800/87/$03.50

© 1987 Elsevier Science Publishers B.V.

tightly bounded, multispecies models will predict a wide spectrum of dynamic behaviour even with relatively few population groups represented (Parrish 1975; Steele 1979). The problem is exasperated as the number of ecosystem components represented in the model increases, since both the potential interaction terms and model degrees-of-freedom increase geometrically. Examples of extremely complex multi-component models are described by Anderson and Ursin (1977) and Laevastu and Favorite (1978). The present nonlinear programming (NLP) approach has been developed to deal effectively with the issue of model parameter uncertainty in applied ecosystem studies. Model parameters are treated as variables and constraint equations are formulated to relate all 6xplicit and implicit information that limits their values. By incorporating these constraint equations directly as elements of the NLP problem statement, the predictive power of the model, i.e., the resolution of the possible response modes, is greatly enhanced. The NLP solution procedures also provide parameter sensitivity results in the form of Lagrange multipliers (Luenberger, 1973). This information can be used to guide field studies to better estimate critical model parameters and to further resolve the ecosystem's dynamic character. System optimization techniques have had limited use in applied ecological problems. The most common application is optimal control theory for questions pertaining to the management of fishery resources (Palm, 1975; Vincent et al., 1975; Agnew, 1979; Ludwig, 1979). Mathematical programming techniques, which include nonlinear as well as linear and dynamic programming, have not been used as extensively (Jaquette, 1972; Jeffers, 1972; Nakamura, 1973; Mendelssohn, 1976). However, Parrish (1975) and the Oak Ridge Systems Ecology Group (1975) have suggested the potential applicability of such procedures. A general class of applied ecosystem problems are formulated as nonlinear programs in the present paper. The NLP approach is then applied to the question of the sardine fishery collapse off the West Coast of North America. The objectives are both to demonstrate the general modeling procedure and to explore the nature of the ecosystem and the reason for the sardine collapse. PROBLEM SCENARIO F O R NLP APPLICATION

The NLP approach will be applied to the problem of predicting the response characteristics of ecosystem populations. While no restrictions need be placed on the nature of the populations or the type of ecosystem, the proposed application is for marine ecosystems and higher trophic-level organisms such as fish and marine mammal species. Predictions are sought for some future perturbations to these populations or to their environment that can be expected to have significant, and typically complex interplay

between ecosystem components. Underlying the scenario is the further practical assumption that knowledge of the ecosystem's dynamic response characteristics is limited. An ecosystem dynamics model can be postulated in some general form, such as coupled differential or difference equations, but data are lacking to adequately resolve the model parameters. The above problem scenario, restated simply as the prediction of marine ecosystem response dynamics under model parameter uncertainty, is the focus for the NLP approach. The dynamics equations will be assumed by hypothesis to adequately represent the first-order features of complex population interactions and processes provided appropriate parameter values are input. The estimated range of model parameter values is assumed to be quite large, spanning both the measurement uncertainty and any inherent uncertainty associated with the ecosystem conceptualization and mathematical structuring. The NLP approach addresses the fundamental dynamic characteristics of the ecosystem and uses deterministic functions to describe them. Stochastic treatments, while offering a potentially better representation of ecosystem characteristics, further add to the estimation problem as both the nature and parameters of its probability distribution must be defined. It is assumed that the stochastic features of the modeled ecosystem can be superimposed on the NLP which determine its fundamental character. Kremer (1983) discusses various stochastic approaches used in ecosystem models. The structure of nonlinear programming is such that extreme (minimum or maximum) response characteristics will be sought. In a practical application, this information might be used to show that worst case results are acceptable or that, in fact, additional data are needed to better estimate parameters and resolve the ecosystem response. In another application, the optimum strategy for maximizing or minimizing a particular population's response might be of interest given that there is some control mechanism available. Further insight into the use of NLP procedures in applied ecosystem problems is provided by the following general development and the subsequent example in which the sardine population collapse off western North America is analyzed. NLP FORMULATION

OF APPLIED ECOSYSTEM PROBLEMS

The perturbed ecosystem problem can be stated as the following general nonlinear program: minimize • (x) subject D( x ) <~0 G(x) = 0 (1) X 0 ~< X > X c¢

I"I

PERTURBED DYNAMICS MODEL

(PARAMETER~k VALUES) J J

~.__....~ ~[

RESPONSE

OPTIMIZATION SEARCH

.l' J i

MODEL PARAMETER Ir I CONSTRAINTS F I

I

OBJECTIVE FUNCTION

Fig. 1. Flow logic showing nonlinear program elements for the perturbed ecosystem problem. Inputs to the NLP problem are represented by the ovals, while computational components are represented by the rectangles. where x is an n-dimensional vector of problem variables; ~ ( x ) is the objective function; D(x) and G(x) are vectors of implicit constraint functions that are distinguished from each other by the type of information encompassed (see below); and x ° and x ~ are vector constants defining explicit bounds on variable space. The problem reduces to a linear program if the objective and constraint functions are all linear. This will not be the case for the practical problems of interest here. The elements of the perturbed ecosystem problem are structured as an N L P in the flow diagram of Fig. 1. A description of these problem elements follows.

System dynamics model Arbitrary functions can be used to represent system response dynamics in the N L P formulation. However, for present purposes, a coupled set of difference equations will be assumed for representing the changes in ecosystem populations over time:

P(t + 1 ) = E[P(t), a, 8]

(2)

where E is a vector function relating numbers or biomass of the population vector, P, at time t + 1 to that at time t. The vector of ecological parame-

ters, a, are distinguished from the vector of perturbation terms, 8, in this representation. Components of the a vector correspond to competitor and predator interaction coefficients; to fecundity, growth and mortality rates; and to other ecological factors that quantify the response of a population to its environment and to other populations. These parameters represent variables in the context of the NLP formulation since there is a large degree of uncertainty associated with their values and, hence, in the resulting ecosystem dynamics. The components of the perturbation vector, 6, represent disturbances to the ecosystem, the effects of which are to be predicted. For example, a perturbation parameter could be defined to represent the catch rate of a new fishery to be introduced to the ecosystem. This parameter could be either estimated or treated as an additional unknown depending on the issue at hand. Before the perturbation, the population dynamics are represented by equation (2) with the vector 6 = 0. After the perturbation, 6 4:0 and the perturbed ecosystem response is determined.

Model parameter constraints Explicit constraints bounding the ecological parameters, a, can generally be stated based on ecological theory, empirical observations, or practical considerations. If the disturbance terms, 8, are treated as variables, similar bounds can be defined. These relationships are expressed as: a ° ~
8

8

(3)

Model parameters can often be compared to other parameters even though their absolute values are not well described. For example, the growth rate of a particular species, while difficult to estimate directly, may be defined as no greater than, say, one-half of another species based on physiological arguments: a 1 ~<0.5a 2

(4)

These relationships can be generalized and written in the implicit form of equation (1) as:

D(a, 8)<~O (5) where D(a, 8) is a vector of such functions. This type of parameter relationship need not necessarily be an inequality as assumed here (see application later). If there is a data base available on historic population levels, a second type of implicit constraint can be added. For simplicity, assume first that

average population levels are estimated from these data and that these levels have remained approximately constant over the time period prior to the perturbation under consideration. This simple scenario of population equilibrium implies the following set of nonlinear constraints from equation (2): i f = E[ff, a, 0]

(6)

where P is the average population value. Again this form can be generalized and represented as in equation (1) by: G(a) = P - E [ P ,

o~, 0] = 0

(7)

Note that this type of constraint does not include perturbation terms, 3. The implicit constraint functions represented by equation (7) will depend on the knowledge of the ecosystem. A varying pre-perturbation population history, if known, would yield a sequence of constraints for each time period. Stability conditions, when related to the mathematical form of the given dynamics model, can also provide implicit parameter constraints for incorporation into the NLP problem. All such constraints provide additional information for resolving ecosystem dynamics, even if the pre-perturbation populations are treated as bounded variables instead of known constants. The implicit constraints described in general by equations (5) and (7) incorporate valuable information about the ecosystem's dynamic characteristics directly into the NLP problem formulation. This is a unique feature of the NLP approach and one that provides it with much of the power in analyzing such problems. The perturbed ecosystem's response range is resolved by forcing the dynamics to behave consistent with relevant historical observations.

Objective function The objective function defines the mathematical criteria for problem solution in terms of some minimum or maximum response to the perturbation. This search for an extreme response further distinguishes the NLP approach and adds the necessary mathematical structure to a problem otherwise having many degrees-of-freedom even when the above model parameter constraints are incorporated. As a practical consideration, minimum or maximum responses are often of fundamental interest in applied ecosystem problems. The issue may be to estimate a worst case condition to be avoided (e.g., impact problems) or a best strategy to follow (e.g., fish harvesting). The result is not actually a prediction of expected response characteristics, but a determination of the potential extremes for that response. The response of the perturbed ecosystem populations can be quantified in any general way, such as individual population levels averaged over some

time period or at the end of a time period. For the example of introducing a new fishery, the response of interest might be the long-term average yield expressed in numbers, biomass, or value of the catch. Multispecies fisheries questions might best be formulated in terms of the combined catch value. Any number of response functions might be investigated depending on the ecosystem issues. Defining the response function as ~, the NLP objective is then to find the minimum (or maximum) value of ~. The response is always measured as some function of the perturbed population histories and determined via a simulation using the dynamics model (equation 2). The NLP objective can be expressed as a function of a and 8: minimize

• (a, 8)

The perturbation parameters are distinguished from the ecological parameters as before because they are not constrained by the known past history of the ecosystem, i.e., the constraints represented by equation (7). SOLUTION OF THE PERTURBED ECOSYSTEMNLP The statement of the perturbed ecosystem problem as a NLP adds mathematical structure to the already complex dynamics equations describing multispecies population interactions. The problem can easily become intractable for dynamics models involving large numbers of coupled differential equations (e.g., Laevastu and Favorte 1978). NLP solution procedures require multiple simulations in order to explore variable space, particularly in gradient search techniques where numerical gradients are continuously updated. If sophisticated integration p r o c e d u r e s (e.g, Runge-Kutta, Adams-Moulton) are needed to perform these simulations, the problem can get rapidly out of hand. One way to improve the basic computational efficiency for complex ecosystem models is to use a difference form for the dynamics equations, provided of course such equations are representative of the ecosystem response characteristics. A general difference model is presented by Atkinson (1986) for representing multispecies fish population dynamics. Essential population features and processes described in his model include spawning, growth, mortalities, age classes, temporal variability, and migration. Atkinson (1980) has considered various NLP solution techniques and has applied two state-of-the-art procedures to the perturbed ecosystem problem: (1) a generalized reduced gradient (GRG2) routine developed by Lasdon et al. (1978), and (2) a sequential simplex search (S 3) routine developed by Atkinson. The capabilities of both software routines were demonstrated for this class of problems and have been used in the analysis described below.

ANALYSIS OF THE SARDINE ECOSYSTEM OFF THE WEST COAST OF NORTH AMERICA

Sardine and 'sardine-like' fishes are found throughout temperate and subtropical coastal waters of the world's oceans, and represent a major fishing resource that has long been under intensive exploitation. These populations have also been found to be among the most vulnerable to fishing pressure as illustrated by the dramatic collapse of the Pacific sardine population off the West Coast of North America during the 1930's and 1940's (Murphy, 1966). The latter situation has been selected for the modeling application because of the sardines' general importance as a world fishery resource as well as the relative availability of applicable data. The problem to be addressed by the NLP approach is that of making a priori predictions of the sardine's rapid decline. The perturbation is assumed to be the large increase in fishing pressure that occurred concurrent with the population collapse. The objective will be to predict bounds (minimum and maximum) on the potential sardine response given only the historic population data prior to the time of increased fishing pressure. Comparisons will be made between these a priori predictions and the actual (estimated) population response.

General ecosystem description The waters of the California Current System flow southward along the Pacific Coast of North America, covering a general region from near the U.S.-Canadian border southward to the tip of Baja California and offshore about 200 n miles (Fig. 2). These waters, which support relatively high levels of biological productivity, can be viewed as an ecological system (Sette, 1969). The sardine represented the largest fish population of the California Current region during the early years of this century. The sardine fishing fleet, which began to expand following World War I, had yearly catches as high as 200 kt (in 1930) while the population remained at near-virgin stock levels. However, beginning in 1932, both the fishing effort and the catch increased significantly, the catch reaching a peak of almost 800 kt by 1936 (Fig. 3). The period from 1932 to 1936 also coincided with the beginning of the major decline of the sardine population as shown by the estimates in Fig. 4. The fishing effort remained roughly constant for the next 10 years or so kt, kiloton = 103 (metric) tonne = 106 kg. n mile, (international) nautical mile = 1.852 km (def).

BRITISH COLUMBIA 50°

45° OREGON ~"

PACIFIC NORTHWEST N. CALIFORNIA S. CALIFORNIA BAJA CALIFORNIA

40°

35°

30°

AVERAGE YEARLY SARDINE CATCH ( 1 9 2 0 - 1950)* = = = =

44K TONS 208K TONS 115K TONS 0 TONS**

- -

- -

25° --

D

SARDINE DISTRIBUTION

[]

FISHING LOCALITI ES

"FROM DATA IN MURPHYI l l S ,

\

**BAJA FISHERY NOT SIGNIFICANT UNTIL AFTER 1950

I 125 °

I 120 °

I 115 °

I 110 °

Fig. 2. Map of the California Current region showing sardine distribution and major fishing localities in the period before 1950 (from Murphy, 1966).

following the peak catch of 1936. The sardine population rebounded sharply during this time due primarily to two large year classes in 1938 and 1939, before a further dramatic reduction occurred in the 1940's. Subsequent decreases in the population have occurred and it is now only a minor component of the ecosystem. A striking feature associated with the sardines' decline is the increase of the anchovy population within the California Current region (Murphy, 1966). Since anchovy represent the major competitor to the sardine with

10

8011-

CATCH

600--

1920

30

40

50

1960

YEAR

Fig. 3. Sardine catch and effort statistics during the major sardine collapse period (from Murphy, 1966). similar food requirements and over-lapping boundaries of their natural habitat, the common explanation for this decrease is that the anchovy replaced the sardine within the trophic structure (Murphy, 1966; Gulland 1971). The fishery for anchovy was quite small during the peak years of the sardine. It consisted of a commercial fleet and a live-bait fleet whose total landings did not exceed 100 kt until just recently (Huppert et al., 1980). Prior to 1950, the total was always below 10 kt with the exception of 1947. Competitors to the sardine and anchovy consist of diverse fish species including jack mackerel and Pacific mackerel as well as numerous invertebrate species which dominate the biomass of the group as a whole (Riffenburgh, 1969). Predators to the sardine and anchovy consist of generally larger, slower growing fish and include such important species as hake and barracuda (Murphy, 1966). Many of the fish populations within the competitor and predator groups have been subjects of intensive sport and commercial fisheries during the period of interest in the 1930's and 1940's. For example, the commercial mackerel fishery flourished during this time, reaching a peak of 73 kt in 1935, before declining rapidly (Parrish and MacCall, 1978). At the base of the food web, the sardine trophic levels is supported by the abundant phytoplankton and zooplankton populations of the California

11

/

°

I

/', I

~)

SARDINE

o 1920

[

I

I

30

40

50

1960

YEAR

Fig. 4. Estimated adult sardine and anchovy populations during the major sardine collapse period (from Murphy, 1966).

Current. Young sardine and anchovy stages feed directly on the phytoplankton while the adults feed primarily on the zooplankton (Huppert et al., 1980). In fact, Murphy (1966) using caloric intake arguments indicates that the sardine was consuming about one-half the total zooplankton production in its peak abundance years. Zooplankton concentrations show significant yearly fluctuations as do the phytoplankton resulting from highly variable environmental influenzes including coastal upwelling and offshore transport. These environmental effects, particularly as they impact the survival of the young vulnerable sardine and anchovy stages, can cause large fluctuations in year class strength (Lasker, 1978; Parrish et al., 1981; Methot, 1983).

Ecosystem dynamics model The present modeling effort focuses on the sardine and anchovy populations and considers the ecological subsystem defined by Riffenburgh (1969) and illustrated in Fig. 5. This obvious simplication of the complex trophic food web of the California Current is intended only to be representative of the important sub-classes of the sardine and anchovy and the aggregated populations that significantly interact with them at various stages in their

12

I I ®,REO,OR" ANc.ovY ® A,c.ow .

k."

J

Fig. 5. Schematic showing interactions between sardine ecosystemgroups as modeled by Riffenburgh(1969). Competitiverelationships are indicatedby the connectinglines with dual arrowheads, while predator-prey relationships are defined by arrows pointing to the predator. development. Details of the ecosystem model assumptions are provided by Riffenburgh (1969). The sardine is divided into three age groups: larvae, yearlings, and adult. The larvae group encompasses all development stages during the first year of life when the sardine is the most vulnerable to environmental and predatory influences. The yearling stage likewise lasts a year, with adulthood consisting of the sardines 2 years old and above. The adult group contains all the sexually mature members of the population. The anchovy population is divided into two age classes, larvae and adults. The competitor and predator groups are assemblages of a wide variety of species which are lumped together for purposes of representing their composite interactions with the sardine and anchovy. The predators are assumed to have additional food sources outside the California Current region a n d / o r ones whose dynamics are not closely coupled to the other ecosystem components. Phytoplankton and zooplankton are represented by implicit carrying capacity terms as to their effect as a sardine trophic-level food source.

Difference equations The biomass dynamics of the seven population groups in Fig. 5 will be described by a set of coupled discrete-time difference equations. The assumed time step is 1 year so that only year-to-year population changes are represented. The simulation period will span the years from 1932 to 1952

13 when the northern sardine subpopulation of the California Current region dominanted the southern subpopulation found off Southern California and Baja California (Murphy, 1966). Regional differences, while ignored for the 1932-1952 time frame, would be important after the northern subpopulation collapsed. The general form of the difference equations is presented below and was taken from a multispecies fish population model presented by Atkinson (1986). Note that subpopulations are defined for each age class of the sardine and anchovy:

Pi(t + 1 ) =

Y'~ F,.j Sj(t)

exp(-ajP(t)) exp(-fljP(t))

jcA(i)

x[1 - Rj exp(-~,j

P(t))]Pj.(t)

where Pi(t+ 1)=biomass of age class subpopulation i at time t + 1; Pj.(t) = biomass of age class subpopulation j at time t; A(i) = set of all age class subpopulations that "communicate" with subpopulation i via fecundity, or aging processes (see below); F, j = fecundity/aging function (see below); Sj(t)= maximum survival/growth rate of subpopulation j during time period t to t + 1; aj = competition vector relating competitor effects on the survival/growth of subpopulation j; flj = predation vector relating predator effects on the survival/growth of subpopulation j; Vj= prey vector relating prey availability effects on the survival/growth of subpopulation j; Rj = starvation mortality rate under conditions of no (explicitly modeled) prey. The set A (i) includes age class subpopulations of the same species group. The fecundity/aging function F,j then relates the biomass transfers between these age classes as follows: {! F,j =

whenjisthematureadultclassandiis the first year (larvae) class when j is preceeding year class to i or j = i is the adult class otherwise

where f is the adult fecundity expressed as the yearly spawned biomass per unit adult biomass. F~j = 1 for the single age classes of the competitor and predator groups with spawning effects included in the survival/growth term. The time-step difference equations representing the coupled dynamics of the seven California Current subsystem populations are given in Table 1. The equations describe the yearly changes in population biomasses as a function of the important first-order processes effecting their dynamics.

14 TABLE 1 Difference equations describing biomass dynamics of the sardine subsystem populations of the California Current region during the period 1932-1952 Population 1 - - s a r d i n e larvae P l ( t + 1) = ):3 S3(t) e x p ( -

a33P 3 -

a35P 5 - ot36P6) exp(

- •37P7)

P3(t)

Population 2 - - sardine yearling P2(t + 1) = S 1 exp( - a l l P 1 -- o~14P4) exp(-/~13P3

-

fllsP5 --/~16P6) P l ( t )

Population 3 - - sardine adult P3(t + 1) = $2 exp( - ot22P2 - - ot25Ps) e x p ( -

1~26P6)

P2(t)

+ S3(t ) e x p ( - ot33P3 - ot35P5 - ot36P6) e x p ( - f157PT) P3(t) Population 4 -

anchovy larvae

P4(t + 1 ) =f5S5 e x p ( - ot53P3 - ot55P5 - ot56P6) e x p ( - f157P7) P5(t) Population 5 - - anchovy adult Ps(t + 1) =

$4 e x p ( - o~41P1 -- ot44P4) e x p ( -

f143P3 - / ~ 4 5 P5 - / ~ 4 6 P 6 )

e4(t)

+ $5 e x p ( - a53P3 - a s s P 5 - a56P6) e x p ( - flsvP7) Ps(t) Population 6 - - competitor group P6(t + 1 ) = $6 e x p ( - a63P3 - a65P5 - a66P6) e x p ( - f167P7) P6(t) Population 7 - - predator group Pv(t + 1) = $7 e x p ( - avTP7) (1 - R 7 e x p ( - ]173P3 - ~ 7 5 P 5 - " / 7 6 P 6 ) ) PT(t)

Age class sub-populations are defined for the sardine (1-3) and the anchovy (4-5). All population variables in exponentials are assumed to be at time t.

The ecological parameters, except for the adult sardine survival term S3(t ), are assumed to be independent of temporal and spatial features during the period of interest from 1932 to 1952. As mentioned previously, explicit stochastic effects are neglected in the N L P approach, but can be superimposed on the more fundamental deterministic response modes under investigation here. Larvae prey terms are ignored in the equations for adult sardine, anchovy, and competitor feeding since the larvae are a secondary food source compared to zooplankton. The S3(t) term is expanded below to explicitly define the fishing rate, 83(t), which is evaluated later as a perturbation function. S3(t ) = $3.0 [1 - 63(t)] (9) $3,0 is the adult sardine maximum survival/growth rate exclusive of fishing mortalities. Fishing rates are included in the other survival/growth terms of the model and are assumed to be either constant a n d / o r minor mortality components.

15 The simplicity of the above mathematical formulation does not allow for representation of all ecological processes and effects operating between populations and their environment. Hopefully, dominant mechanisms are represented for the problem at hand, and when the estimated range of parameter values is substituted into the equations, the simulations span the potential ecosystem response characteristics.

Ecological parameter estimates The dynamics model described in Table 1 consists of a large number of ecological parameters, most of which can not be estimated with any degree of certainty. This parameter uncertainty reflects the general scenario developed previously for applying the nonlinear programming approach to perturbed ecosystem problems. The ecological parameters used in the sardine ecosystem model are categorized for purposes of estimation as: (a) parameters that are (approximately) known constants; (b) parameters that are treated as independent variables; (c) parameters that are treated as dependent variables. There are only two assumed knowns, maximum adult sardine and anchovy survival/growth rates ($3,0 and $5) both of which can be reasonably estimated from the extensive data in the literature (Clark and Phillips, 1952; Murphy, 1967; MacCall, 1970; Huppert et al., 1980). The estimates are $3,o = 1.40 for sardine (exclusive of fishing mortalities) and S 5 = 1.20 for anchovy (including fishing mortalities). The independent variable category consists of those parameters thought to be the most critical to the dynamic response characteristics of the ecosystem, but with limited available data to estimate their values. The independent variables are defined in Table 2 and estimated bounds on their values are presented. Two carrying capacity terms are included in Table 2, one for the sardine trophic level (K 2) and one for the predator trophic level (K3). Estimates by Riffenburgh (1969) and Murphy (1966) indicate that a sardine trophic-level biomass of the order of 10000 kilotons could be supported by the average zooplankton crop under ideal conditions (no predatory losses). The sardine trophic level could then, in turn, support of the order of 1000 kt at the predator level assuming atrophic efficiency of 10%. These data are reflected in the estimates in Table 2. The range of values for the survival/growth parameters ($1, $2, $4, $6, and $7) were based on relative estimates given in Riffenburgh (1969) as compared with adult sardine survival/growth. The estimated uncertainty is greatest for the larvae stages of sardine and anchovy because of both the difficulty in making survival/growth measurements and the inherent varia-

16 TABLE 2 Estimated bounds on the ecological parameters treated as independent variables Population or trophic level

Parameter type

Symbol

Parameter range

Sardine trophic level

Carrying capacity

K2

1 x 10 4_4 × 10 4

Predator trophic level

Carrying capacity

K3

1 X 10 3-4 × 10 3

1-- Sardine larvae

Survival/growth Competition Predation

S1 an fl13

1.0-7.4 0.1/K 2- 1 . 0 / K 2 10 4_10- 3

2-- Sardine yearling

Survival/growth

$2

1.5-2.2

3--Sardine adult

Competition Predation Fecundity

a35 fl37

f3

0-0.3/K 2 0-0.6/K3 0.2-0.6

4--Anchovy larvae

Survival/growth Competition Predation

$4 a44 fl45

0.5-2.7 0.1/K 2-1.0/K 2 10 -- 4--10- 3

5--Anchovy adult

Competition Predation Fecundity

~53 f5

0-0.6/K2 0-0-6/K3 0.2-0.6

6--Competitor group

Survival/growth Predation

$6 fl67

1.3-2.0 0-0.6/K3

7-- Predator group

Survival/growth Mortality Prey

$7 R7 '),73

1.05-1.40 0.1-0.6 0-5.0/K 2

fl57

Carrying capacity (K) is in units of kt of biomass. Survival/growth (S), fecundity (f), and interaction parameters (a, t , and 7) are specific rates per year, while mortality (R) is a fractional value per year.

bility in the e n v i r o n m e n t a l influences. T h e c o m p e t i t o r biomass, w h i c h is d o m i n a t e d b y i n v e r t e b r a t e s with relatively short life cycles ( R i f f e n b u r g h , 1969), is e x p e c t e d to increase m o r e rapidly t h a n the a n c h o v y to fill a n y void left b y the decreasing sardine p o p u l a t i o n . Also, the relative p r e d a t o r s u r v i v a l / g r o w t h estimates indicate that this g r o u p r e s p o n d s m o r e slowly to c h a n g e t h a n either the a n c h o v y or c o m p e t i t o r group. T h e c o m p e t i t i o n terms in T a b l e 2 are a s s u m e d to be d o m i n a t e d b y f o o d s u p p l y limitations a n d can t h e r e f o r e b e r e f e r e n c e d to trophic-level c a r r y i n g capacity. F o r the classic case of self c o m p e t i t i o n , the V e r h u l s t - P e a r l logistic e q u a t i o n (Pielou, 1977) relates the c o m p e t i t i o n t e r m to g r o w t h rate (equal to

17 the natural logarithm of S in our difference equation terminology) and carrying capacity as follows: In S ct- K (10) The range for self-competition terms (Oql and ~44) was selected to correspond to the estimated uncertainty in the survival/growth terms. Competition parameters between populations (~3s and as3 ) were allowed to vary over a range from zero to the maximum self-competition values. The range of values for the predation terms effecting sardine trophic-level populations (fl37, fl57, and /367) were estimated from the applicable data given in Murphy (1966) and Riffenburgh (1969). The predation terms for sardine and anchovy adults feeding on their own larval stages (fl13 and fl45) were derived from Gulland's (1971) correlation of sardine larvae survival rate versus total adult biomass. Although the relative contribution of adult predation versus larval competition cannot be distinguished from these data, the predatory effect was bounded as between 10 -3 and 10 4 for both sardine and anchovy larvae. All predation terms in Table 2 were expected to have major implications on the perturbed dynamics as populations and mixes of populations change. The prey parameter in the predator group's dynamic equation (~'73) relates the value of sardine as a food source, while the starvation mortality parameter (R7) defines the general importance of all sardine trophic-level populations. With such a diverse predator group and given little meaningful data to estimate these effects (e.g., Riffenburgh, 1969), values ranging from one indicating little or no influence to one that reflects a significant dependence were established. The remaining two ecological variables in Table 2 are the sardine and anchovy fecundity parameters (f3 and fs). The spawned biomass is estimated to be between 0.2 and 0.6 of the adult spawning biomass based on Riffenburgh's (1969) estimate of 0.36, which was derived from estimates of the weight and number of eggs laid by the female per season. Sardine and anchovy should both have approximately the same fecundity ratio since the anchovy lays twice as many eggs but they are only one-half the weight. The parameters classified as dependent variables are defined in Table 3 and expressed in terms of the independent variables. The dependent variable category contains those parameters that could be reasonably defined as functions of other parameters a n d / o r that appeared to be less critical to the ecosystem response dynamics. These relationships represent an equality form of the implicit parameter constraints defined previously by equation (5). Perhaps the most significant of the dependent variable relationships are those given for the competitor group. These relationships define relatively

18 TABLE 3 Functional relationships for the ecological parameters treated as dependent variables Population

Parameter type

Symbol

Functional relationship 0.5 a 11 0"5•13

1-- Sardine larvae

Competition Predation Predation

o~14 ill5 fl16

2-- Sardine yearling

Competition Competition Predation

a 22 a25 ~26

In S 2 / K 2 0.5 In S2//K 2 0.5 In $2//K2

3-- Sardine adult

Competition Competition

a33 0~36

0.3/K 2

4--Anchovy larvae

Competition Predation Predation

a41 fl43 /346

0.50/44 0"5fl45 0.1fl45

5--Anchovy adult

Competition Competition

a55 a56

0.6/K 2

6-- Comptitor group

Competition Competition Competition

o/63 a65 ot66

2.0 In $6/K2 2.0 In S6/K 2 2.0 In $6/K2

7--Predator g r o u p

Competition Prey Prey

a77 3'75 3'76

0.5 In $7//K3 3'73 0.5 3'73

0"1fl13

0.1/K2

0.2//K2

large competitive exponents compared to the analogous sardine and anchovy terms, indicating that the competitors will be ineffective competitors against them in the long run (Murphy, 1966; Riffenburgh, 1969). However, the competitors will tend to respond more rapidly to any void left b y decreasing sardine or anchovy populations in the ecosystem based on previous survival/growth rate estimates.

Sardine fishing perturbation The dramatic decline of the sardine population during the 20-year period between 1932 and 1952 has generally been associated with the large concurrent increase in fishing pressure. This fishing pressure increase is defined as the perturbation function in the present N L P problem scenario. Although environmentally related fluctuations in year class recruitment success have also been postulated as the cause of the sardine fishery collapse, Atkinson (1980) indicates that this would not, in fact, have happened without the increased fishing pressure. Certaintly, the sardine's decline would not have occurred in such a short time period.

19

1.0

ESTIMATED FISHING RATE DERIVED FROM MURPHY'S (1966) DATA

?. I-

. . . . . . -.~

~< 0.5 0

/ \ - - . _ . _ . / ~ _.^. = \ ,,] i

"'--''/

I<

~ 0

/

1930

"

FISHING RATE MODEL

I

L

I

I

I

34

38

42

46

50

YEAR

Fig. 6. Model of sardine fishing rate, 83(t ), used as the perturbation function in the NLP problem formulation.

The actual (estimated) sardine fishing rate from 1932 to 1952 is presented in Fig. 6 in terms of 83( 0, which is defined as the fraction of the total adult population caught during a fishing year. These data were derived from the estimates presented previously in Figs. 3 and 4. The postulated or modeled fishing rate is also shown in Fig. 6 and reflects the general trend of the actual (estimated) increase that occurred during the 1932-1952 time frame. The postulated rate increases about four-fold from a fishing rate of 63 --0.10 before 1932 to a rate of ~3 ~" 0.40 by 1937. This large fractional catch rate is assumed to remain constant over the remainder of the period to 1952, Initial conditions The sardine population before 1932 were at high, perhaps virgin levels of about 4000 kt based on estimates shown in Fig. 4. The anchovy population is estimated to be about 1000 kt during this same period. Direct estimates of competitor or predator biomass are not available for the period before or after 1932, but approximate levels can be derived based on trophic level values. If the total biomass at a trophic level remains relatively constant independent of species composition changes as generally proposed (Sette, 1969; Steele, 1979), then the competitor group biomass can be found by subtracting the sardine and anchovy biomasses from the total. This gives an approximate value of 3000 kt for the pre-1932 competitor group based on Riffenburgh's (1969) trophic-level estimate of 8000 kt. Riffenburgh's trophic level estimate for predators was of the order of 2000 kt. Estimates for the sardine and anchovy larvae and yearling stages are also available from Riffenburgh (1969). The estimates of pre-1932 biomasses are summarized below. These represent the biomasses at the start of 1932 which is defined to begin in the

20 summer just after the main spring spawning season of the sardine and anchovy: sardine larvae

F 1 = 1600 kt

sardine yearling

fiR = 300 kt

sardine adult

53 = 4000 kt

anchovy larvae

54 = 400 kt

anchovy adult

~ = 1000 kt

competitor group

fi6 = 3000 kt

predator group

57 = 2000 kt

The initial population biomasses are treated as known constants even though the limited available data does not allow their accurate determination. However, the estimates are more accurate than the majority of the ecological parameter estimates in Tables 2 and 3. Furthermore, it is convenient to maintain fixed reference values for making relative comparisons after the ecosystem perturbation effects the populations.

Implicit population constraint functions The ecosystem dynamics model introduced above allows for a wide range of dynamic responses to the given perturbation since most of the parameters are ill-defined and treated as variables with large uncertainty bounds. In fact, if no further condition is placed on the value of these parameters, almost any response could be predicted due to the increased fishing rates of the 1930's and 1940's. For example, a simulation case was run using certain selected parameter values within the range of Table 2 and resulted in the sardine population actually increasing. While logic alone would tend to discount this behavior, additional mathematical structure needs to be introduced in order to force the model to behave 'properly'. This takes the form of implicit population constraints associated with the pre-perturbation conditions. The ecosystem populations are assumed to be in an approximate equilibrium state prior to the fishing perturbation in 1932 based on the limited but consistently high sardine population estimates. Functional constraints are defined by substituting the known (estimated) initial populations into both sides of the seven model equations presented in Table 1:

PI(t+I)=PI(t)=P 1 P2(I+ I ) = P 2 ( I ) = P 2 "

P (t+

m

(11)

21 By forcing the ecosystem dynamics model to reflect an equilibrium state, the parameters are coupled in a significant fashion. Future (post-1932) response characteristics of the ecosystem are limited to those perturbed modes which are also consistent with the equilibrium state. Stability conditions are ignored here since it is not clear whether the sardine-anchovy subsystem is naturally stable. For example, Soutar's (1966) data on fish scales deposited in coastal sediments over the past 1000 years indicate that the sardine has undergone large changes in abundance even when there was essentially no fishing pressure. The implication is that natural causes, perhaps effecting year class spawning success, can swing what may be a tenuous population balance between the sardine and anchovy.

Objective function The response of the ecosystem to the post-1932 fishing perturbation will be measured in terms of the adult sardine population. The interest is to predict the longer term effects of the fishing increase and to bound the resultant sardine population change. The response will be quantified as the adult sardine population average over the final 15 years of the 20-year simulation period after the assumed sardine fishing rate has reached the maximum of 83 -- 0.40. This response function is defined by: 1952

=~

E

P(t)

(12)

t = 1937

with P(t) determined from a simulation using the dynamics model in Table 1 with the perturbation function defined in Fig. 6. Two objective functions will be evaluated to bound the potential population response: (1) lower bound: minimize (2) upper bound: maximize This completes the NLP problem description. The computed results will predict the extremes of the sardine population change due to increased fishing pressure during 1932-1952. The large parameter uncertainty in the dynamics model will be constrained by the assumed pre-perturbation equilibrium condition.

Computational results and analyses Computations were performed to solve the two NLP ecosystem problems corresponding to minimum and maximum objective cases. FORTRAN computer codes implementing the GRG2 procedure (Lasdon et al., 1978) and the S 3 procedure (Atkinson, 1980) were both used in the analysis. The

22

MINIMUM SARDINE RESPONSE -- MAXIMUM SARDINE RESPONSE

4000

~

~

ESTIMATED FROM SARDINE DATA BY MURPHY (1966)

3OOO

f

.p, ,ANCHOVY

2000

1000

ESATEDROJ

ANCHOVY DATA BY MUR.PHY (~1966) ,

o

I 1932

I

I

I

36

,

~ ' ~ - -. ,, ~ ' " ~

I

40

I

I

44

I 48

52

YEAR (a) SARDINES AND ANCHOVIES

8000 m ~--

MINIMUM SARDINE RESPONSE MAXIMUM SARDINE RESPONSE

6OOO D

COMPETITOR v <

4000

O PREDATOR GROUP 2000

0

L 1932

I 36

t

I

I

40

I 44

I

I

I

48

52

YEAR

(b) COMPETITOR AND PREDATOR GROUPS

Fig. 7. Minimum and maximum sardine response calculations for the period 1932-1952 showing adult sardine and anchovy population trajectories and competitor and predator population trajectories.

calculated results of the two codes were found to be in close agreement and it was concluded that near-optimal results had been found for the minimum and maximum bounding cases. The S 3 calculations are presented and described below. The predicted minimum and maximum sardine response curves are compared in Fig. 7a with the actual response estimated by Murphy (1966) during the 1932-1952 time period. Discounting fluctuations which can be related to yearly larval survival rates (Atkinson, 1980), the predictions appear to bound both anchovy and sardine data very well. The sardine

23 population is predicted to be less than 1000 kt by 1952 and to continue to decline unless the fishing rate or other assumed conditions changed. The anchovy is predicted to increase to between 2000 and 3500 kt by 1952, while the competitor group increases at a somewhat faster rate to over 4000 kt (Fig. 7b). Both fill the ecological void left by the declining sardine. Meanwhile, the predator group is predicted to decrease slightly as its food supply is less desirable, i.e., components of the competitor group have replaced part of the sardine population. Comparing the predicted minimum and maximum sardine trajectories with each other indicates that they exhibit similar characteristics and closely bound the actual (estimated) population history. One reason is, perhaps, that unrealistic limitations were assumed for the range of model parameters. The knowledge of the ecosystem in 1932, when such predictions would be of most interest, would probably not have allowed the tight bounding of the parameter values or the authoritative description of the ecosystem structure that present hindsight allows. On the other hand, the assumed pre-perturbation population equilibrium constraints were shown by Atkinson (1980) to strongly limit the possible model response characteristics. Widely different behavior resulted when these implicit constraints were removed. The potential advantages of the NLP formulation for these situations is evident. The variables and objective function for the minimum and maximum response cases are presented in Table 4. An analysis of these data indicates that the maximum sardine response case results when sardine larvae predation by sardine adults (parameter /~13) is relatively high and the sardine larvae survival/growth rate ($1) is near the assumed upper bound. The resilience of the young sardine stages are improved because of the relatively strong density dependence and intrinsic growth capability. A significant number of young sardines can survive even as the adult population decreases; the population as a whole is less effected by increased fishing pressure. When density dependence is low and with sardine larvae survival rate at the lower bound, the opposite occurs and the minimum bound conditions are defined. Parameters associated with the competitor and predator groups appear to have only secondary effects and the trajectories of these groups are about the same for both the minimum and maximum cases (Fig. 7b). Parameter sensitivity information can be obtained from the Lagrange multipliers which were computed and output by the NLP solution procedures. Normalized multipliers are defined at optimality (minimum or maxim u m condition) by: ~ i - ~(~i/~i)

(13)

24 TABLE 4 Minimum and maximum sardine population calculations: variables and objective Optimization results

Minimum case

Maximum case

Sardine trophic level carry cap., K 2 Predator trophic level carry cap., K 3

1.671 × 104 1.824 × 103

1.996 X 104 1.973 X 103

S. larvae survival/growth, $1 S. larvae competition, all S. larvae predation, 1313

1.000

O.IO00/K2 3.483;410 .4

7.389 0.1000/K 2 7.648×10 4

S. yearling survival/growth, $2

1.711

2.200

S. adult competition, a35 S. adult predation,/337 S. adult fecundity, f3

0.2211/K 2 0.1364/K 3 0.4542

0.1953/K 2 0.2040/K 3 0.4715

A. larvae survival/growth, S 4 A. larvae competition, 0t44 A. larvae predation,/345

0.637 0.1000/K 2 1.546 × 10 -4

0.500 0.1000/K 2 3.878 × 10-4

A. adult competition, a53 A. adult predation, /357 A. adult fecundity, f5

0.5034/K2 0.2175/K 3 0.4747

0.3828/K2 0.1853/K 3 0.4255

Competitor survival/growth, $6 Competitor predation,/367

1.609 0.0132/K3

1.626 0.0950/K3

Predator survival/growth, $7 predator mortality, R 7 Predator prey, Y73

1.224 0.4523

1.229 0.4906

4.257/K2

5.000/K2

537.54

1596.78

Variables

Objective function

Objective is averaged sardine adult biomass (see text for definition); variable units are defined in Table 2.

where ~ri is the normalized Lagrange multiplier associated with a tight (equality) bound on parameter ai; ~* is the optimal objective value; and if, is a reference value. The reference values are used so that the magnitude of the normalized multipliers provide a direct relative measure of each parameter's importance. A summary of the normalized Lagrange multipliers are given in Table 5 for the minimum and maximum response cases. Multipliers are defined only for those ecological parameters at their lower or upper bounds. If the bounds on the other parameters are changed such that they become tight, then associated multipliers would be calculated and sensitivity information would be obtained.

25 TABLE 5 Minimum and maximum population calculations: Lagrange multipliers Case

Population

Parameter type

Symbol

Bound

Min

Sardine larvae

Survival/growth Competition Competition

$1 all ~44

Lower Lower Lower

Survival/growth Competition Survival/growth Survival/growth Competition Prey

$1 all $2 $4 a44 "/73

Upper Lower Upper Lower Lower Upper

Anchovy larvae Max

Sardine larvae Sardine yearling Sardine adult Predator

Normalized Lagrange multiplier - 672.6 - 19.2 - 65.7 -

839.0 243.0 151.8 244.6 - 90.8 - 492.3

The largest Lagrange multipliers in Table 5 are associated with the sardine larvae survival/growth term ($1) for both minimum and maximum cases. The other relatively significant value corresponds to the sardine prey term in the predator equation (~73) for the maximum case only. The competition terms seem to be less important except for the sardine larvae self-competition term (all) in the maximum response case. Since lower b o u n d sardine predictions would perhaps be of most interest, the sensitivity results point to better resolving the sardine larval year s u r v i v a l / g r o w t h parameter in order to improve predictive capability. H o w ever, this factor has obvious stochastic features that make it difficult to determine and perhaps inappropriate to represent as an average deterministic value in the dynamics equations. A second generation analysis would be to develop a stochastic representation. CONCLUSIONS A general class of applied ecosystem problems was modeled as a nonlinear program. Applications include problems where a given ecosystem is perturbed from some approximately known state and the minimum or maximum response characteristics are to be predicted or o p t i m u m strategies are to be determined. The N L P approach can be used to effectively treat uncertainty in the dynamics model parameters and to resolve response modes to the degree possible using all implicit information available from historical ecosystem records. The special N L P solution techniques also provide valuable sensitivity information (Lagrange multipliers) as to the critical ecological factors causing the predicted response.

26 The N L P modeling approach was demonstrated via the analysis of the sardine ecosystem off western North America. The calculations showed that relatively complex ecosystem problems can be formulated and solved used available state-of-the-art N L P software packages. Predictions tightly b o u n d e d the observed ecosystem response for the period of sardine population collapse during the 1930's and 1940's. System response was closely constrained by the assumed condition of pre-perturbation population equilibrium which implied relatively stringent relationships between variable model parameters. By including this information directly in the problem formulation, the predictive power of the N L P approach appears to be substantially increased over other approaches. ACKNOWLEDGMENTS This work was based on part of a dissertation submitted in partial satisfaction of the requirements for the Ph.D degree at the University of California, Los Angeles. Dr. S. Jacobsen, chairman of the dissertation committee, provided guidance throughout these studies. Dr. J. Kremer (Department of Biological Sciences, University of Southern California) reviewed early versions of this paper and made helpful comments. REFERENCES Agnew, T.T., 1979. Optimal exploitation of a fishery employing a non-linear harvesting function. Ecol. Modelling, 6: 47-57. Anderson, K.P. and Ursin, E., 1977. A multispecies extension to the Beverton and Holt theory of fishing, with accounts of phosphorus circulation and primary production. Medd. Dan. Fisk. Havunders., 7: 319-435. Atkinson, C.A., 1980. Analysis of perturbed dynamic systems under parameter uncertainty--a nonlinear programming approach with applications to marine ecosystems. Ph.D. Thesis, University of California, Los Angeles, CA, 188 pp. Atkinson, C.A., 1986. Discrete-time difference model for simulating interacting fish populations. Fish. Bull., 84: 535-548. Beverton, R.J.H. and Holt, S.H., 1957. On the dynamics of exploited fish populations. Minist. Agric. Fish. Food (G.B.). Ser. II, 18: 21-135. Clark, F.N. and Phillips, J.B., 1952. The northern anchovy in the California fishery. Calif. Fish Game, 38: 189-207. Food and Agricultural Organization of the United Nations, 1978. Some scientific problems of multispecies fisheries. FAO Fish. Tech. Pap. 181, 42 pp. Goodall, D.W., 1972. Building and testing ecosystem models. In: J.N.R. Jeffers (Editor), Mathematical Models in Ecology. Oxford Press/Blackwell, Oxford, pp. 173-183. Gulland, J.A., 1971. Ecological aspects of fishing research. Adv. Ecol. Res., 7: 115-176. Huppert, D.D., MacCall, A.D., Stauffer, G.F., Parker, K.R., McMillan, J.A. and Frey, H.W., 1980. California's northem anchovy fishery: biological and economic basis for fishery management. NOAA-TN-NMFS-SWFC-1, NMFS Southwest Fishery Center, La Jolla, CA.

27 Jaquette, D.L., 1972. Mathematical models for controlling growing biological populations: a survey. Oper. Res., 20: 1142-1151. Jeffers, J.R.N., 1972. Discussion of the application of mathematical models to a specific ecological problem. In: J.N.R. Jeffers (Editor), Mathematical Models in Ecology. Oxford Press/Blackwell, Oxford, pp. 345-354. Kremer, J., 1983. Ecological implications of parameter uncertainty in stochastic simulation. Ecol. Modelling, 18: 187-207. Laevastu, T. and Favorite, F., 1978. Fluctuations in Pacific herring stock in the Eastern Bering Sea as revealed by ecosystem model (DYNUMES III). Pap. 31, Symposium on the Biological Basis of Pelagic Stock Management, International Council for the Exploration of the Sea. Lasdon, L.S., Warren, A.D., Jain, A. and Ratner, M., 1978. Design and testing of a generalized reduced gradient code for nonlinear programming. ACM Trans. Math. Software, 4: 34-50. Lasker, R., 1978. The relation between oceanographic conditions and larval anchovy food in the California Current: identification of factors contributing to recruitment. Rapp. P.-V. Reun. Cons. Int. Explor. Mer, 173: 212-230. Ludwig, D., 1979. Optimal harvesting of a randomly fluctuating resource. I: Applications of perturbation methods. SIAM J. Appl. Math., 37: 166-183. Luenberger, D.G., 1973. Introduction to Linear and Nonlinear Programming. AddisonWesley, Reading, MA, pp. 240-278. MacCall, A.D., 1979. Population estimates for the waning years of the Pacific sardine fishery. Calif. Coop. Oceanic Fish. Invest. Rep., 20: 72-82. May, R.M., Beddington, J.R., Clark, C.W., Holt, S.J. and Laws, R.M., 1979. Management of multispecies fisheries. Science, 205: 267-277. Mendelssohn, R., 1976. Optimal problems associated with a Leslie matrix. Am. Nat., 110: 339-349. Methot, R.D., 1983. Seasonal variations in survival of larval northern anchovy, Engraulis mordax, estimated from the age distribution of juveniles. Fish. Bull., 81: 141-150. Murphy, G.I., 1966. Population biology of the Pacific sardine. Proc. Calif. Acad. Sci., 34: 1-84. Murphy, G.I., 1967. Vital statistics of the Pacific sardine and population consequences. Ecology, 48: 731-736. Nakamura, M., 1973. Some programming problems in population projection. Oper. Res., 21: 1048-1062. Oak Ridge Systems Ecology Group, 1975. Dynamic ecosystem models: progress and challenges. In: S.A. Levin (Editor), Ecosystem Analysis and Prediction. SIAM, Philadelphia, PA, pp. 280-296. Palm, W.J., 1975. Fishery regulation via optimal control theory. Fish. Bull., 73: 830-837. Parrish, J.D., 1975. Marine trophic interactions by dynamic simulation of fish species. Fish. Bull., 73: 695-716. Parrish, R.H. and MacCall, A.D., 1978. Climatic variations and exploitation in the Pacific mackerel fishery. Calif. Dep. Fish and Game, Fish Bull., 167: 11-49. Parrish, R.H., Nelson, C.S. and Bakun, A., 1981. Transport mechanisms and reproductive success of fishes in the California Current. Biol. Oceano., 1: 175-203. Pielou, E.C., 1977. Mathematical Ecology. Wiley, New York, pp. 8-110. Ricker, W.E., 1958. Handbook of computations for biological statistics of fish populations. Fish. Res. Board Can. Bull. 119, 300 pp.

28 Riffenburgh, R.H., 1969. A stochastic model of inter-population dynamics in marine ecology. J. Fish. Res. Board Can., 26: 2843-2880. Saila, S.B. and Parrish, J.D., 1972. Exploitation effects upon interspecific relationships in marine ecosystems. Fish. Bull., 70: 383-393. Sette, O.E., 1969. A perspective of a multi-species fishery Calif. Coop. Oceanic Fish. Invest. Rep., 13: 81-87. Soutar, A., 1966. The accumulation of fish debris in certain California coastal sediments. Calif. Coop. Oceanic Fish. Invest. Rep., 11: 136-139. Steele, J.H., 1979. Some problems in the management of marine resources. Appl. Biol., 4: 103-140. Vincent, T.L., Lee, C.S., Pulliam, H.R. and Everett, L.G., 1975. Applications of optimal control to the modeling and management of ecosystems. Simulation, 24: 64-72.