J, rheor. Biol. (1977) 69, 321-344

‘I’heory and Application of the Measurement of Structure and Determination of Function for Laboratory Ecosystems/ROBERT J. MULHOLLAND

AND

CAROLYNE

M. GOWDY

Center -for Systems Science, Oklahoma State University, Stillwater, Oklahoma 74074, U.S.A. (Received 24 May 1976, and in revisedform 6 May 1977) This paper concerns the development of an algorithm for the measurement of ecosystem structure designed for the purpose of computing ecosystem function. The use of the algorithm is based upon a new deterministic theory of sampling which gives bounds on the sample period for data collection. These bounds, which depend upon the arithmetic precision of the data, the dimension of the ecosystem and its maximum turnover, can be calculated before experiments are begun. The application of the sampling theory and identification algorithm developed is primarily intended for laboratory ecosystems. However, extension of these results to natural ecosystems is considered in detail. Several examples are provided for both laboratory and natural ecosystems illustrating the use of the results presented in this paper.

1. Introduction During

the past decade

the development

of state

models

for

ecological

systems has progressed to a high degree of sophistication. Recent surveys edited by Levin (1975) and Patten (1974, 1976) give some indication of the technical

level

achieved

in the application

of systems analysis

to ecology.

As a result, systems ecology has become a viable “bisociation” (Patten, 1971, p. xiv), firmly based upon the fundamental principles which underlie all natural science and the mathematical principles of system theory. For many years laboratory ecosystems, including microcosms and controlled semi-natural or artikial systems on a small scale, have served as useful models for the dynamical behavior of natural ecological systems. In light of the recent results obtained in systems ecology, it is perhaps useful to consider again the analysis of laboratory ecosystems, thus coupling 7 Supported in part by the National E:NG 7505341.

Science Foundation 321

through the research grant

322

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

research in biological modeling (microcosms) to that of mathematical modeling, the ultimate objective being to simulate and predict the behavior of ecosystems in response to natural and perturbed environmental conditions. In this paper the idea is presented that system analysis and mathematical modeling can be employed as a guide for the measurement of ecological structure in laboratory microcosms. This paper concerns the development of a sampling theory for identifying state models of ecological systems. The identification algorithm presented is not new, but rather it is a basic principle of system theory which has been discussed in the engineering literature by Lee (1964) and Sherif & Wu (1974). The associated sampling theory proposed is, however, unique in its approach to the data analysis problem. Unlike most engineering problems, the sampling frequency is not generally the only limiting factor in the analysis of ecosystems. In other words, it is possible to sample both too slowly and too rapidly. It is not sufficient to merely fix the sampling frequency according to the Nyquist rate (Carlson, 1975, p. 160), and then decide the arithmetic precision of the quantization based upon required accuracy or fidelity criteria. For ecosystems the Nyquist rate is rarely known, and furthermore the sampling frequency critically depends upon the quantization of the samples. It is shown later that upper and lower bounds on the sample period for the identification of time invariant models for ecosystems can be expressed in terms of the arithmetic precision of the samples, the dimension of the system, and some intuitive bounds on the systems dynamics. These bounds provide a means for determining the sampling period before data collection for system identification begins. While the sampling theory proposed is based upon the specific identification scheme developed by Lee (1964), the research results are believed to provide a strong tie between theory and experiment for ecological systems. Other related methods of systems identification for compartment models exist, including schemes presented by Halfon (1974, 1975), Leary & Skog (1972), Milanese & Molino (1973, and Cobelli & Romanin-Jacur (1976), but their application is limited because of the absence of a general sampling theory to guide data collection. There are specific theories of sampling based upon test input signals designed for system identification. For example, Ng & Goodwin (1976) have considered a sampling strategy based upon methods of experimental design as surveyed by Mehra (1974); however, these methods rely for the most part on the Nyquist theory with no consideration given to limiting effects inherent in the system. Therefore, it is hoped that the proposed general sampling theory will provide a basis for the implementation of specialized, and perhaps more effective, identification schemes (see Astrom

ECOLOGICAL

SAMPLING

THEORY

323

& Eykhoff, 1971; Nieman, Fisher & Seborg, 1971; and Bellman & Astrom, 1970). It is proposed that the substantial advances in the mathematical modeling of ecological systems be applied to the analysis of microcosm experiments. To this end the following discussion presents an outline of a theory for the measurement of structure and determination of function using tracer analysis for laboratory ecosystems. It is hoped this theory can be extended to future applications in the sampling, measurement and analysis of large scale natural ecosystems.

2. Modeling

Ecosystem Structure and Function

According to Odum (1962) the discipline of ecology is based upon the study of structure and function in ecosystems, i.e. nature. The composition of the biological community, the quantity and distribution of abiotic materials and the range of environmental conditions all constitute ecological structure. Ecological function within the ecosystem is defined in part by the throughput of biological energy, the rate of material cycling, and environmental regulation and reaction. These statements imply that structure is implicity a static concept, while function is dynamic being dependent upon time as a variable. However, both ecosystem structure and function are ultimately time dependent because they are forced by environmental inputs which are time varying. At first, ecological research was concerned with structure using what could be called the descriptive approach and leading from mere species lists to’ more mathematical realizations of species diversity involving, for example, Shannon’s formula. The functional approach came later, dealing with studies ranging from production to community metabolism to total energy budgets for ecosystems. The search for relationships between structure and function under natural and perturbed environmental conditions presently dominates much of ecosystem research (Smith, 1975). Odum (1971, p. 17) gives the best known model as: throughput

= turnover x content.

(1)

This equation defines ecological function to be linearly proportional to ecological structure, with the constant of proportionality equal to turnover. In steady-state systems, turnover equals the reciprocal of the time required for throughput to completely replace compartmental content. An example of the use of equation (1) in ecosystem analysis is provided by the compartment modeling principle. A compartmental model is defined by considering the component parts of an ecosystem as lumped entities, i.e.

324

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

the compartments, which are able to receive and distribute energy, biomass, or materials and which accumulate these quantities at a rate proportional to the net balance of the inflows and outflows. Each compartment in the ecosystem is a storage element for the quantities flowing into and out of it. When a conservation of matter or energy law is applied to the system, a rate equation for each compartment is obtained. If the compartmental outflows are constrained to be linearly proportional to the amount of stored matter or energy in the donor (or source) compartment, then the so-called linear donor controlled model results. This type of compartment model represents the ecological condition of ultimate resource limitation (Patten, 1975). The mathematical model is realized by the following expression : ii = ~

(i = 1, . . . , n),

afjXj-a*fXj+uf

j=l

(2)

where there are n compartments in the ecosystem, the prime notation indicates summation for all j # i, U,~ is the rate coefficient for transfer from compartment j to compartment i, Xi represents the content of compartment i and Zi its time rate of change, Ui expresses any possible exogenous inputs to compartment i, and aii

=

i;

aji+aOi9

(3)

where U,i indicates a transfer from compartment i out of the ecosystem to the environment. The compartment model described by equation (2) has the properties of aij 2 0 for i #j, uii > 0, and Ui 2 0, from which it can be proven that initially positive compartmental contents, x,(O) > 0 (i = I, . . ., n), implies positive contents for all future time, xi(t) > 0 (i = 1, . . ., n) for all t 2 0 (Mulholland & Keener, 1974). It can also be shown that Uii represents the turnover for compartment i, or aii = l/~i

(i = I ). * .) n),

(4) where Ti is the turnover time for compartment i. There are essentially two interpretations of turnover involving either steady-state or transient ecosystem operating conditions. In steady-state each compartment in the ecosystem is balanced with respect to inputs and outputs, and the content of each compartment is constant. Hence, z?* = 0 for all t, and from equation (3) II

affXf

=

c' I-1

UfjXj+Ui.

(5)

ECOLOGICAL

In equation (5), the right side partment i as the sum of the (summation). In steady-state, turn both equal throughput. becomes : aiixi=zi

SAMPLING

THEORY

325

of the equality expressed total influx to comexogenous input (Ui) and endogenous inputs compartmental input equals output which in Denoting throughput by zI, equation (5) (i=

l,....n).

(6)

which is clearly the same as equation (I) when a,, is interpreted as the compartmental turnover. If compartment i is assumed to be isolated from the rest of the ecosystem. then equation (2) becomes: ki = -Uiixi. Note that both exogenous and endogenous inputs to compartment zero. The solution of equation (7) is: xdt) = xi(O) e-‘jT1,

(7)

i are now (8)

where Xi(O) is the compartment content at t = 0, and ai, is given by equation (4). Equation (8) prescribes the compartment content at time t = Ti to be 36.8 % of the initial value. The turnover time Ti also prescribes the half-life T,, of the contents of the isolated compartment to be: To = Ti In 2.

(9)

The steady-state and transient conditions described often do not represent natural ecological phenomena involving the throughput and content of energy, biomass, or material. However, when radioactive isotope labels (tracers) are superimposed upon the material flows in ecosystems, these conditions are often observed when equation (2) is taken as the material balance equation for the label. The basic assumption is that equation (2) holds equally for the material as for its label so that the rate coefficients (aij) of material transfers are identical. This implies the material and its isotope label are ecologically equivalent. Mathematically, the use of isotope labels corresponds to a linearizing assumption mainly because the superimposed isotope levels are generally small in magnitude when compared to the material transfers within the ec:osystem. Thus, even though the ecosystem structure and function are known to obey a steady-state non-linear relationship, the introduction of the isotope represents a perturbation for which the linear approximation (2) holds true (Patten, 1975). Through the use of radioactive tracers, experiments can be designed to measure the rate coefficients of equation (2) given time series data for the

326

R.

.I.

MULHOLLAND

AND

C.

M.

GOWDY

compartment contents and inputs. Compartment content plus the associated rate coefficients determine material transfers as defined by the assumed model for ecosystem structure and function (1) or equivalently (2). These ecosystems experiments differ from the tracer methods of physiology and biochemistry in that it is not always possible to inject labeled material directly into all compartments within the ecosystem. Indeed, natural inputs to ecosystems are typically limited to one or two compartments. This fact requires the development of different methods of experiments and sampling for ecosystem analysis. A consideration of these methods constitutes the subject of the following discussion. 3. Identification

Problem

For the purposes of this discussion, equation (2) will be accepted as a local model for ecosystem structure and function. Using vector-matrix notation equation (2) can be written as: k = Ax+Bu,

(10)

where the state x is an n-dimensional vector of the compartmental contents, the coefficient matrix A = (Q,j) is an ordered array of the compartmental rate coefficients, the input u defines the ecosystem inputs as an m-dimensional vector, and the input connection matrix B = (b,j) is of n rows and m columns with entry bij equal to one if input Uj enters compartment i, if not, the entry is zero. In equation (IO), the inputs u and their entry into the ecosystem compartments, the B matrix, are known. The identification problem is defined for (10) as a method to determine the entries in the A matrix of rate coefficients given sufficient measurements of the state x. The tracer method of biochemistry and physiology provides a solution to the identification problem when inputs to all n compartments are available (Atkins, 1969), i.e. when B = I, the nth order identity matrix. For ecosystems m < n generally, therefore the tracer method cannot be fully employed. Consider the following experiment in which a constant (known) input is applied to the ecosystem modeled by equation (10). After sufficient time has elapsed, a steady-state is observed, as indicated by Z = 0. Denote this steady-state by the vector xS, which is the solution of: 0 = Ax,+Bu.

(11)

In the experiment under discussion, the input is assumed to be a constant rate of infusion of labeled material. It is also assumed the steady-state can be measured using standard radiological methods. On a new time scale,

ECOLOGICAL

SAMPLING

THEORY

377

the experiment is re-started at time t = 0 by removing all inputs to the ecosystem. Thus, before 1 = 0 the system is in steady-state x(0) = x,, while for t > 0 the state of the system is governed by the solution of: ;ir = Ax,

(12)

blecause u = 0 for I > 0. The labeled contents of the compartments of t!he ecosystem start at t = 0 with values determined by x, and decrease exponentially as t becomes large as all the labeled material flushes from the system along the pathways for material transfer. It is assumed that input u which generates the steady-state x, can bc chosen so that x, does not lie in a proper subspace of the n-dimensional state space of solutions of equation (10). This assumption is satisfied by systems (10) which are completely controllable, that is systems capable of being dlriven by inputs to any state in state space. Johnson (1976) shows that open systems are generally controllable if the matrix of the intercompartrnental rate coefficients (which exclude environmental transfers) is of rank fir-I. Despite the fact the state of the system, under the experimental conditions described, is governed by the differential equation (12) the problem is basically one of algebra. This is because the solution of equation (12) is prescribed by: x(t) = 4(t)x(O), where the matrix exponential

(13)

is denoted by the n x )z matrix : 4(t)

= exp (At).

(14)

and called the state transition matrix. Thus, any future state x(t) at t > 0 is defined by the initial state x(O) through the transformation 4(t) defined by equation (14). Consider equation (13) with t = T, for which, X(T) = 4(T)X(O).

(15)

Now let t = 2r, giving: x(2z) = @(2z)x(O) = aqz)x(O)

= cb(T)X(T).

(16)

Thus, for any positive integer i, it is clear that: x(h) = 4 x[(i-

l)r]

(i = 1, 2, . . .),

(17) of this difference

where Q, = @p(r) and (13) with t = k is the solution equation. A continuous record of the state x(r) is not required in order to determine @, which through equation (14) gives A. Since Q, is an n x n constant matrix,

328

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

contraint equations numbering n2 are required to prescribe the entries in this matrix. The problem reduces to finding A from a finite number of discrete and regularly spaced samples of the state, given at times t = 0, 7, 22, . . .) kz where T is the sample period. Substitution of these state samples into equation (17) gives: x(z) = @x(O) x(2r) = @x(z) (18) x(h)

= cDxc(k:

1)z-J.

I This linear algebraic system contains n2 unknowns (all 4v entries of 0) and k independent vector equations of 12components each, for a total of kn equations. For k < n no unique solution of system (18) exists, and the system is over determined for k > n. Thus, k = n or (n+ 1) discrete samples of the state are necessary and sufficient to determine Cp. And, from equation (14): A = (l/z) In a, (19) which identifies the rate coefficient matrix. It is possible to re-write system (18) in the following algebraic form : x, = @Xl,

more standard (20)

where X, is an IZx k matrix formed by columns from the first k state samples, and Xz is similarly sample period.

X, = {x(O), ~(4, . . ., x[(k- %I), (21) obtained from the next k state samples delayed by one X, = (x(x), x(2r), . . ., x(kz)}.

(22)

For k = IZ in equation (20) the solution for CDis obtained by inverting the matrix X,. This assumes n linearly independent state samples which is indeed the case when the ecosystem model defined by equation (10) is completely controllable. Thus, the identification problem is solved by: A = (l/z) In (X,X;‘) (23) under these stated conditions. It is mathematically necessary and sufficient to make (n-l- 1) discrete state measurements in order to solve for Cp. At least (n+ 1) measurements are required for a unique solution, however on occasion it may be desirable to obtain more than (n+ 1) state samples. For example, when experimental

ECOLOGICAL

SAMPLING

THEORY

329

error is suspected in the state data, it is necessary to increase the dimension of the data set beyond (n+ 1). Under these conditions k > it, and there are more equations than unknowns in (18). A Gaussian least squares estimate for the solution of (20) can then be obtained by multiplying (on the right) by the matrix transpose of X,, denoted by XT. This results in : x2x;

= @x1x;.

(24)

The model identified by the solution of equation (24) is then given by: A = (l/r) In [@,XT)(X,XT)-‘1.

(25)

This model produces a minimum in the error defined as the square of the difference between predicted state vaiues and measured state vaIues on the time interval t = 0 to t = kz (Lee, 1964). It should be noted that some care is required in the application of equation (119) to ecosystem analysis. The identification scheme discussed is unconstrained with respect to the entries in the A matrix, that is equation (3) does not constrain the columns of the identified matrix. Also, equation (1) is an approximate linear relationship, which is not generally viewed as being universally applicable to ecosystem analysis (Blesdoe, 1976). Thus, transfer matrices with confused ecological meaning are possible as outcomes of the computation defined by equation (19). However, the ecological validity of tlhe identified A matrix always can be judged by whether or not the constraint (3) holds. The problems associated with unconstrained identification are known (Halfon, 1974), and constrained identification methods which for the most part overcome these problems exist, e.g. Halfon (1975). Indeed, the single dose and constant infusion methods of biochemistry and physiology (Shipley & Clark, 1972) are both constrained by equation (3). The obvious dlisadvantage of the constrained identification methods is that they force all ecosystems to obey equation (1). A detailed comparison of the applicability of the constrained versus unconstrained identification methods to ecosystem analysis is not available. However, both methods are hampered by the lack of a suitable sampling theory which is the subject of the remaining portion of this paper. 4. Sampling Problem This section is concerned with experimental design and limitations of the theory. Problems of precision and accuracy of the state samples are of prime concern. Because the periodic state samples are obtained from radiological measurement instruments, these data are of finite arithmetic precision. In other words, each sample x(kz) is a number rounded to a fixed number of

330

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

decimal places. This round-off error affects the sample period r, or the rapidity with which the ecosystem is sampled, as required for identification. The sample period is also affected by the maximum and minimum compartmental turnovers within the ecosystem. These turnovers indicate the extent of the ecosystem dynamic behavior, hence their relationship to the sample period is clear. In that which follows, the relationships between the dynamic range of compartmental turnovers, the finite arithmetic precision of state samples, and the sample period are explored. These relationships are to be used as a guide for choosing the proper sample period for the solution of the identification problem. The state transition matrix of (14) is given by: a)(t) = LY(t)L-1,

(26)

Y(t) = exp (At)

(27)

where :

and L is an invertible n x n matrix of the eigenvectors of A, and A is a diagonal n x n matrix of the eigenvalues, so that: L-‘AL

= A.

08)

The n eigenvalues of A are defined as the roots of the polynominal: IIA--A( = 0.

(29)

These eigenvalues, which describe the ecosystem dynamics, are related to the compartmental turnovers. This relationship, in part, follows from a theorem by the Russian mathematician Gerschgorin: the eigenvalues &, . . ., ;1, of the matrix A = (Q) lie on one of the circular regions of the complex A-plane described by :

The right-side of (30) is the sum of the absolute values of all the entries in the jth column of A except the prinicpal diagonal entry - Ujj. Let (31) Then, (30) is the closed circular region in the I-plane of radius rj and center -ajj. It should be noted that the theorem is also true for row sums (Wilf, 1962, p. 39). In the following development it is assumed the eigenvalues of the A matrix are real. This assumption may not be necessary but it certainly

ECOLOGICAL

SAMPLING

THEORY

3311

sirnplifies the analysis of the sampling problem. For example, the controllability of the continuous system (10) with real eigenvalues implies the controllability of the discrete system (18) regardless of the sample period value r (Kalman, 1963). The conditions for the eigenvalues of A to be real have been presented by Hearon (1963) and discussed in detail by Thron (1972). Shipley & Clark (1972) and Atkins (1969) along with many other authors of books on the tracer method accept these conditions for real eigenvalues as not being unusual. The A-matrices of ecosystems described by system (2) are column diagonshy dominant, that is ajj > 0 and lajjl 2 if

Jaijl = rj

(j = 1, . . . , n).

i=l

Let dj = Iujj]-rj

(j = 1, . . ., n).

(33)

Then, 6j 2 0 defines the distance between the jth Gerschgorin circle and the Im {1”} axis of the I-plane as shown in Fig. I. As shown, these circles

FIG.

1. The jth Gerschgorin circle plotted in the complex I-plane.

lie in the left-half of the A-plane. This follows from the fact that -ujj (the center of the circle) is a negative real number which is in magnitude greater than or equal to the radius rj. Ifaj =Oforallj= 1, . . .. n, then the ecosystem is closed with respect toI the material flows defined. Under this condition, all Gerschgorin circles intersect d = 0, the A-matrix is singular, and hence the set of eigenvalues includes zero. Very few model ecosystems, even laboratory microcosms,

332

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

are truly closed. Therefore, at least one of the terms 6, will be assumed to be nonzero, which implies the solution of equation (12) is a bounded function of time and x(t) + 0, at t + co, because all the eigenvalues of A are negative real numbers (Hearon, 1963). That is, equation (12) is asymptotically stable. This enables the set of eigenvalues of A to be ordered and bounded as follows : -c! < ?q < 3.2 < . . . < 1. < -6.

(34)

The case of repeated eigenvalues is ignored. The eigenvalue bounds are: (35)

and 6 = min

l 0 implies xi(t) > 0 for all t 2 0. Assuming Cij 2 0 for all j = I, . . . n yields: exp (-C(t) I Xi(t) < exp (-st), (38) where each exponent has been replaced by equations (35) and (36) and the constraint

has been used. The collection of time series data for each compartment in the ecosystem takes place in two distinct steps, in which first time and then amplitude are

ECOLOGICAL

SAMPLING

THEORY

333

quantized. The samples of the state are obtained only at discrete times. These samples are periodically spaced apart by r time units giving a data base defined by {x(O), x(r), . . ., x[(k+ l)r]}. It has been shown the minimum data base for solution of the identification problem is (n+ 1) state sa.mples or k = n. Next, the finite precision with which the measurements of state are made must be accounted for in terms of a mathematical relationship. The assumption is that x,(t) and its samples are fixed point numbers which are less than or equaI to one. These samples, arising from radiological measurement equipment, are data with no more than a fixed number of digits to the right of the decimal point. Let IOVd, where d is a positive integer,

Fro. 2. Functional relationship for a quantizer which truncates and where x represents the input variable and y the output.

be the smallest number in the data set. Then, the quantizing operation on the samples is described by the function y = f(x) as plotted in Fig. 2, or by the relations : y=o for 0 I x < 10Sd, y = 1o-d for 10Wd I x < 2x lOed, y = 2x10ed for 2~10~~ Ix < 3x10Vd, etc. These quantized samples assume a pure truncation scheme for the finite precision arithmetic. Depending upon the measurement equipment used, other arithmetic schemes may be necessary. For example, it is common to round 0.478554 to either O-4786 or 0.4785 as a four place number depending upon whether or not the last digit is even (or odd). The round-off error fol such a scheme is IEI I 0-5x lOed while the proposed truncation scheme

334

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

gives an error of 0 5 (El < 10Td. Since the samples are positive numbers, the pure truncation scheme was selected for reasons of symmetry to represent the quantizer. The relationship between the sample period, the quantization of the samples, and eigenvalue bounds is of interest. As has been pointed out, it is possible for the sample period to be too small. This occurs when x,(t) and x,(t+r) have the same arithmetic precision, that is when .u,(t+r) is truncated to d digits it equals the truncated value of xi(t). Since (38) gives an upper bound for xi(t), a measure of the minimum sampling time, denoted by r-, is given by: 1 - 10md = exp (-6~~) 2 s,(r-), (40) or solving for the minimum

sampling period yields:

z- = (- l/6) In (1 - 10-3.

(41) If z > r- the state samples are independent, but an upper bound on z exists. Suppose z is chosen large enough so that x&), the last sample in the data set, is truncated to zero by the quantizer, then the data set is incomplete and the identification problem cannot be solved. Thus, the maximum sample period, denoted by rf , must be such that the (n+ 1)th sample is at least equal to the smallest number in the finite precision arithmetic, that is .~i(nZ’) 2 10Td. But, equation (38) gives a lower bound for xi(t), thus the maximum sample period is easily computed from: 10bd = exp (+ctnr+)

I x,(nr+),

(42)

which gives : 7+ = (-l/cm) In (10-1. (43) As an example of the application of these equations, consider a system with eigenvalue bounds of CI = 190 and 6 = 2. Suppose samples of the state of this system can be made with four place precision, so d = 4. Then, from equations (41) and (43) z- = 5 x 10m5 and z+ = 0.012 are obtained. A reasonable choice for z is the geometric mean of rf and z-, or z = (z+z-)+,

(4) choice is

which for this example gives z = OGOO78. A more practical obtained by rounding z to O*OOl. Equation (44) gives a relationship for choosing the sampling period for the identification problem which is valid as long as rt- < t+, or

6/m > In (1- lo-31 In (10-3. (45) The eigenvalue bounds (S and or)and the number of ecosystem compartments (n) are fixed constants, hence a computation yielding r- > z + implies

ECOLOGICAL

SAMPLING

THEORY

335

the precision of the measurements is too low to solve the identification problem. This computation can be made prior to the collection of data, thus validating the sampling program. For the example under consideration, elquation (45) holds for all integers d 2 3. Therefore, radiological measurement equipment with at least three place accuracy for fixed point numbers is required for sampling the state of this example system. In order to compute the bounds z+ and T- for the sample period r it is necessary to know the eigenvalue bounds c1and 6. However, computations of these bounds using equations (35) and (36) requires knowledge of the entries in A, which are unknown. Based upon intuition, experience, or isolated laboratory experiments ecologists tend to know the turnovers for the species present in an ecosystem. Such knowledge is at least accurate to an order of magnitude. Thus, a relationship between the turnovers and the eigenvalues would enable the computation of the eigenvalue bounds, and in turn determine rf and z-. Consider again the compartment model defined by equation (2) where the tlurnover time Ti is defined by equation (4). From the ecologically unrealistic assumption that no coupling exists between the n compartments of the ecosystem, the eigenvalues can be shown to equal the compartmental turnovers. Consider equation (8) as the solution of equation (2) under these conditions. The eigenvalues are I+ = T,yl for i = 1. . ., n. For this case the Gerschgorin circles are degenerate, with all radii equal zero, consisting of the n points given by TL1 for i = 1, . . ., n. From (31) the Gerschgorin circle radii are clearly determined by the pattern and amount of coupling between ecosystem compartments. The Gerschgorin circle centers are equal to the diagonal elements of A, which regardless of intercompartmental coupling equal the turnovers as defined by (4). Since all the Gerschgorin circles are known to lie in the left-half of the J-plane, the center magnitude must be greater than or equal the center radius : Ti’ = lUjjl 2 rj, (46) for each circle j = I, . , ., n. Let T, 1 be the maximum turnover, then it is clear that all eigenvalues of the system lie in the circle of center and radius lr;l. Hence, 2T;’ 2 c( and -2T; ’ I 2, < AI < . . < il,t < 0. Using 3’0 ’ as an estimate for CIin equation (43) yields : z+ = (- 1/2nT;‘)

In (lOmd).

(47 )

Given the number of compartments in the ecosystem, the maximum turnover, and the precision with which the measurements of the system state are to be made, equation (47) prescribes the maximum sample period. 22 ‘T.B.

336

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

The minimum sample period is not easily computed without the eigenvalue bound of equation (36). Fortunately, it is the maximum sample period bound which is the more important extremal. Assuming sufficient arithmetic precision of the state samples to prescribe r- < z+, it is merely necessary to choose r I z+ in order to solve the identification problem. Exact knowledge of z- is of little use when d is large enough. Earlier the solution of the identification problem was explained to begin with an experiment which initially places the ecosystem in steady-state. The observation of this steady-state is obtained by noting pi = 0 for all i= I,..., n, or by condition (11). If both observations are made, sufficient data then exists to determine the steady-state condition and 6j for each compartment. In other words these data determine the compartmental lossesfO, and steady levels jzj, SO that 6, = fo//Zj forj = 1, . . ., n. The value 6 is then determined by equation (35) which in turn prescribes z- through equation (41). The measurement of Xj in steady-state is required to determine xj(0), the first sample for compartment j, thus the additional measurement offO, enables the computation of r-, which may be worth the additional effort in certain cases. 5. Applications As an elementary illustration of the application of the identification algorithm and sampling theory consider the recent study of phosphorus kinetics in freshwater microcosms provided by Sebetich (1975). This research concerns the transfer of radiophosphorus among lake water, diatom, snail and sand compartments. Using a standard technique pioneered by Patten & Witkamp (1967), the analog computer is employed to identify the transfer rates and to simulate the dynamic behavior of microcosms of varying complexity. The single and combined effects of diatoms, snails and sand on the movement of radiophosphorus through the system are noted. Sebetich presents time series data of radiophosphorus concentrations within each compartment of the microcosm system. These data are compared with analog computer simulations of a mathematical model of fixed structure and free parameters corresponding to the transfer rates of radiophosphorus. By manipulating the free model parameters a fit between the computer simulations and the time series data is obtained by trail and error. It is this approximate identification scheme which the results of this paper propose to replace. Table 1 gives the maximum turnover, the corresponding turnover time, and the maximum sample period, assuming precision d = 3, for the seven microcosm systems considered by Sebetich. Because the microcosms are

ECOLOGICAL

SAMPLING

THEORY

337

TABLE 1 Computed sampling times for the compartmental radiophosphorus microcosms qf varying complexity (Sebetich, 1975) Compartments Water-diatom Water-snail Water-sand Water-diatom-snail Water-diatom-sand Water-snail-sand Water-diatom-snail-sand

To-1

(hr-1)

04f589 O-0043 0.0633

a.0379 0~1ooo 0*0611 0.1035

To CM 14.5 233 15-8 264

IO.0 16.4

9.66

in

7+ 0-4 25-06 401-61 27.28 30.38 11.51 18.84 8.34

presumed to be closed systems no minimum sample period is computed. CIonsider the system in which diatoms (Nitzschia palea) take up phosphorus from lake water. Three samples spaced at intervals of 25 hours (maximum) are sufficient to prescribe the compartment model for the system. Data presented by Sebetich (1975, Fig. 7) is variously spaced from one hour to twenty-five hour intervals. On a time scale of 50 hours with a sample period (r) of ten hours, samples were extracted from these data and several associated compartment models were constructed using (23). The best of these three-point models indicated a mean square error of O-0574 when used to reproduce the entire 50 hour data set. The Gaussian least squares algorithm of (25) was used next to construct a compartmental model from all six data points on the 50 hour time scale. The mean square error incurred by this rnodel in reproducing the data base was ONl43, a significant, but expected, improvement in accuracy. The identification algorithm of (25) produced the following six-point compartment model : 1, = - 0*063x, +0*0054x, i2 = 0.052x,-0*016x, i.

(48)

Sebetich (1975) does not present details of the transfers produced by the analog computer identification scheme, but he does give values for the turnovers of O-069 and O*Olg for the diatoms and water respectively, which should be compared with O-063 and O-016 as obtained from (48). It is interesting to note that the model of (48) is not closed with respect to the phosphorus transfers. This is not unexpected mainly because a more accurate model for the uptake of phosphorus by algae includes dissolved inorganic,

338

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

organic and dissolved organic phosphorus compartments (Watt & Hayes, 1963). Table 2 provides a survey of computed sampling times for various aquatic and terrestrial microcosms based upon an assumed three place decimal precision for measured data. It should be noted that with the exception of Azotobacter the maximum sample period z+ is less than the reciprocal of the maximum turnover T,‘. Thus, without the detailed information necessary to employ equation (47) for the computation of r+, a rough estimate for the sample period is provided by TO. TARLE

2

Computed sampling times for elemental dynamics within selected laboratory ecosystems Ecosystem (type) Artificial pond Aquarium Aquarium Sea water Lake water Azotobacter cultures Terrestrial microcosm Liriodendron forest Terrestrial microcosm

Element transferred

(2) OG91 0.200 0.104 0.17 0.24

IO.99 540 9.66 5.87 4.17

0.0128

77.9

89.9

Visser et (I/. (I 973)

0.0313

31.9

22.1

Patten & Witkamp (1967)

ON)185

540

O+IO967

103

4.74 3.45 8.34 6.76 3.60

311 71.5

Whittaker (1961) Whittaker (1961) Sebetich (1975) Watt & Hayes (1963) Lean (1973)

Olson (1965) Witkamp & Frank (1970)

The sampling theory developed for implementing the identification algorithm presented in this paper is based upon the analysis of laboratory ecosystems. The basic idea being that such systems are isolated from the cyclic environmental inputs present in nature. As a first step in ecosystem modeling the cyclic behavior of environmental parameters is often ignored in favor of analysis with average inputs. This results in open systems with time invariant relationships between structure and function which have as linear approximations models of the form considered in equation (10). For such systems the sampling theory developed in this paper has application to the design of field experiments. Table 3 provides a comparison of sample periods for terrestrial (old fields and forests) and aquatic (lakes and springs)

Bog Lake Mendota Pond Lake Springs Spring Spring

Post oak-blackjack Mixed deciduous Oak-pine Rain Forest Liriodendron forest

Forest

Cedar Lake Lago Frains Silver Cone Root

Aquatic

Grazing food chain a b Arthropod food chain a b Vegetation a b

Old field

Ecosystem (type) 7+

0.957E-1 0.207EO

0.902E I 0.416El

0.302E2 0.811E2 0.437EO 0.380E3 064682

0.143E-1 0.266E-2 0.263El O.l30E-2 0.2OOE-2

0.834EO 0.517EO O.l04E-2 04@E-3 0.329E-2 0.620E-3 0.25OE-2

0.523E-2 0.742E-2

(Yr)

0.220E3 0.15583

0.138El 0.167El 0.370E3 0.130E4 0.210E3 0.351E3 0.138E3

3

O.lOOE-1 0.2OOE-2 O.l27E-1 0.210E-1 -

0,269E-2 0.1428-2 O.lOOE-2 0.48E-3 0.550E-3 0.58OE-3

0.391 E-2 0.391E-2

0.217E-3 0.256E-4

0.848E-3 0.556EA

O.l20E-1 0.230E-2 0.183EO 0.2OOE-2

0.474E-3 O,377E-1 O.l02E-2 044OE-3 O.l25E-2 0.58OE-3 0.120E-2

O.l93E-1 0.285E-1

0.107E-2 0.436E-3

0.273E-2 0.422E-3

times .for ecoq*sten? energy or biomass

0.879E-2 0.321 E-2

sampling

0.131E3 0.215E3

Computed

TABLE

et al. (1969)

(1960)

(1967) er al. (1970)

Reference

Johnson & Risser (1974) Satchel1 (1971) Woodwell (1970) Odum (1970) Shugart et al. (1976)

Lindeman (1942) Lindeman (1942) Emanuel & Mulholland Saunders (1971) Odum (1957) Tilly (1968) Teal (1957)

Galley Kelly

Menhinick van Hook

Golley (1960) Pearson (1964)

&namics

(1975)

340

R.

I.

MULHOLLAND

AND

C.

M.

GOWDY

natural systems. From the maximum turnovers, number of compartments and an arithmetic precision of three decimal places, the maximum sample period for each system is presented. Also, except for the Frains lake and Liriodendron forest studies, which represent closed material cycles, the minimum sample periods are computed. The sample period r of Table 3 is taken as the geometric mean of z+ and z- according to equation (44). The rain forest data provided by H. T. Odum (1970) presents a problem with d = 3, in that z+ < z-. This problem arises because the model includes a microbe compartment with an extremely high turnover typical of tropical rain forests. This is a severe limitation on the experimental procedure requiring four place data precision (d = 4) for which a sample of period r = 0.002 yr results. Table 3 is based in part upon data provided by O’Neill (1971). The application of the sampling theory developed in this paper to natural systems as depicted in Table 3 may shed new light on the study of ecosystem development. The value for the sample period provides a single number as a measure for the dynamics of the structure and function relationship at a particular stage of ecological succession. Some knowledge of the successional stage for each system in Table 3 is required for such a comparison. This represents a study unto itself beyond the intent of this paper. However, it should be noted that the sample period depends upon the maximum turnover which is believed to be correlated with the stages of ecosystem development (Odum, 1971, p. 252). This correlation is evident in the experimental data produced by Whittaker (1961) on the transfer of radiophosphorus within aquaria. Table 4 summarizes Whittaker’s data to show that in the developmental stages ecosystems require more rapid sampling techniques, while in the mature stage the system structure and function are such that a longer sample period will suffice. Again, a more detailed analysis of these relationships is indicated by these preliminary results.

TABLE

4

Computed sampling times for the compartmental radiophosphorus transfer in aquaria (Whittaker, 1961) Time period following inoculation

TO_ 1 (hr _ ‘)

Initial (O-7 hr) Intermediate (7-63 hr) Final (546 days)

r+ (hr)

TO 0x3

__.

.~ 0.313 0*0855 OGO91

3.19 11.7 109.9

._-~-.-_ 2.21 8.08 75.9

ECOLOGICAL

SAMPLING

THEORY

341

6. Summary and Conclusions

This paper has dealt with the development of an algorithm for the measurement of ecosystem structure in order to determine ecosystem function. The implementation of the algorithm is based upon a deterministic theory of sampling which gives bounds on the sample period for data collection. These bounds, which depend upon the arithmetic data precision, the dimension of the system and its maximum turnover, are calculated before the experiment or field work is initiated. The application of the sampling theory and identification algorithm is primarily intended for laboratory ecosystems, however when cyclic environmental effects are considered only as average inputs these results apply equally to the analysis of natural systems. Several examples of both laboratory and natural ecosystems illustrating the use of the sampling theory and identification algorithm have been presented. However, the real usefulness of the theory remains to be discovered in the design of future ecological experiments. Because of the interactive nature of ecosystems, the turnover of species within a system differs from that observed for isolated species. Observations of ecological systems depend upon these interactions, however the system dynamics can be estimated from the isolated species behavior. Referring to equation (47), which estimates the sample period, the term T; i defines the turnover of the fastest species in the system, n is a measure of the ecosystem size and hence of all its possible interactions, and d gives the precision with which measurements can be made. Since the sample period is inversely proportional to n, the interactive nature of the ecosystem is clear. The problem of extremely small sample periods required for very large ecosystems can be to a certain extent compensated for by increasing the precision of the measuring instruments. But, fundamental limits on precision exist, necessitating simultaneous measurement and analysis within subsystems of lower dimension. It should be noted that Smith (1975) using similar concepts including the characterization of fast and slow turnovers, argues for theoretical relationships between evolution and ecosystem science. Thus, the sampling theory presented in this paper provides a potential tie between theory and experiment at the ecosystem level. At the present time a large number of models exist for widely differing ecological systems (see Table 3). It is often true that these models are developed using incomplete data sets, therefore the problem of model validation arises (Caswell, 1976). The sampling theory presented in this paper provides a meaningful way to compare model output with measured time series data. The model dynamics are captured completely when the sample period is chosen according to equation (44). Of course, as indicated in equation (47).

342

R.

J.

MULHOLLAND

AND

C.

M.

COWDY

trade-offs between sample period length and data arithmetic precision are possible. The results presented in this paper are founded upon a linear model for relating ecosystem structure and function. Therefore, these results are ideally suited to the tracer analysis of ecosystems wherein the tracer element represents a small amplitude perturbation superimposed upon the material flows within the system. And, it follows that the application of the identification algorithm must result in a linear model for ecosystem function as computed from measurements of structure. For various trace loadings of the system, the tracer analysis described produces linear structure and function relationships

as local approximations

to the global nonlinear

relation-

ship which governs ecosystem behavior. If many such local determinations of ecosystem function can be obtained from structure measurements, then it is certainly possible to construct a table of these linear structure and function values which approximate the global (non-linear) relationship. By using the technique of descriptors of ecosystem attributes proposed by Smith (1970), it may be possible to construct a global non-linear model from this table relating ecosystem structure and function. Thus, even though the proposed identification scheme is firmly based upon the concepts of linear analysis, the technique developed, at least in theory, can be applied to the construction of a non-linear model which accurately describes ecosystem dynamics. The generality of this approach to systems modelling is intriguing, but its true usefulness with respect to applications remains an open question. The authors wish to thank the following people for their many helpful suggestions: Dr W. R. Emanuel, Environmental Sciences, Oak Ridge National Laboratory; Dr K. W. Thornton, Waterways Experiment Station, U.S. Army Corps of Engineers; W. M. Sharp, Biochemistry Department, Oklahoma State University; and Dr J. W. Falco, Environmental Research Laboratory, U.S. Environmental Protection Agency.

REFERENCES ASTROM, K. J. & EYKHOFF, P. (1971). A~fomatira 7, 123. ATKINS, G. L. (1969). Multicompartment Models for Biological BELLMAN, R. & ASTROM, K. J. (1970). Math. Biasci. 7, 329. BLEDSOE, L. J. (1976). In Systems Analysis and Simulation in

Systems. London:Methuen.

Ecology (B. C. Patten,ed.), Vol. IV, pp. 283-298.NewYork: AcademicPress. CARLSON, A. B. (1975).Communication Systems, 2ndedn.New York: McGraw-Hill. CASWELL, H. (1976). In Systems Analysis and Simulation in Ecology (B. C. Patten,ed.), Vol. IV, pp. 313-325.New York: AcademicPress. COBELLI, C. & ROMANIN-JACUR, EMANUEL, W. R. & MULHOLLAND,

98.

G. (1976). Math. Biosci. 30, 139. R. J. (1975). I.E.E.E. Trans. Autom. Control AC-20,

ECOLOGICAL

SAMPLING

‘THEORY

343

GOLLEY, F. B. (1960). Ecol. Monogr. 30, 187. HALFON, E. (1974). In Proc. First International Congress of Ecology (A. J. Cave, Conf. Sec.), pp. 262-265. Wageningen, The Netherlands: Centre for Agricultural Pub. and Dot. HALM)N, E. (1975). Simulation 38, 149. HEARON, J. Z. (1963). Ann. N. Y. Acad. Sci. 108, 36. JOHNSON, JOHNSON, KALMAN, KELLY,

F. L. & RISSER,

P. G. (1974).

Ecology

55, 1246.

L. E. (1976). Math. Eiosci. 30, 181. R. E. (1963). SIAM J. Control 1, 153. J. M., OPSTRUP, P. A., ORSON, J. S., AUERFJACH,

Models of Seasonal Primary Production in Eastern Ecosystems. Oak Ridge, Tenn.: ORNL Tech. Rpt. LEAN, D. R. S. (1973). Science, N. Y. 179, 678. LEARY, LEE, R.

R. A. &

53,

C. K.

Identificarion,

SKOG, K. (1972). Ecology (1964). Optimal Estimation,

S. I. & VAN

Tennessee

DYNE,

Festuca

4310.

969.

Cambridge,

and Control.

chusetts: M.I.T. Press. LEVIN, S. A. (ed.) (1975). Ecosystem Analysis and Prediction. Publications. LINDEMAN, R. L. (1942). Ecology 23, 399. MEHRA, R. K. (1974). I.E.E.E. Trans. Autonr. Control AC-19, 753.

NG,

T.

NIEMAN, ODUM, ODUM. ODUM, ODUM,

Philadelphia:

E. F. (1967). Ecof. Monogr. 37, 255. MOLINO, G. (1975). Math. Biosci. 26, 175. R. J. & KENNER, M. S. (1974). J. theor. Biol. 43, 105. S. 8c GOODWIN, G. C. (1976). Int. J. Control 23,459. R. E., FISHER, D. G. & &BORG, D. E. (1971). Int. J. Control

MENHINICK, MILANESE, MULHOLLAND,

M. (1969).

G.

and Andropogon

MassaSIAM

M. &

E. P. Jap. J. Ecol. 12, 108. E. P. (1971). Fundamentals of Eeolog,v,

3rd

edn.

Philadelphia:

13, 209.

Saunders.

Ecol. Mongr. 27, 55. (1970). In a Tropical Rain

H. T. (1957).

H. T. Forest (H. T. Odum, ed.), pp. 1191-1281. Washington, D.C. : U.S. Atomic Energy Commission. (h.SON, J. S. (1965). Hlth Phys. 11, 1385. O’NEILL, R. V. (1971). Examples of Ecologica/ Transfer Matrices. Oak Ridge, Term.: ORNL Tech. Rpt. IBP-71-3. PATTEN. B. C. (ed.) (1971). Systems Analysis and Simulation in Ecology. Vol. I. New York: Academic Press. PATTEN, B. C. (ed.) (1974). Systems Analysis and Simulation in Ecology, Vol. III. New York: Academic Press. PATTEN, B. C. (1975). Am. Nat. 109, 529. PATTEN, B. C. (ed.) (1976). Systems Analysis und Simulation in Ecology, Vol. IV. New York : Academic Press. PATTEN, B. C. & WITKAMP, M. (1967). Eco1og.v 48, 8 13. PEARSON, 0. P. (1964). 1. Mammal. 45, 177. SATCHELL, J. E, (1971). In Proc. Symp. Productivit), o/‘the Forest Ecosystems qf the World. Paris: UNESCO. G. W.

SAUNDERS,

(1971).

In

The Structure

and Function

of Fresh-Water

Microbiul

Com-

(J. Cairns, Jr., ed.), pp. 31-45. Charlottesville: University Virginia Press. SEBETJCH, M. J. (1975). Ecology 56, 1262. SHERIF, W. & Wu, M. Y. (1974). Int. J. Control 19, 185. SHIPLEY, R. A. & CLARK, R. E. (1972). Tracer Method for in vivo Kinetics. New York: Academic Press. munities

SHUGART,

H.

H.,

REICHLE,

D.

E., EDWARDS,

N.

T. & KERCHER,

J. R. (1976).

Ecology

57,

99. SMITH,

F.

E. (1970).

In

Analysis

of Temperate

pp. 7-18. New York: Springer-Verlag. SMITH, F. E. (1975). Bull. ecol. Sot. Am. 56, 2. TEAL, J. M. (1957). Ecol. Monogr. 27, 283.

Forest

Ecosystems

(D. E. Reichle, ed.),

344

R.

J.

MULHOLLAND

AND

C.

M.

GOWDY

THRON, C. D. (1972). Bull. Math. Biophys. 34, 277. TILLY, L. J. (1968). Ecol. Monogr. 38, 169. VAN HOOK, R. I., REICHLE, D. E. & AUERBACH, S. I. (1970). Energy and Nutrient Dynamics of Predator and Prey Arthropod Populations in a Grassland Ecosystem. Oak Ridge, Tenn. : ORNL Tech. Rpt. 4509. VISSER, S. A., WITKAMP, M. & DAHLMAN, R. C. (1973). PI. Soil 38, I. WATT, W. D. & HAYES, F. R. (1963). Limnol. Oceanogr. 8,276. WHITTAKER, R. H. (1961). Ecol. Mongr. 31, 157. WILF, H. S. (1962). Mathematics for the Physical Sciences. New York: Wiley. WITKAMP, M. & FRANK, M. L. (1970). Ecology 51,465. WOODWELL, G. M. (1970). Scient. Am. 223, 64.