It is possible to rewrite system (3.10) in the following more standard
algebraic form

X = &X, (3.11)
where X% is an n x k matrix formed by columns from the first k state samples,
X, = {x(0), x(t), ... . xf(k-1)t]}, (3.12)

and X is similarly obtained from the next k state samples delayed by one
sample period,

Xy = (x(t), x(2t), ... , x(ke)}. (3.13)

For k = n in (3.10) the solution for ® is obtained by inverting the matrix X,.
This assumes n linearly independent state samples which is indeed the case
when the ecosystem model defined by (3.1) applies. Thus, the identification
problem is solved by

A = (1/t) 1n(%X,7!) (3.14)

under these stated conditions.

It is mathematically necessary and sufficient to make (n+1) discrete
state measurements in order to solve for 6. 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 significant experi-
mental error is suspected in the state data, it may be of use to increase the
dimension of the data set beyond (n+1). Under these conditions k > n, and
there are more equations than unknowns in (3.10). A Gaussian least squares
estimate for the solution of A can then be obtained (Lee, 1964).

SIMULATED EXPERIMENT

Early in the project, it was surmized that problems might arise in get-
ting the data from the EPA planned microcosm experiments which were designed
to test the Mulholland-Gowdy sampling theory. So, a simulation of experimen-
tal data was devised using an ecosystem model. It is difficult to theoretic-
ally analyze the propagation of data noise through a nonlinear algorithm.
Since the identification algorithm is nonlinear, a simulation was used to
examine the effects of data noise errors (on the algorithm). In the simula-
tions, data noise was propagated through the algorithm to obtain the statis-
tics of individual parameters within the identified model.

The identification algorithm, defined by (3.11), requires discrete samples
of the system state: x = x(kt). Actually, the ensemble of these samples is
itself a sample function for a stochastic process. Because of noise inherent
in the microcosm processes and errors introduced by the measurement techniques
and instruments, the analysis of the basic experiment (defined in Table 1)
cannot be wholly based upon deterministic systems theory, but it must include
a statistical analysis. The statistical interpretation of the basic microcosm
experiment relies upon the fact that the initial and final steady-states are
random variables with measurable statistics. These steady-states are connec-
ted by a sample function of a stochastic process which is a statistical repre-

10


Theorem 5.1

A system given by (5.1) in which measurements are obtained from % accumu-
lation pools is controllable if the test matrix Pee (5.11) of the system is
such that 1

(i) The upper partition Py is a reduced controllability matrix PL

with 2k < n-1, where n is the rank of A and k is the reduced
order of the controllability matrix, and

a k .

(ii) CA*P, is of rank 2.

Proof: vk
Rewrite the indirect controllability matrix as

k n-1

Py Bs ABs ...5 AkB;  AKTBS 0.5 at lp 6.16)
Pp. = = 5.16
nz) Os KBs 2.05 KARTS KAKBs 20.5 kaP@B
Now if the matrix ia is reduced controllable, then Pe is of rank n where
= [B: AB: ...: ak
Le [B; AB; ...; A’B]. (5.17)
Since Py is of rank n, then Pa is of rank n and
E
Py = [I; 0] (5.18)

where I is the n-th order identity matrix. Next multiply the matrix Po by
cak, yielding
pp = Coakes cak*?, ... ca*kpy. (5.19)

If 2k < n-1, then 2k is at most n-2 so that if Pr is of rank 2, then

w= teakp. caktla. - cqne2
Py = [CASBs CAN"'Bs ...3 CA

8] (5.20)

is of rank 2. (Recall that the rank of a product is less than or equal to the
rank of either matrix. Therefore, the matrix C is at most of rank 2.) Notice
that Po can be rewritten as

Po = [0; CB; ...3 cak-1p, wea ed (5.21)

and if Pe is of rank n, then

p, £ fos 1 (5.22)

where I' is the %-th order identity matrix. Therefore

um
[Rn
lo

= = I" (5.23)
—1c Pp I

%
lo

31


Definition 5.2

A linear system is said to be completely output controllable if it is
possible to drive any initial output at ty to any final output in a finite
time interval.

Definition 5.3

__ A linear system is said to be completely observable if for every initial
time t, and some finite time Ty t, <t<T, every initial state at t, can be

determined from the knowledge of the output on t,<t<t.

Rigorous treatments on the properties of controllability and observabil-
ity can be found in system theory literature (Kalman, 1963; Luenberger, 1963;
Lee, 1964; Ogata, 1967; and Chen, 1970). In the following, these properties
will be considered as they apply to the previously discussed methods of
identification for compartmental models.

DIRECT METHOD

Assume that a single dose of tracer may be directly inserted via special
routes into system compartments and that tracer behavior can be experimentally
monitored in specified compartments when the system is switched to a tracer-
free environment. The controllability and observability of such a system has
been explicitly discussed by Bellman and Astrom (1970) and Johnson (1976).

A system described by (5.1) is controllable if and only if the column
vectors of the n x nm test matrix

= [B: AB: AB: 2 nel
Po = [Bs ABs ATB; ...3 AD BJ (5.8)
span the n-dimensional space on which x(t) is defined. In other words, P.
must be of rank n if the system is to be controllable (Kalman, 1963).

Definition 5.4

A system is open (closed) when there is some (no) exchange of material
from the system either to the environment or to another subsystem. Using the
properties of the controllability matrix Poe Johnson (1976) shows that closed

systems are uncontrollable regardless of their structure. Open systems are
also examined and are found to be controllable, in general, if the flow matrix
A is of rank n or (n-1) and at least one flow between the system and its
environment is controlled. Necessary conditions for complete controllability
can also be derived from physical reasoning. Cobelli and Romanin-Jacur (1976)
adopt a structural point of view by relating controllability and observability
of compartmental systems to their experimental design. A system is shown to
be completely controllable only if there exists at least one path from some
input to every compartment. In addition, it should be noted that the speci-
fication of bounds on the controlled elements implies that a system otherwise
completely controllable must now have a restricted range of controllability in
state space.

The concept of observability of a dynamic system is associated with the
processing of data obtained from measurements on the system. Analogous to the

27


SECTION 6
CONSTANT INFUSION

A constant infusion input has occasionally been used in tracer methods to
label compartments to a high specific activity. The constant influx of tracer
was probably first incorporated by Hevesey and Hahn (1940) to obtain high con-

centrations of 32>_phosphate in plasma. Hevesey and Hahn, however, were not
interested in analyzing the data to determine model parameters. Constant in-
fusion is also used when there are problems with incomplete and non-instantan-
eous mixing within compartments, usually due to spatial homogeneities. Ship-
ley and Clark (1972) present textbook cases in which the method and its usage
are discussed. Assuming that the chronic infusion of tracer does not disturb
the behavior of the recipient compartments, then the constant input level of
Table 1 gives rise to a steady-state defined by Xe

The constant infusion method as discussed in Section 3 is based upon the
compartment model defined by (3.1) as

x = Ax + Bu. (6.1)

Linearly independent inputs generate an ensemble of linearly independent
steady-state samples of the system state x from which the coefficients of the
matrix of rate coefficients may be determined.

When in steady-state, the system of (6.1) can be described by

QO = Ax + Bu. (6.2)

Recall that x is an n-vector of compartmental states, u is an m-vector of
specified inputs, and B is an n x m matrix where

bij =44

n
and >> by, = 1.
moO

Let Ups Ups sss Uy Span RTM Assuming that the compartmental system is open

so that al exists, then Xy> Xoo res Xp Span at most an m-dimensional subspace
of R” where
=], .... mM. (6.3)

Ax = -Buy k

36


SECTION 3
SAMPLING THEORY AND APPLICATIONS

MODEL IDENTIFICATION

The determination of a model relating structure and function of a system
from experimental time series observations of the state is known as system
identification. It includes the estimation of system parameters for the ob-
servable components within the biological system. One of the fundamental
problems of system identification is whether or not the parameters can be
determined uniquely (Lee, 1964). If unique estimates cannot be obtained (the
system is not identifiable), then the experimental procedures should be al-
tered. Determination of structure from function results from an input-output
tracer experiment. The different identification schemes that are available
can be classified according to the basic elements of the problem and the input
signals, the earliest methods being those based on frequency, step and impulse
responses. Special inputs rather than those of normal operation are necessary
for these schemes, so the techniques are termed off-line identification and
therefore apply only to linear stationary processes where input/output rela-
tions which hold for one set of inputs are good for all inputs (Graupe, 1972).

For the purposes of tracer analysis the local model is a linear donor
controlled biological compartment model defined by

x = Ax + Bu (3.1)

where the state x is an n-dimensional vector of the compartmental contents,
the coefficient matrix A = (a; 5) 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 = (b55) is of n rows and m columns
with entry bij equal to one if input uy enters compartment i, if not, the

entry is zero.

In equation (3.1), the inputs u and their entry into the ecosystem com-
partments, the B matrix, are known. The solution of the identification prob-
lem for (3.1) requires a method to determine the entries in A, containing the

rate coefficients, given sufficient measurements of the state x.
BASIC EXPERIMENT

The use of radioactive tracers is prone to varying experimental proce-
dures, especially with the advent of highly sophisticated instruments and new
radioisotopes. Practical experiments must be designed to coexist with devel-
oped theoretical techniques. A standard method has developed using three
basic experiments to assess the environmental impact on aquatic microcosms of


u(t)
ary ——_——_—
1 H
t i)
4a | \
1 ‘
Us a ‘——$<$— ‘
' 1
—eon 1
1 t
uy h | cen |
| 1 " L 1 —_
0 yy T, 13 Ty Tye
Figure 10. Diagnostic input for experimental design.
x(t)
XxX toe -+
! .
! ; ‘
14
Xeaa { Poe po
i419 a
rotopoy te
oot ya te!
a a >
Ti 4 | | T; t
Figure 11. Response to diagnostic input for experimental design.

46


sufficient, but not necessary, conditions for the structural identifiability
of the system when direct measurements can be obtained from every compartment.
There is, however, no need to investigate the observability of the system as
this concept plays no role in solving the problems of constant infusion iden-
tifiability.

Although the conclusions on identifiability for the constant infusion
method are somewhat limited by system structure and the lack of necessary con-
ditions, it should be pointed out that the prevalent notion that "“identifica-
tion of dynamic parameters on the basis of measurements is only possible if
the measurements are taken when the system is in a transient state" (Graupe,
1976, p. 7) is erroneous, because complete sets of parameters can be identi-

fied when a system is in steady-state under conditions which are often met in
practice.

44


in validation of the assumed nonlinear model structure.

To complete the nonlinear analysis, it should be noted that the N (local)
linear models which result from the diagnostic staircase input can be used to
compute the material flows between each compartment as they depend upon the
compartment levels. Such tables of global measurements of ecological struc-
ture and function have been used by Smith (1970) to identify the model struc-
ture for an ecosystem. Thus, it may be possible to use the local models to
validate the nonlinear model structure, and by doing so provide an important
connection between the linear and nonlinear model identification methods.

A pre-experiment to determine the parameters of the Mulholland-Gowdy sam-
pling theory has been run at the Environmental Research Laboratory in Athens,
Georgia, for a microcosm consisting of water, algae, snails and sediment com-
partments. A slug of radiolabelled methoxychlor was introduced into the water
of the laboratory ecosystem. Samples of the toxin and radioactivity were
obtained in order to determine estimates for degradation and compartmental
transfers. These estimates define a and 6, the maximum turnover and minimum
exogenous turnover. By analyzing the methoxychlor samples, some idea of the
arithmetic precision of the data is obtained. All of these estimates for the
sampling theory parameters determine an experimental design, based upon the
diagnostic input of Figure 10, for the particular microcosm under study. It
is hoped the experimental design will be tested in the laboratory sometime in
the future.

50


sta’ = ts ay! (8.2)

where S is a similarity transformation having the eigenvalues of % as its T
column vectors. Equation (8.2) is solved successively for the columns of A’,
again using Gaussian elimination. Note that the column vectors must be line

early independent for the existence of S” -l_ This follows from the assumption
that the eigenvalues of 6 are distinct.

It should be mentioned that EIVECT derives the complex matrix A for a
real matrix @ using the complex principal value natural logarithm in (8.1).
Also, all manipulations are in double precision. The input format for N, K
and the two flags is I4, T is F10.0, and the state data is 8F10.0. The output
is available in full double precision, but the current version of the program
truncates to seven digits as given by 8D14.7 format.

The code, as submitted to the project officer, is well documented by
comment cards, and the particular deck submitted will include data input cards
for Smith's (1970) ecosystem model. An output from this deck, as run on the
Oklahoma State University computer is included in the user package. A com-
plete listing of the code which implements the identification algorithm is
also in the user package.

The user package, as described, fulfills the grant obligations. However,
because the interactive version of the program was of use to the researchers
during the grant period, a listing of this code is also included. This list-
ing indicates the changes which have been made to the batch program in order
to implement it as an interactive program. Depending on the computer opera-
ting system, it is possible to load the batch program version as a data set
and then edit it to conform to the interactive mode. This was the procedure
used successfully by the researchers at Oklahoma State University.

54


443+ 433 = 0, (6.27)

which along with equations (6.22) and (6.23) and one equation drawn from
(6.24) defines six independent equations in the six unknown entires in A.

If the compartment(s) without an exogenous input is (are) not closed, it
may be difficult, if not impossible, to identify the compartment model from
steady-state measurements alone. In such cases, it may be useful to define
the maximum identifiable sub-model of the system. The identifiable sub-model
parameter matrix Anax consists of all transfers associated only with compart-

ments which are either closed or have exogenous inputs. The dimension of the
sub-model is (mtc) where c is equal to the number of compartments without an
input or a flow to the environment. For example, consider the system of Fig-
ure 9 with an additional flow from compartment 3 to the environment. The
matrix of rate, coefficients for maximum identifiable sub-model is

a
aA = | 1 . (6.28)

The identifiability of this sub-model depends only on the measurement of
steady-states. In order to identify the full system model an experimental
design which incorporates material flow measurements must be implemented.

As a practical matter, it is often true that the closure of a compartment
indicates its inaccessability with respect to exogenous inputs. Thus, the
case when (mtc) = n is important in the applications of the theory of struc-
tural identifiability. For this case the following theorem gives sufficient
conditions for identifiability.

Theorem 6.4
An open compartment model with m inputs is identifiable provided:

1. The pair (A, B) is controllable.
2. A maximum of (m-1) endogenous inputs enter each of the n compart-

ments.
3. The (n-m) compartments without exogenous inputs are closed.
Proof:
Because (A, B) is controllable, the row rank of xt is m, and
T = =
xX a = zy k=1, 2, ...,TM (6.29)

has a unique solution for the me unknowns in the n-vectors ay Also, the
controllability of (A, B) along with the fact that A is open implies that

T - -
X' aban =O K=1, 25-005 € (6.30)

where c=n-m, gives c(m-1) equations in cm unknowns, These results follow from
the fact that controllability implies an input to (6.1) can produce an arbi-
os n
trary steady-state which is not constrained to lie in a proper subspace of R.
Also, since A is open without closed proper subsystems, the resultant coeffic-

42


TABLE 3. CHI-SQUARED TEST FOR ayy WITH
= 2 =
N = 400 AND 17 DEGREES OF FREEDOM x 1730.05 * 27.59

a2

t i 6, a(chi)
.005 = -16.024 1.39 x 1072 32.4
006 = -15.973 7.90 x 1079 31.4
.007 — -18.964 4.98 x 10°? 36.2
.008 — -16.057 3.43 x 103 35.3
009 -16.039 2.51 x 107? 41.7
.010 --16.010 1.95 x 10° 41.2
.011-16.010 1.57 x 1073 31.
.012-15.983 1.31 x 1079 25.7
.013 -16.022 1.12 x 1073 23.6
1020 -16.016 6.38 x 1077 8.2
.030 -15.992 6.39 x 1074 7.0
040 -16.050 9.89 x 1074 13.7
050 -15.997 1.86 x 1073 25.2
1060 -15.965 3.95 x 1079 10.8
.070  -16.007 9.67 x 107° 15.1
.080 -15.780 2.02 x 107% 20.9
1090 -16.249 2.14 x 107 6.8
.100 -16.190 3.40 x 1072 14.5
110 -16.178 5.57 x 1072 19.7
1120. -16.368 9.33 x 107% 55.7

14


where I" is an (n#2) x (n+2) identity matrix. Thus, the indirect controll-
ability matrix Pie is of rank (n+2) which is a necessary and sufficient con-
dition for controllability.

Two special cases of the indirect controllability matrix are to be inves-

tigated.
Case 1. If B and CB are nonsingular, then the matrix Pi, is of rank 2 +n.

This result follows because ie can be regularly partitioned into a triangular
block matrix

|o

pr, £ a (5.24)
PF 5.
‘ic 0 cB

which is nonsingular if and only if each diagonal box is nonsingluar (Parker
and Eaves, 1960).

Case 2. If C is the n x n accumulative flow matrix of (5.1) which is used in
the indirect algorithm, C is nonsingular because all flows are on the diagon-
al. Now if the matrix A is multiplied on either side by a nonsingular matrix,
the product has the same rank of A. That is, the rank of CA isn. It
follows that the rank of cak isn and if P) is of rank n, then cake, is of
rank n. In theorem (5.1), the matrix cak P, had to be of the rank C. There-

fore, in the special case where there is an accumulation pool for each com-
partment, this condition is automatically satisfied.

Theorem 5.2

For indirect identification, there must be at least two inputs to the
compartmental system.
Proof:

The number of accumulation pools within a system is assumed to be less
than or equal to the number of compartments in the system. That is, 2<n.
In addition, the number of rows in Pi, must be at most equal to the number of

compartments n times the number of inputs m, or
R+n<nm. (5,25)

Now divide both sides of (5.25) by n, yielding

L
Leaecn, (5.26)

but since & <n, then 2 <1 and unless & = 0 (which is the direct method),

m < 2 because there can be no partial inputs. . ; _
~ Measurements of the accumulation pools are to be obtained for identifica-

tion of the system of equations (5.10). In order for the system to be
able in the sense of Kalman (1963), it is necessary and sufficient that the

32


FOREWORD

Environmental protection efforts are increasingly directed towards
preventing adverse health and ecological effects associated with compounds
of natural or human origin. As part of this Laboratory's research on the
occurrence, movement, transformation, impact, and control of environmental
contaminants, the Environmental Systems Branch studies complexes of environ-
mental processes that affect the transport, transformation, degradation,
and impact of pollutants or other materials in soil and water.

This research includes the development and testing of mathematical
models that can predict general behavior of toxic organic materials in
terrestrial and aquatic ecosystems. The sampling theory and related systems
analyses in this report combine model development with preliminary testing
using small laboratory ecosystems. This method represents an explicit
procedure for designing laboratory ecosystem experiments that provide
appropriate tests of model validity.

David W. Duttweiler

Director

Environmental Research Laboratory
Athens, Georgia


SECTION 5
IDENTIFIABILITY

Before implementation of the identification algorithm, conditions should
be investigated that will insure the uniqueness of a numerical estimation of
the associated parameters. A priori knowledge of which and/or how many param-
eters of the flow matrix A can be determined by means of the chosen input-out-
put experiment is termed identifiability (Cobelli, 1976). If the internal
couplings cannot be uniquely determined, changes in the mathematical model
structure or the experiment itself must be modified. Therefore, an identifi-
ability analysis may be considered as a necesseary step preceding identifica-
tion.

Often, identification is derived from numerical values of a few known
parameters and is reduced to a problem of state estimation in which it is
natural to deal with identifiability in a probabilistic manner. Identifiabil-
ity in such stochastic systems may be defined by the existence of an estimator
(Tse, 1973, 1978) or by the uniqueness of the probability distribution func-
tion (Bowden, 1973). A survey of similar investigations is provided by Reid
(1977). In the study of biological systems, identification often depends upon
@ priori knowledge of the model structure (at least the order) and it is
appropriate to use a deterministic canonical system approach to the determin-
ation of model parameters. Other approaches to the problem include paramet-
rization (Glover, 1974) and parameter sensitivity analysis (Reid, 1977), which
are both closely related to transfer function methods.

An explicit discussion of the problem of structural identfiability in the
frequency domain is given by Bellman and Astrom (1970) where they propose the
use of the transfer function to investigate identifiability. Consider the
system

X(t)
y(t)
which is the linear compartment model defined by (3.1) with a single dose

input u(t), where y(t) is an m-vector of measurements of the state x and C is
ap Xn connectivity matrix of the form

uw
{>
Ix

(t) + B u(t) (5.1a)
(t) (5.1b)

"
ike)
[><

Ci. = (5.2)

and


sentation of the sequence of transition states.
. The stochastic process which represents the transition states has been
simulated by a deterministic variable (z,) perturbed by a uniformly distrib-

uted random variable (wy)

y= Zy + We (3.15)

for each k = 0, 1, ... , mn. The motivation for this representation derives
from knowing (by measurement) the system state to an arithmetic precision of

1074, Then, the random variable We is uniformly distributed from 0 to 104

with a probability density function given by Figure 2. This perturbation
simulates noise superimposed upon the truncated state.

P (a)

108 __

— nel

104
Figure 2. Probability density function for We

The solution of (3.11) for the identified model shows in (3.14) a depen-
dence upon the sample period. Indeed, there probably exists an optimal sample
period which depends upon the statistics of the experiment. In order to test
for this optimal sample period and in general study how noise propagates
through the identification algorithm to the entries in the mode? A-matrix, a
Monte Carlo experiment was decided upon. The exact solution of the model,
correct to d-digits, has been obtained. The state is then truncated to d-
digits and corrupted by noise according to (3.15). Simulated noise measure-
ment data, defined by (3.15) are passed through the identification algorithm
many times until the statistics of the entires in the A-matrix are measure-
able. The following example will clarify this scheme.

The ecosystem model chosen to represent, through computer simulations,
the effect of measurement noise on the identification algorithm and sampling
theory is a hypothetical ecosystem first described by Smith (1970). Smith's
system is not a model of a real system, but is a model built upon realistic
numerical parameters and a state response which mimics an aquatic ecosystem
operating under a phosphorus loading. Very simply, Smith's system consists of
three compartments: water (x1)» an aquatic plant population (x9) and an

herbivore population (x3). The only ecosystem input (uy) is the phosphorus

contained in the water flowing through the system. The outputs are phosphorus
from the water compartment and from the herbivore compartment due to emigra-
tion or emergence. The plants exchange phosphorus with the water compartment,

im


Zt

10 T T TOIT T T TTT TTY] T

O8F

06

>
Oo

Pw

047

02+

1 or Dui iil L

fe) fl
OO! .003 Ol O03 l
Tv

Figure 4. Spectral graph of the estimated variance for ayy of Smith's model.


the mean and variance for Wye It was then assumed that 400 computer runs of

the identification algorithm were sufficient.
For each sample period, the random number Wy, was added to the simulated

states 2, to form the measured data Xy for the identification algorithm. A
total of 400 models were generated for each t. Using the estimates given by

N
ua = (1/N) x a; (3.17)

and

, N
Gy? = (1/N-1) 2% (a; - u,)? (3.18)

for N = 400, the mean and variance for the aj, and aj, entries in the A-matrix
of (3.16) were determined.

A chi-squared fitness test was performed to test for normality of the
distributions for ay and ayo: The results are given in Tables 3 and 4.

Spectral plots of fi, and 6,” for aj, and a), appear in Figures 3 through
6. For this particular model ecosystem, the empirical data suggest bounds of

0.007 and 0.07 for the sample period t as compared tot = 5 x 1078 and 1* =

0.02 as computed from the theory for d = 5. Bounds for t cannot be specified
more closely because of the statistical nature of the data. For example, the
estimated mean for ayy is always within acceptable limits of the true value of

a F -16 (see Figure 3) regardless of the sample period and independent of

the estimated variance. Therefore, the error inherent within the identifica-
tion algorithm cannot be exactly specified. However, the spectral plots in-
dicate a relationship exists between the sample variance and the sample per-
jod. It is possible that future research efforts will disclose a statistical
theory for this relationship.

Finally, a comment on the apparent disagreement between the theoretical
sample period bounds and those derived empirically from this specific example
seems to be in order. The theoretical bounds given in (3.4) and (3.6) are
general bounds for all compartmental systems with maximum turnover of a/2 = 95
and minimum non-zero exogenous turnover of 6 = 2, and that these bounds are
exceeded for a specific model does not invalidate the theory. Indeed, a

sample period chosen less than = 0.02, say t = 0.01, appears to be perfect-
ly acceptable for model identification.

13


for all i=1,

«++ 5» P. Let Y(s) be the Laplace transfi f
that of u(t). Then 4s) P nsform of y(t) and U(s)

Y(s) = G(s) U(s) (5.3)
where G(s) is the transfer matrix given by
G(s) = C(st- ay! (5.4)

and I is the n-th order identity matrix. The problem of identifying the flow
matrix A from the experimental data may be reduced to that of identifying the
coefficients of the transfer matrix, each coefficient being a combination of
the elements of A. Therefore, the A matrix is identifiable when this non-
linear system of equations has a unique solution.

Lin and Yu (1977) considered the single input and single output special
case of (5.1) in canonical form, where A is the companion matrix associated
with the characteristic polynomial ~

m m-1
Sota, s +t... tajstay = 0 (5.5)
and
Cc = [cqc, wee cyl (5.6)
B) = [00... 1. (5.7)

This leads to a scalar input/outout formulation for (5.3) in which G(s) con-
tains 2n unknown parameters (n poles, n-1 zeros, and a gain factor). The
identification of these parameters is demonstrated by a time series analysis
of the output variable given by
Ayt Ae
y(t) = ye ae

A numerical algorithm called Prony's method (Hildebrand, 1974) is applied to
y(t) in order to resolve the parameters 4G and Y% of the exponential curve

fit. The G and Y; parameters are then shown to uniquely determine the 2n

unknown parameters of G(s).

The precise definition of identifiability provided by Bellman and Astrom
(1970) is equivalent to saying that the parameters of A are identifiable if
and only if any change in the values of the coefficients causes a change in
the measured quantities y(t). Thus, a necessary condition for identifiability
is that a system be controllable and observable in the sense of Kalman (1963).

Definition 5.1
A linear system is said to be completely state controllable if for any

initial state at ty it can be driven to any other state in the state space at
te=T.

26


Because the eigenvalues of the system described by (3.1) are bounded by
Gerschgorin circles of known radii (Hearon, 1963), it is not difficult to show
that transitions between steady-states, as in Table 1, are bounded by two ex-
ponential factors (see Mulholland and Gowdy, 1977):

exp(-at) < x, (t) < exp(-st) (3.2)

where a/2 = 1, is the maximum turnover and 6 is the minimum non-zero trans-

fer out of the system. The selection of the sample period t for data collec-
tion is now possible. Using (3.2), it is easily shown that t is bounded by

TST < vv, where v is the maximum sample period and t is the minimum,

both of which are determined by the arithmetic precision of the state samples.
It is possible for the sample period to be too small. This occurs when

x; (t) and x, (ttt) have the same arithmetic precision; that is, when x, (t+)

is truncated to d digits, it equals the truncated value of x, (t). Since (3.2)
gives an upper bound for x, (t), a measure of the minimum sampling time, de-
noted by t , is given by

1-104 = exp(-6x7) <x, (7) (3.3)

or solving for the minimum sampling period yields

w= (-1/8) In (1-074). (3.4)
If t >t the state samples are independent, but an upper bound on t
exists. Suppose t is chosen large enough so that x, (nt), the last sample in

the data set, is truncated to zero, then the data set is incomplete and the
identification problem cannot be solved. Thus, the maximum sample period,

denoted by v, must be such that the (n+1)-th sample is at least equal to the
smallest number in the finite precision arithmetic; that, x, (nt*) > 104, But

(3.2) gives a lower bound for x, (t)s thus, the maximum sample period is easily
computed from

1074

exp(-ant”) < x,(nt"), (3.5)

which gives
+
T

(-172nt,7!) in (1074). (3.6)

Based upon intuition, experience, isolated laboratory experiments, and/or
estimates from the literature, biologists tend to know the turnovers for the
compartments present in a system. Such knowledge is at least accurate to an
order of magnitude. Thus, the maximum sample period is easily determined from
(3.6).


EIVECT )

READ:
Tv
PRINT:
2
| i
CALL Ait di
EIGENP A Tat?
PRINT: PRINT:
Wo Len A
Tij wit LeN
|
dj Cin w)/e

Figure 12. Flow chart for computer algorithm.

52


ACKNOWLEDGMENTS

It is a pleasure to acknowledge the constant support and sustained inter-
est in this research program provided by Dr. William R. Emanuel of the Envi-
ronmental Sciences Division of Oak Ridge National Laboratory. Specifically.
Dr. Emanuel provided the FORTRAN code for the subroutine EISPACK which proved
to be extremely useful for eigenvalue/eigenvector calculations.

It is often difficult to distinguish between basic and applied research.
And, with respect to the research described in this report, it is often hard
to discern the part funded by the EPA grant from that funded as basic research
in systems theory by NSF Grant ENG75-05341-A01.

viii


Odb

-40db

-60db

| ee

T —T

20 Log |; (f) |

Figure 8.

30 100
f(1/days)

Example Bode plot.

23

300 1000


SECTION 4
FREQUENCY DOMAIN TECHNIQUES

; For both ecosystem modeling and experimental data analysis, the tech-
niques of frequency domain analysis, using the Fourier series and transform,
are becoming more popular with ecologists, e.g., see Platt and Denman (1975),
Waide et al. (1974), Shugart (1978), Emanuel et al. (1978), and Van Voris et
al. (1978). This is intriguing because the common method of analysis for
determining the sample period for a communication system involves the frequen-
cy domain concept of bandwidth, expressed in terms of the Nyquist theory.
Therefore, it would be useful to have a system bandwidth defined for the sam-
pling and quantizing precesses used in the solution of the identification
problem.

A transfer function can be defined for linear systems with one input u(t)
and one output x(t). When the single input equals the unit impulse, also
called the Dirac delta function 6(t), the output is called the impulse re-
sponse and is denoted by h(t). The transfer function for the system is then
defined as the Fourier integral transform of the impulse response:

H(f) = of h(t) a MEae (4.1)
where H(f) is the transfer function, f = w/2m is the frequency independent

variable measured in reciprocal time units and i = v-T.
A sequence of n transfer functions for an ecosystem can be defined by

choosing, in turn, each compartment x (j = 1, ... , n) of the state vector x
as an output. The n transfer functions are then calculated by using
= 0 -iwt
H(f) = 2? x; (t)eTM"at (4.2)

for all j=l, ...,n.

Transfer functions reflect the dynamics of a system through the concept
of bandwidth, loosely defined as the range of frequencies which the system
described by H(f) passes in a relatively unchanged form. Carlson (1975) notes
that the rate at which a system can vary stored energy is determined by the
frequency response to perturbations as measured by the bandwidth. Thus, for
ecosystems, bandwidth and compartmental turnover appear to be related.

The so-called bandlimited systems have transfer functions which satisfy

H(f) = 0, for all f > B. (4.3)
For such systems the bandwidth is easily defined to be equal to B. However,

most systems, including ecosystems, are not bandlimited, so other methods of
defining the bandwidth are employed. One such method, called by Carlson

20


CONTENTS

Foreword . 2. 2. ee ee Vii
Abstract . 2. 1 ee ee ee ee ee iv
Figures 2... ee ee ee ee vi
Tables 2... . ee eee ee vii
Acknowledgments ..........0..-05008. viii
1. Introduction ..........2..0.-. 1
2. Conclusions and Recommendations ..... 3
3. Sampling Theory and Applications .... 4
4. Frequency Domain Techniques ....... 20
5. Identifiability...........0.0- 25
6. Constant Infusion... .......04. 36
7. Experimental Design... 1... 2... 0 45
8. Computer Program .........20- 51
References... 1 ee ee eee ee ee ee 55


ABSTRACT

This research considered the application of systems analysis to the study
of laboratory ecosystems. The work concerned the development of a methodology
which was shown to be useful in the design of laboratory experiments, the pro-
cessing and interpretation of the results of these experiments, the develop-
ment of model structures for microcosms based upon experimental data, the
testing of the predictive capabilities of evaluative models, and the analysis
of measurement errors in the presence of noise. Development of the experi-
mental design methodology was based upon a non-Nyquist sampling theory which
gave a priori computable bounds on the sample period for data collection.

Results of the research project included a definition for the bandwidth
of an ecosystem in the frequency domain, an indication of the effect of mea-
surement noise on the sample period for data acquisition, a scheme for reduc-
ing this noise, and a computer code for model identification based upon the
experimental design developed. The project also produced the conditions under
which two new schemes for model identification were applicable to the con-
struction of linear compartment models. These new identification schemes were
based upon system input perturbations and state observations which were often
readily implemented in the laboratory.

The research objective was to produce a scheme for the evaluation of
microcosm models used for the assessment of the environmental transport,
transformation, degradation, and fate of toxic materials in aquatic systems.
An hypothesis was advanced that the evaluative scheme will be developed in the
near future into a procedure by which models are tested against laboratory
data and refined until they realistically represent perturbed natural eco-
systems receiving toxic inputs.

The report was submitted in fulfillment of Research Grant R-805564010 by
Oklahoma State University under the sponsorship of the U.S. Environmental Pro-
tection Agency. This report covers the period October 23, 1977, to October
24, 1978, and the work was completed as of November 1, 1978.

iv


IDENTIFICATION ALGORITHM

Consider the experiment of part (d) of Table 1 in which a constant
(known) input is applied to an ecosystem modeled by (3.1). After sufficient
time has elapsed, a steady-state is observed, as indicated by x = 0. Denote
this steady-state by the vector Xoo which is the solution of

O = Ax, + Bu. (3.7)
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, the exper-
iment is re-started at time t = 0 by removing all inputs to the ecosystem.
Thus, before t = 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 “s

m= Ax (3.8)

because u = 0 for t > 0. The labeled contents of the compartments of the

ecosystem start at t = 0 with values determined by X, and decrease exponen-

tially as t becomes large as all the labeled material flushes from the system
along the pathways for material transfer.
Discretized in time the solution of (3.8) is governed by

x (it) = & x[(i-1)r] (i= 1,2, ...), (3.9)

where & = exp(At) with tT <t<t.
A continuous record of the state x(t) is not required in order to deter-
mine &, which by definition gives A. Since ® is an n x n constant matrix, con-

straint equations numbering ne are required to prescribe the entires 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, T, 2T, ... , kt
where t is the sample period. Substitution of these state samples into (3.9)
gives

x(t) = & x(0)
x(2t) = @ x(t)
x(kt) = & x[(k-1)t]. (3.10)

This linear algebraic system contains n2 unknowns (all the entries of @) and k
independent vector equations of n components each, fora total of kn equations.
For k < n no unique solution of (3.10) exists, and the system is overdeter-
mined for k > n. Thus k = n or (nt1) discrete samples of the state are neces-
sary and sufficient to determine 6.


TABLES

Number Page
1 Basic Microcosm Experiments... 6. 1 ee ee eee ee ee ee 6
2 Truncated State Samples for t=0.01 ........25.+02- 12

3 Chi-Squared Test for any with N = 400 and 17 Degrees of
2 =
Freedom x 1730.05 = 27.59 ee 14

4 Chi-Squared Test for arp with N = 400 and 17 Degrees of
2 =
Freedom y’ 1730.05 7 27.59 ee ee eee 15

vii


TABLE 4. CHI-SQUARED TEST FOR aj WITH
N = 400 AND 17 DEGREES OF FREEDOM 17:0 05 = 27-59

T i, 5, a(chi)
005 5.147 2.655 27.2
-006 4.793 1.194 23.3
007 4.541 0.627 20.6
.008 5.578 0.372 19.2
-009 5.404 0.241 21.2
-010 5.035 0.168 18.2
-011 5.184 0.124 16.5
-012 4.888 0.096 16.9
-013 5.165 0.077 13.9
-020 5.113 0.033 5.7
-030 4.947 1.028 20.8
-040 5.322 0.040 14.3
-050 4.984 0.076 22.8
-060 4.786 0.145 14.4
-070 5.046 0.345 12.5
080 3.680 0.701 19.5
-090 6.516 0.781 9.1
+100 6.160 1.262 12.9
-110 6.081 2.096 7.7
120 7.306 3.485 45.1

15


United States Environmental Research EPA.600 3-78-104
Environmental Protection Laboratory December 1978
Agency Athens GA 30605

SEPA New Sampling
Theory for
Measuring
Ecosystem Structure

KM Sy


partments as long as they do not form a decoupled subsystem. The A-matrix for
an open system is characterized by a nonzero determinant. This result implies

that in the special case when the inputs wy (k=1, ..., 0) span Rk, the resul-

tant steady-states X (k=1, ..., n) also span R" and xl exists. Hence, the
matrix equation (6.5) can be solved by multiplying on the right by xt, giving

A= Wk V Input access to every compartment in the system and the ability to
measure all compartmental steady-states then imply structural identifiability.
These sufficient conditions are only rarely met in practice, requiring for
biochemical systems the organic synthesis of tracer material for all chemical
species present in the system. For ecological systems, only a few natural
input pathways to the system exist with many compartments being essentially
inaccessable. Therefore, the remainder of this chapter concerns the develop-
ment of conditions which relate identifiability to system structure with re-
stricted compartmental access.

Theorem 6.1

A compartmental system is identifiable only if there exists at least one
path from. some constant input to every compartment.
Proof:

By contradiction, if there is no path from any input to one or more com-
partments, then for all k, equation (6.3) can be partitioned as

Ay Ag | | Seq Wy
= (6.13)
2 As} L ke 2
Hence, the following relationships arise:
Ay Xq * Ap Xo = Wk (6.14)

But Ay! exists, that is, there are no closed subsystems, and the vector x, is

identically null. Thus, regardless of the values of the elements comprising
Ag, the model parameters contained in A, and Ay cannot be identified by sol-

ving equations of the form (6.10). Thus, structural identifiability depends
on the ability of inputs to reach every compartment.

Theorem 6.2

A system is identifiable by the constant infusion method only if at least
one outlet to the environment is reachable from every compartment along any
path.

Proof:
By contradiction, if no system outlets can be reached from compartment m,

then this compartment and all compartments coupled to it continuously accumu-
late material without reaching a steady-state. This, of course, prohibits

identification using (6.5).

38


experiments provided a rigorous theory for sampling ecosystem structure is
employed for data acquisition. Within the arithmetic precision of the experi-
mental data acquired from the microcosm, such a theory for sampling ecosystem
structure must completely capture the inherent dynamics of the mathematical
model. Thus, the use of system analysis is employed to help insure the
results of microcosm experiments are useful in testing the prediction of
ecological fate and effect made by mathematical models.

The idea that systems analysis has application to microcosm research is
not new. System applications usually focus on the use of models to explore
potentially interesting microcosm experiments, the testing of models against
microcosm data, and the use of basic model structure to define what variables
are of interest and hence what should be measured. The last potential use
mentioned may hint at the most important application for system analysis, i.e.
its use in the design of experiments. Mulholland and Gowdy (1977) presented,
for perhaps the first time, the concept that systems analysis and modeling can
be used to develop a quantitative theory for measuring ecological structure in
laboratory microcosm. The idea is to use systems analysis, in much the same
way statistics is now used, to design microcosm experiments. Starting with a
theory used to produce bounds on the temporal sampling period for measuring
compartment level variables in a microcosm, research during the past year has
produced a design for input function which places in evidence the complete
dynamic response resulting from a system level perturbation. The resultant
experimental design provides for the collection of sufficient data to identify
the parameters of two models types: (1) a sequence of (local) linear models
with a priori unconstrained structure and (2) a nonlinear model of postulated,
but not quantified, structure. The basic input function design also allows
for the collection of independent data useful in testing the responses of the
resultant models to perturbations.

In summary, research has disclosed an experimental design for laboratory
ecosystems which provides a method for introducing the appropriate system
level perturbations, introduced through natural system input ports, in order
to elicit full system dynamic response. The design also provides for ampli-
tude calibration and periodic time sampling which captures the full dynamic
state response. The state response data which result from the experimental
design establish a test for system linearity with respect to multiple pertur-
bations, a data base for constructing local linear models, an independent data
base for nonlinear model construction, and still another data base for testing
model response.


SECTION 1
INTRODUCTION

This research report concerns an investigation into the application of
systems analysis to the study of laboratory ecosystems. Despite the recent
resurgence of interest in microcosms, there appears to be no clear understand-
ing of how to best use system theory and mathematical modeling in this impor-
tant research area. However, there is apparently some agreement as to the use
of laboratory ecosystems in ecological research. At least with respect to
general research objectives important to the U.S. Environmental Protection
Agency, the following reasons for microcosm research are mentioned by Draggan
(1976), Brockway, et al. (1978) and Van Voris, et al. (1978): (1) to develop
a screening method capable of assessing the environmental fate and effect of
chemicals, (2) to study basic ecosystem processes under controlled environmen-
tal conditions, (3) to study ecosystem level variables as measures for ecolog-
ical complexity, and (4) to test or verify mathematical models.

The appeal for microcosm studies is presumably based upon the replicabil-
ity of experiments which can be performed. This replicable nature of micro-
cosm experiments is generally assumed to result from the control which can be
exerted over the laboratory environment. Thus, an ecologist concerned with
the study of ecosystem function can employ the hypothesis that microcosms sub-
jected to similar physical and chemical conditions, even when inoculated with
samples derived from different sources, will develop into functionally similar
systems, see Brockway, et al. (1978). Two quantitative tools are required for
work based upon this hypothesis of functional invariance: (1) a rigorous ex-
perimental design for measuring ecosystem structure and (2) a theoretical con-
struct which relates ecosystem structure and function. The research reported
on here presents the idea that system theory and mathematical modeling can
provide both of these tools.

There are also diverse opinions on how to interpret or apply the results
of microcosm research. Many ecologists view microcosms as “stripped down" or
“simplified" biological models which are intended to simulate some specific
natural process (Odum, 1971, p. 21). On the other hand, Brockway, et al.
(1978) view laboratory ecosystems neither as biological models nor as analogs
of natural ecosystems. Instead, they take what may be a more realistic view-
point that microcosms are ecosystems in their own right which have their own
community structure, mineral cycling, and pathways for energy flow typical of
all natural systems. It can be expected that toxic material would disrupt
such integrated microcosm processes in a manner similar to any other ecosystem
of related structure and function. Thus, a mathematical model which is de-
signed to predict the fate and effect of a substance in any ecosystem which
represents a class of systems defined by specific functional attributes should
at least be capable of similar predictions for a microcosm belonging to the
same class. It follows that models developed for screening can be tested in
the laboratory against replicable results derived from perturbed microcosm

1


SECTION 8
COMPUTER PROGRAM

A digital computer user package implementing the compartment model iden-
tification algorithm has been developed. Programmed in standard FORTRAN, the
code which implements the identification algorithm is available in two forms,
batch and interactive, both of which run efficiently on a standard IBM 370
operating system. The code is independent of any special library routines,
and as such is completely transportable to any IBM System/360 or System/370
series computer. The interactive program can be implemented on any IBM
System/370 time share operating system (TSO).

The software consists of a main program and one subroutine, EIVECT, the
flow diagrams of which are shown in Figures 12 and 13.

The program begins by reading the number N of compartments in the bio-
logical system being analyzed, the number K of measurements obtained from a
single experiment, and the pre-calculated sample period T. Two flags are then
set which indicate row or column data input and printing options, respective-
ly. The program then moves immediately to a loop that reads and stores the
sample data array in matrix form, i.e., X(t). Since N+] discrete measurements
of the state are necessary and sufficient for the determination of the state
transition matrix , the program terminates for K < N+]. If K = N+1, two
matrices of coefficients, X and Xos are formed from the N+] state samples as

described by (3.12) and (3.13). The solution for the state transition matrix
@ is acquired by inverting matrix X, of (3.12). Gaussian elimination with

partial pivoting on a matrix scale is used. Upon computing the N x N matrix
subroutine EIVECT is called for the identification of the matrix of rate
coefficients A. The indices I and L of the main flow diagram are utilized in
the formation of a loop in which several A matrices may be computed for a
comparative analysis if K > N+].

Subroutine EIVECT solves (3.14) for the matrix of transfer constraints A.
The EISPACK (IBM/360 version by Smith, et al., 1976) is utilized to obtain the
eigenvalues wy (i = 1, ..., N) and the eigenvectors of ®. The EISPACK uses
QR transformation and decomposition, an explanation of which is given by Wil-
kinson (1965) and Wilkinson and Reinsch (1971).

The matrix can be reduced to a diagonal matrix @ in which the eigen-

values 0; (i= Ty > N) of & appear on the diagonal. The eigenvalues ;
(i = 1, ..., N) of A are determined by
d~ = (1/T) a, a> (8.1)

1

and the reduced diagonal form A of the matrix of transfers A is formed. The
matrix of rate coefficients is readily identified by

5]


; os Chr alcls .2c3 aT)! 1 ch tos uly ately 2.5 caTynm2gt
i 7 : . 5.31
“di 1 : 0 : 0 (5.31)

is of the rank p + 2 where I is an &-th order identity matrix. Define the
upper partition to be

Pay = [Os Pos Pel (5.32)

and notice that Pe is the observability matrix La for the direct algorithm
while Pe is the observability matrix Pig for the indirect scheme for identifi-

cation. The system is completely observable if the conditions for direct
observability are satisfied because Pe would be of rank p, making the system

efrio P
Pay = 6 (5.33)
0 1 0

of rank p-where I' and I are p-th order and £-th order identity matrices,
respectively. Also, if there are as many (or more) accumulation pools as
directly measureable states (p), and the conditions for indirect observa-
bility hold, then the system is of the form

Io

ol

Pai (5.34)

Io dh?

and the desired rank is obtained. In addition, if the two test matrices
(direct and indirect) combine to form a matrix P, of rank p, the system
becomes a

Io

P. I' 0
Pag 7 7) E I> = (5.35)
Io ool

and the necessary and sufficient conditions for identifiability are fulfilled.
Consider next a system which has as its output the sum of observations.
In this combined case, the measurement vector for the system is

z= fC I |x (5.36)
y

and the system is said to be observable in the sense of Kalman (1963) if and
only if the test matrix

po = cols alcls 2.3 (ATI 'CT] + Cos xUs aTKTs «5 (AT)TM?K] (5.37)

is of rank p. Thus, the system is observable if the sum of the observability
matrices for the direct and indirect algorithms is of rank p. However, two

34


can be used to reconstruct continuous time functions of system responses. The
two inequalities derived are incommensurate, however assuming 6 = 2B as a rule
of thumb, the Nyquist result implies that ecosystems should be sampled with a
sample period less than or equal to the reciprocal of the minimum exogenous
turnover. Another rule of thumb is obtained by assuming equality in (4.8) and
H,(0)a = 1 in (4.6) giving 2B/a = 1. Since a/2 = Ths the bandwidth of an
ecosystem is then equal to the maximum compartmental turnover.

For systems which are not bandlimited the bandwidth concept is at least
elusive. The bandwidth of an ecosystem has not been clearly defined, there-
fore it may be useful to consider the results presented in the light of an
example. As an elementary example of the frequency domain techniques pre-
sented, consider again Smith's hypothetical ecosystem defined by (3.16).

The transfer function H.(f) defines the frequency response of compartment

j with respect to the input which enters compartment one. The Bode plot for
the amplitude response is obtained from

Logarithmic gain = 20 Togy gIH5 (FI (4.9)

which is measured in units of decibels (db). Plots of (4.9) for the three
compartmental transfer functions of Smith's model are given in Figure 8. These
Bode plots of the system transfer functions, each of which contains several
poles and zeros, are obtained by adding the individual plots due to each pole
and zero. Furthermore, this procedure is simplified by using only the asymp-
totic approximations to these curves, as shown in Figure 8. These Bode plots
describe the frequency response of the example ecosystem.

From (3.16) simple calculations yield 6 = 2 and a = 190. Therefore,
using t = Yt*t~ and (3.6) and (3.4), the sample period is t = 0.00078 day for
d = 4, which can be rounded to t = 0.001 day.

As for the ecosystem bandwidth, the inequality 1 < 2B/S gives B >
1(1/day). Figure 8 indicates the bandwidth to be approximately equal to 100
(1/day). Equating T= Ts results in B < 500 (1/day). Thus, the inequality

1 < 2B/8 yields very little useful information regarding the bandwidth of this
ecosystem, but the rule of thumb, derived as 2B/a = 1, gives B = 95 (1/day)
which seems reasonable.

Consider Figure 8 and measure the logarithmic gain at f = 95 (1/day) for
each compartmental transfer function. Note that |H, (95) | > |Hp (95) | >

|H3(95) |. which implies B, > B, > B, as an ordering of the bandwidths for the

transfer functions of the three compartments. This ordering in turn implies

Ty < Ty < T3 for the three compartmental turnover times where T. = 1/a...

The ordering of the turnover times given is ecologically reasonable. In this
open (flow through) system, the phosphorus dynamics of the water compartment

are fast, more rapid than the uptake and transfer of phosphorus by the plants.
Also, the dynamics of phosphorus in the plants should be faster than for ani-
mals, as indicated. All this tends to add ecological intuition to the inter-
pretation of the mathematics of determining the ecosystem bandwidth using B =

To

22


The preceding concepts of structural identifiability based upon the data
derived from constant infusion experiments parallel those already known for
the single dose method (see Section 5). However, the solution of the identi-
fication problem of a constant input applied to a system is further con-
strained by the static nature of the resultant steady-state. Indeed, as will
be shown, the dynamic systems concept of observability plays no role in sol-
ving this algebraic problem.

Theorem 6.3

For m exogenous inputs u to an identifiable compartmental system, there

can be at most (m-1) endogenous inputs to any single compartment k.
Proof:

For m inputs into a system, equation (6.10) produces m independent equa-
tions in the n unknown parameters of aye But at most m unknowns can be deter-

mined from the solution of m independent equations. Recall that the k-th
entry in the n-vector ay represents the turnover of the k-th compartment. If

the compartment turnover is to be estimated, then at most m-1 of the endogen-
ous inputs to the k-th compartment (the remaining entries in ay) can be iden-

tified. Thus, for m inputs to a system, there can be at most m-1 endogenous
inputs to a single compartment.

The (n-m) compartments r without exogenous inputs produce m dependent
equations of the form

x'a, = 0 (6.16)

where a. is defined as in equation (6.10). Therefore, (n-m) additional equa~
tions are needed for the solution of (6.5). The necessary equations may be
obtained from each of the (n-m) compartments a, in one of the following ways:
1. Turnover of x, must be known from an independent measurement, or from
a priori knowldege of the compartmental turnover, giving the r-th
entry in aye

2. Flow out of the compartment to the environment must be determined for
an open compartment.

3. A closed compartment offers an additional independent equation given
by the nul] column sum

n
> ai, = 0 (6.17)
iF

where aap (i=1, ..., n) define the components of a.

Consider the three compartment system of Figure 9 which gives an A-matrix
of unknown system parameters of the form

39


(1975) the effective bandwidth, uses the equal area concept. Consider Figure
7 which plots and compares the amplitude of a typical transfer function (the
solid curve) with that of a rectangular window (the dashed curve). The rec-
tangular window is clearly bandlimited with bandwidth of B. The effective
bandwidth for H(f) is equal to B under the constraint:

2B | H(0) | = _s° | H(f) | af (4.4)

which amounts to finding B to give equal area under the two curves.

| H(F) |

roc rote s-crrrra4

i

!

{

; +

~B 0 B
Figure 7. Effective bandwidth construction.

It is now possible to compute the bandwidth for the ecosystem transfer

functions defined by (4.2). First, recall that x,(t) is amplitude limited by
(3.2). Next, let f = 0 in (4.2), giving 4

H;(0) = of x; (t)dt. (4.5)
Substituting (3.2) into (4.5), gives the following inequality

Va < H;(0) < 16 (4.6)
which holds for all j = 1, ... , n. In order to complete the computation of B

from (4.4), it is necessary to evaluate the integral which appears in this
equation. Estimating this integral is possible through the definition of the
inverse Fourier transform:

x(t) = f° Hy (feltar, (4.7)
Taking absolute values of (4.7) for all t > 0 and j = 1, ... , n, yields
Tos Lf” | Hy(F) | af. (4.8)

It is not possible to fully combine the inequalities (4.6) and (4.7) into
(4.4) to estimate the bandwidth, the best result being 1 < 2B/6. Thus, the
bandwidth of an ecosystem depends upon turnover, but not in the ordinary
sense. This result implies the bandwidth B to be greater than or equal to
one-half the minimum exogenous turnover, that is the transfer rate of material
out of the system to the environment.

For bandlimited frequency domain representations, the Nyquist sampling
theorem gives TS B/2, where Tw, is the sample period from which measured data

21


é

TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)

1, REPORT NO. 2. 3. RECIPIENT’S ACCESSION NO.
EPA-600/3-78-104
4. TITLE AND SUBTITLE 5. REPORT DATE

New Sampling Theory for Measuring Ecosystem Structure December 1978 issuing date

Is. PERFORMING ORGANIZATION CODE

7. AUTHOR(S) 8. PERFORMING ORGANIZATION REPORT NO,
Robert J. Mulholland
I9. PERFORMING ORGANIZATION NAME AND ADDRESS 10, PROGRAM ELEMENT NO.
1LA760

School of Electrical Engineering
Oklahoma State University

71. CONTRACT/GRANT NO.

Stillwater, OK 74074 R-805564010
12, SPONSORING AGENCY NAME AND ADDRESS 13. TYPE OF REPORT AND PERIOD COVERED
Environmental Research Laboratory - Athens, GA Final, 10/77-10/78
Office of Research and Development 14. SPONSORING AGENCY CODE
U.S. Environmental Protection Agency EPA/600/01

Athens, GA 30605

15, SUPPLEMENTARY NOTES.

16. ABSTRACT

This research considered the application of systems analysis to the study
of laboratory ecosystems. The work concerned the development of a methodology
which was shown to be useful in the design of laboratory experiments, the pro-
cessing and interpretation of the results of these experiments, the development
of model structures for microcosms based upon experimental data, the testing
of the predictive capabilities of evaluative models, and the analysis of measure-
ment errors in the presence of noise. Development of the experimental design
methodology was based upon a non-Nyquist sampling theory which gave a priort
computable bounds on the sample period for data collection.

Results of the research project included a definition for the bandwidth
of an ecosystem in the frequency domain, an indication of the effect of measure-
ment noise on the sample period for data acquisition, a scheme for reducing
this noise, and a computer code for model identification based upon the experi-
mental design developed. The project also produced the conditions under which
two new schemes for model identification were applicable to the construction
of linear compartment models. These new identification schemes were based
upon system input perturbations and state observations which were often readily
implemented in the laboratory.

17. KEY WORDS AND DOCUMENT ANALYSIS.
a. DESCRIPTORS. Ib. IDENTIFIERS/OPEN ENOED TERMS |c. COSATI Field/Group
Systems analysis
Mathematical models 62B
Environments 68D
72E
78. DISTRIBUTION STATEMENT 19. SECURITY CLASS (This Report) 21, NO, OF PAGES
Unclassified 66
Release to Public 20. SECURITY CLASS (This page) 22. PRICE
Unclassified

EPA Form 2220-1 (9-73) 58 ‘US. GOVERRMENT PRINTING OFFICE. 1979 ~657-060/1556



Number

a
12
13

FIGURES

Page.

An example two compartment system with precursor ......... 7
Probability density function for Wee eee W
Spectral graph of the estimated mean value for ayy of Smith's

model 2. ee 16
Spectral graph of the estimated variance for any of Smith's

model 2... ee 7
Spectral graph of the estimated mean value for alo of Smith's

model 2 wwe ee ee ee ee 18
Spectral graph of the estimated variance for a2 of Smith's

Model 2 ww we 19
Effective bandwidth construction... 2... 1 ee eee eee 21
Example Bode plot . 2. 1. ee ee te ee ee ee ee 23
Hypothetical compartment system .. 1... ee ee ee ee 40
Diagnostic input for experimental design... 2... 2. ee eee 46
Response to diagnostic input for experimental design. ...... 46
Flow chart for computer algorithm ... 1... ee ee eee eee 52
Flow chart for subroutine EIVECT .. 2... 2 ee ee ee ee ee 53

vi


which is transferred to the herbivores through grazing.

_ From data presented by Smith, the rate coefficients can be computed,
which results in

x F 5x, + 5x3 - 16x, + uy

14x, - 95x, (3.16)

u

90x, - 14x

*3 2 3
as a (local) linear donor controlled model defined near the prescribed input
of uy = 100 mgP/day. The time scale of the model is defined by the rate co-

efficients which are given in terms of units of reciprocal days.
With uy = 100, equation (3.16) was simulated on the IBM-370 using the

CSMP simulation language package. State samples correct to 4 1/2 decimal
places, i.e. Xo takes initial values with 5 significant digits but it eventu-

ally declines to 4 signigicant digits, were obtained for 20 different sample
periods ranging between 0.005 and 0.120. It should be noted that (3.16) with

a = 190 and 6 = 2 has sample period bounds of t = 5 x 1076 and t = 0.0202

ford = 5, and t =5 x 1079 and t* = 0.016 ford =4. Fort = 0.01, the four
state samples are given in Table 2. These samples illustrate the particular
basic experiment performed where the initial steady-state entrained by the
constant input is the initial condition. The input is then removed and the
state declines, according to the system dynamics, monotonically to zero.

TABLE 2. TRUNCATED STATE SAMPLES FOR t = 0.01

k x] Xo X3

0 9.5000 1.4000 9.0000
1 8.5748 1.3506 8.9844
2 7.7807 1.2549 8.9052
3 7.0939 1.1519 8.7513

A pseudo-random number generator was then used to drive the Monte Carlo
simulations for model identification. The density function for We uniformly

distributed from 0 to 1074, gives a mean y= 5 x 107° and a variance of
ae = 8.33 x 19710, Using standard statistical tests with a level of signif-
icance of 0.05, it was determined that 400 samples were required to estimate

12


system, which results in a steady-state x being observed. Then another con-
stant input Uy # uy is applied, and its corresponding steady-state Xo is also

observed. Substituting these steady-states into (7.13), while noting that a
steady-state is defined by the time derivative of x being identically zero,
two simultaneous algebraic equations result:

0

-ax, + bx, + u

2
to (7.14)
2
2

+ Uy.

Because the inputs and states are measured, and hence known, equation (7.14)
represents two linear algebraic equations in the unknowns a and b. Further-
More, a unique solution of (7.14) for a and b exists if and only if Xx # Xp

Thus, a model of known structure can be identified from steady-state measure-
ments. However, since an infinite number of nonlinear models could fit the
steady-state data, either confidence in the model structure must be present or
some validation scheme for the model is required.

To demonstrate the multiplicity of nonlinear models which can be identi-
fied from a fixed steady-state data set, consider the model structure:

x = cx + dxk tu (7.15)

0 = -aX, + bx

where c and d are unknown parameters and k > 3 is any integer. The same
steady-state data set used for identifying (7.13) now gives

0 = “CX, + dx tu

1
(7.16)

Nw ae

0 = “CX_ + dx. + Uy

which for each fixed k > 3 yields a unique solution as long as x) ? Xp. Some

knowledge of the ecological process under study and its theoretical formula-
tion could help to decide which model structure best represents the system.
Also, each model structure could be compared with time series data for the
transient response of the system in order to decide the best model fit to the
data.

The scheme presented for nonlinear identification of the parameters of
models of known structure generalizes. For (7.13) the model consists of a
single equation with two unknown parameters which are obtained from two steady
state measurements. In general, a model consisting of n compartments and M
unknown parameters requires N steady-states measurements, where

nN = M. (7.17)

Thus, in Figure 10, the N levels in the diagnostic input provide a data set
capable of determining M = nN unknown parameters in a nonlinear model of known
structure which represents a microcosm of n compartments. At time Tye a

large (positive) step input is again applied to the system to re-establish the
initial state. The time series which results from this step input is of use

49


Hevessey, G. and L. Hahn. 1940. Turnover of Lecithin, Cephalin and Sphingo-
myelin. Biol. Meddr. 15:5.

Hildebrand, F. B. 1956. Introduction to Numerical Analysis. New York:
McGraw-Hil1.

Johnson, L. E. 1976. Control and Controllability in Compartmental Systems.
Math. Biosci. 30:181.

Kalman, R. E. 1963. Mathematical Description of Linear Dynamical Systems.
SIAM J. Control. 1:153.

Lee, R. C. K. 1964. Optimal Estimation, Identification, and Control. Cam-
bridge, MA: MIT Press.

Lin, P. L. and Y. P. Yu. 1977. Identification of Linear Time-Invariant
Systems Using Exponential Curve Fitting. Proc. IEEE. 65:1508.

Luenberger, D. G. 1964. Observing the State of a Linear System. IEEE Trans.
Mil. Elec. MIL-8:74.

Mulholland, R. J. and C. M. Gowdy. 1977. Theory and Application of the Mea-
surement of Structure and Determination of Function for Laboratory Eco-
systems. J. Theor, Biol. 69:321.

Odum, E. P. 1971. Fundamentals of Ecology. Third edition. Philadelphia:
Saunders.

Ogata, K. 1967. State Space Analysis of Control Systems. Engelwood Cliffs,
NJ: Prentice-Hall.

Platt, T. and K. L. Denman. 1975. Spectral Analysis in Ecology. Ann. Rev.
Ecol. System. 6:189.

Reid, J. G. 1977. Structural Identifiability in Linear Time-Invariant Sys-
tems. IEEE Trans. Autom, Control. AC-22:242.

Shipley, R. A. and R. E. Clark. 1972. Tracer Method for In Vivo Kinetics.
New York: Academic Press.

Shugart, H. H., ed. 1978. Time Series and Ecological Processes. Philadel-
phia: SIAM Publ.

Smith, F. E. 1970. Analysis of Ecosystems. In; Analysis of Temperate For-
est Ecosystems, D. E. Reichle, ed. New York: Springer-Verlag, pp. 7-18.

Smith, W. D. and R. R. Mohler. 1976. Necessary and Sufficient Condition in
the Tracer Determination of Compartmental System Order. J. Theor. Biol.
57:1.

Thron, C. D. 1972. Structure and Kinetic Behavior of Linear Multicompartment
Systems. Bull. Math. Biophys. 34:277.

56


INDIRECT METHOD

The necessary conditions given for controllability and observability of a
system for the direct identification algorithm imply that any model of a real
system with one or more compartments acting as stores for the system may not
be identifiable. With the indirect scheme, the observer states y forma set
of flows disjoint from the internal couplings, each accumulation compartment
being a singular trap. In this section, criteria are investigated for con-
trollability and observability related to compartmental structure of a system
having & accumulation pools (@ <n) from which all observations are to be
made. The local model of (3.1) is again incorporated along with the linear
differential observer for the accumulation pools:

y= KK (5.10a)
which can be combined to give
x AO} {x B
= + u (5.10b)
v KO} ty o

from which the total system behavior may be explored. Now define the measure-
ment vector z as follows:

z= fo J

I<

(5.10c)
y

Measurements are to be taken from each of the n observer states. The system
described by (5.10) is controllable if and only if the partitioned matrix

Py 2 3 qn-l

{oo
>
a
>

[w
(>
[oo

vat2y (5.11)

sic P

Py 0

5
iE
5

spans the (2 + n) state space on which the system is defined, where the two
n x nm submatrices Py and Po are called the upper and lower partitions, re-

spectively, of the controllability matrix. Now if Pi. is of rank 2n, then P,
and Po must both be of full rank. Note that the upper partition of the con-

trollability matrix is the same test matrix as that for the direct method (see
equation (5.8)). Thus, as with the direct method, if there is no path from

any input to compartment i, then the i-th row of the controllability matrix is
identically null, as no input enters compartment i. Next consider the lower
partition of the matrix Page If, in addition to the aforementioned condition,

there is no accumulation compartment that can be reached from a given compart-
ment i, then every element of the i-th row of the controllability matrix is

29


null as no pool is influenced by the state variable associated with the i-th
compartment.

For a clearer understanding of the observer portion of the controllabili-
ty matrix, consider the dynamic states y to be continuous observations rather
than accumulation pools. Then from Kalman (1963) and the definition of out-
put controllability, the system is said to be completely output controllable
if and only if the composite matrix

= . - KAeR: . yan-lp.
Poe «= ~(EKBs KABs KATBs ...3 KA” Bs

|o

J (5.12)

is of the same rank as the row order of K, where D is the connectivity matrix
which postulates a direct transmission from u to y. The existence of a direct
input into any accumulation pool is not assumed to exist, that is, D=0. In
fact, if the same compartment is available for insertion and observation of
tracer, identification requires that the associated flow matrix be irreducible
(Smith, 1976), whereas the flow matrix for the indirect identification scheme
is reducible. Thus, tracer must be injected into the internal system x and
not into the observer y. If the rank of the 2 x nm matrix

Py = [Os KBs KABs KA’Bs ...3 KA'?B] (5.13)

is of rank 2, then the rank of the g x nm matrix Poe given by (5.12) is neces-

sarily 2 Notice that state controllability is neither necessary nor suffi-
cient for output controllability.

The discussion of indirect controllability to this point has relied upon
the knowledge of the rank of Pao which is not customarily known when discuss-

ing identifiability. Conditions must be placed upon the submatrices Py and Po
such that the indirect controllability matrix will be of rank 2 +n, giving a
necessary condition for identifiability.

Definition 5.6

A system is said to be reduced controllable if and only if the rank of
the matrix Pe where

= [B; ABs AB; ...; ak"1B] (5.14)

Py

is n. The smallest possible integer k for which Py is of rank n will be

called the reduced order. Sufficient conditions for the controllability
matrix P.. to be of the desired rank follow.

Definition 5.7

Two m x n matrices are defined to be equivalent if and only if one may be
transformed into the other by nonsingluar elementary row and column transfor-
mations. If A and B are equivalent, we write

AEB (5.15)

30


SECTION 7
EXPERIMENTAL DESIGN

The experimental design for microcosms receiving toxic material inputs
developed by this research project focuses upon deciding what can be measured
and what models support the measured data. The basic design requires pre-
experiment measurements in order to determine the parameters of the Mulhol land-
Gowdy sampling theory for system identification and modeling. Some idea of
the time scale of the system must be obtained for a useful definition of the
flows and material transfers between compartments. Of course, the number of
compartments in the system should be discerned as a first step.

The Mulholland-Gowdy sampling theory requires estimates for the maximum
compartmental turnover and the minimum exogenous transfer rate. The maximum
turnover describes how fast the toxic material moves from one compartment to
another, while the minimum exogenous turnover bounds the toxicant degradation
rate for each compartment. Pre-experiment measurements are also needed to
establish the precision of the measured data, e.g., the number of significant
digits in the data output of the laboratory measuring equipment. The periodic
sampling interval t for time series state data collection resulting from a
ests experiment, as described by Table 1(e), is then defined by (3.4) and
3.6).

Many diagnostic inputs, of the form of basic experiments (Table 1), are
required to perturb the laboratory system in order to elicit the transient
state response required for data collection and system identification. Thus,
the basis for the experimental design is a diagnostic input of the form shown
in Figure 10. It is assumed the laboratory system is of a flow-through type,
open to its environment, so that a constant inflow of toxic material Uy

applied to the natural microcosm inputs establishes a steady-state. This
initial system steady-state is observed by noting the constant levels of
toxin within each compartment. The inflow Uy is a chronic level of toxin

carefully chosen to have minimal effects on the biota. The time required to
establish the initial steady-state X, response to uy is approximately equal to

six (6) time constants, obtained by multiplying the reciprocal of the minimum
exogenous turnover by six (6). Assuming the flow-through process takes place
at a rate greater than degradation, the settling time estimate produces 6/6,
which essentially accounts for the exponential behavior of degradation in the
transient state. At time Ty the input uy is reduced by Au, where

uy = Uy > Au, (7.1)

so that the microcosm seeks a new steady-state Xx}. In order to identify a

45


gL

7.0

6.0r

T | T | T TT

40

3.0

peo i n 1 pau sriiil i 1 i

.OO|

Figure 5.

"003 Ol 03 1 3
T

Spectral graph of the estimated mean value for Alp of Smith's model.


to insure the system is in steady-state before the next input transition takes
place.

The Mulholland-Gowdy sampling theory has been used to design an experi-
ment for studying laboratory ecosystems. Experimental design parameters (a,
6, d, etc.) can all be obtained from pre-experiments involving the microcosm
or its isolated components. The result of an experiment based upon the diag-
nostic input of Figure 10 is a data set which is capable of constructing a
series of (local) linear models for the microcosm. Furthermore, the validity
of these linear models can be tested using the measured data. Two tests exist
for model validation. First, for linear models, the input level ratios must
equal their respective steady-state response level ratios:

Yi4q/Uy = X54 /%y (7.12)

for all i = 0, 1, ..., N.. Of course, practical laboratory experiments will
result in statistical estimates for these steady-states, and standard statis-
tical tests can be used to determine the validity of (7.12). The second test
is qualitative, noting that linear compartment models produce monotone tran-
sition states in response to basic experiments in which transfers between
steady-states are observed. If the samples of the microcosm state during any
(Ty4y - T,) time period do not decrease monotonically from steady-state Xj to

XGa> then the linear compartment modeling hypothesis is incorrect. Probably

all microcosms will fail this test for linearity when Au is too large. Thus,
local behavior of the system and the ability to describe this behavior by
linear models are believed to be related.

The validity of the experimental design, as based upon the linear com-
partment modeling objective, depends upon the equality of the steady-state
ratios and the monotonic behavior of the transition between steady-states.
When these validation tests fail, the microcosm modeling objective must be
based upon nonlinear system theory.

If a great deal of intuition exists regarding the microcosm processes,
then it may be possible to characterize these processes mathematically within
a nonlinear model of known structure. In terms of the general concepts of
systems theory, it is extremely difficult to identify nonlinear models using
transient state data. The technique presented in Section 3 of this report
applies only to linear systems and because it relies on a general solution for
the linear model, the scheme cannot be generalized to nonlinear identifica-
tion. However, another approach for obtaining the parameters of a nonlinear
model of known structure is possible using only steady-state measurements.

The proposed method of nonlinear steady-state model identification is
best illustrated by example. Consider the very elementary single compartment
system with input u and state x. From pre-experiment observations, scientific
intuition, or theoretical constructs regarding the microcosm processes, the
structure for the model is assumed to be of the form

kos eax t bx? tu (7.13)

where a and b are unknown parameters which characterize the model structure to
the particular ecological processes. For b = 0 the system is linear and the
parameter a is the compartmental turnover. The identification of the param-
eters a and b proceeds as follows. Let a constant input uy be applied to the

48


Tse, E. 1973. Information Matrix and Local Identifiability of Parameters.
In: Proc. Joint Autom. Control Conf., Columbus, Ohio. pp. 611-619.

Tse, E. 1978. A Quantitative Measure of Identifiability. IEEE Trans. Sys.
Man Cyber. SMC-8:1.

Van Voris, P., R. V. O'Neill, H. H. Shugart, and W. R. Emanuel. 1978.
Functional Complexity and Ecosystem Stability: An Experimental Approach.
ORNL/TM-6199, Oak Ridge National Laboratory, Oak Ridge, TN.

Waide, J. B. 1974. A Linear Systems Analysis of the Calcium Cycle in a
Forested Watershed Ecosystem. In: Progress in Theoretical Biology,
R. Rosen and F. M. Shell, ed. New York: Academic Press, pp. 261-345.

Wilkinson, J. H, 1963, The Algebraic Eigenvalue Problem. London: Oxford
University Press.

Wilkinson, J. H. and C. Reinsch. 1971. Linear Algebra. New York: Springer-
Verlag.

57


DISCLAIMER

This report has been reviewed by the Environmental Research Laboratory,
U.S. Environmental Protection Agency, Athens, Goergia, and approved for
publication. Approval does not signify that the contents necessarily reflect
the views and policies of the U.S. Environmental Protection Agency. nor does

mention of trade names or commercial products constitute endorsement or
recommendation for use.

ii


Now set

~Buy = Ww k=l, ..., m (6.4)
which results in
AX = W (6.5)
where the n x m matrices
X= [xys Xpo sees Xd (6.6)
and
Wo= [Wy wos e+ Wa] (6.7)

are both of rank m. Taking the transpose of each side of equation (6.5) gives
T AT T

XA = W. (6.8)
Next define the n-vector a, to be the k-th column of al such that
al = [a,, a a] (6.9)
A Ay» Apr sees ay :

The vectors ay define the unknown system parameters, i.e., the rate coeffic-
jents of the compartment model. Thus,

T = =
x a = zy k=1, ..., 0 (6.10)
where 2 is the m-vector of the k-th column of
ze Ww (6.11)
where
Zo= [zys Zps wes 2] (6.12)

and rank (Z) =m. It should be noted that the vectors ay and zy of (6.10) are
completely determined by the compartmental inputs. That is, ay represents the
k-th row of A giving the endogenous inputs to the compartment and zy defines

the possibility of an exogenous input to compartment k.

Identifiability related to the system structure for a multi-input/multi-
output system where all states can be measured directly is to be considered
first. The system is structurally identifiable when the number of parameters
that can be estimated is not smaller than the number of known system parame-
ters in the algebraic equation (6.5).

Only compartment models which are open to their environment are consid-
ered in this chapter because closed systems with constant rate inputs accumu-
late material at a constant rate, never reaching steady-state. This condition
does not prohibit the closure of an individual compartment or groups of com-

37


matrices of different orders cannot be added, so in order for the test matrix
Pe to be defined, the number of direct observations, p, must equal the number

of accumulation pools, &.

35


toxic material inputs. As shown in Table 1, the first test input represents a
slug of toxic material, the second considers a constant chronic input level of
toxin, and the third forms the superposition of the aforementioned experi-
ments. The microcosm responses to the three basic inputs are also shown in
Table 1. The toxic slug of part (a) elicits a transient response with initial
condition Xj,» and the constant input level of part (b) gives rise to a steady-

state defined by Xe. By assuming a linear model for the microcosm system,

part (c) of Table 1 illustrates the resultant superposition principle. It
should be noted that the basic microcosm experiments contain information re-
garding the validity of the linear modeling assumption. For linear models,
the microcosm responses to the inputs of Table 1, parts (a) and (b), should
sum to give the response to the input of part (c). If this superposition is
in fact true, the single dose method [the input of part (a)] can be effective-
ly used to identify the numerical model for the microcosm. Therefore, the
validity of the linear state mode? (3.1) should be clear from the outcomes of
the basic experiments of Table 1. Generally, a local linear model can be
identified for small perturbations (xy is small compared to x) caused by im-

pulse inputs at a given system loading defined by Uy: A local model defined

in this way will be dependent upon the system loading.

Also in Table 1, two new experiments are proposed for linear model iden-
tification using tracers or toxic material flows as microcosm inputs. These
inputs, based upon step functions, have been designed for ease of application
in the laboratory and mathematical convenience for system identification.
These basic experimental inputs lead to a sampling theory which is essential
to the successful application of standard identification schemes. This re-
sult, in turn, prescribes an experimental design for the a priori determin-
ation of the period for data collection in the laboratory.

The basic experiment of part (e) of Table 1 has been designed as an iden-
tification technique for local linear models. Linearity can be tested in the
laboratory by observing response superposition as result of adding two micro-
cosm inputs.

As previously noted, the basic experiments (d) and (e) of Table 1 have
been designed for ease of application in the laboratory and to aid in the de-
velopment of a sampling theory. It turns out that under certain conditions
these basic experiments also avoid the problem of tracer analysis of having
input access to each compartment. However, the equivalence between the inputs
of parts (c) and (e) of Table 1 is essential in eliminating the so-called pre-
cursor problem. In Figure 1 an example of a two compartment model with a
single dose input into the precursor compartment is shown in part (a) along
with the input and compartment response plots in part (b).

The precursor problem deals with the impossibility of estimating the max-
imum value (x) of the response of compartment Xy to the impulse input to com-

partment XxX: In order to develop an effective laboratory data collection

strategy, it is necessary to have known rules for amplitude and time discret-

ization. Because Xn is known before the experiments is completed, it is very

difficult in ordinary tracer analysis to develop a rule for amplitude quanti-
zation. Given a sample period, data on microcosm structure is collected at
discrete times and the values of these data are quantized by the finite arith-
metic precision of the measuring instruments. However, the quantization of

5


linear model for the laboratory ecosystem, n samples of the transient state
must be resolved between the steady-states Xo and xy: Thus, the input change

Au must be large enough to produce this amplitude resolution. For linear com-
partment models, equation (3.2) prescribes that the transient state x(t) is
bounded by

x(t) > x, + (xj-xy Je" (7.2)
Linear system theory also gives
(ug/uy) = (xg/x,)s (7.3)
and since
x, =X, > aK (7.4)

it is easily shown that (7.1), (7.3) and (7.4) combine to give
Ax = (xo/U) Au. (7.5)
Assuming a sampling period of t, it is clear that x(nt), the n-th sample taken

after Ty. must be different than Xp> j.e.,

d

x(nt) > xt 107 (7.6)

within the finite arithmetic precision of the measurements. Thus, equation
(7.2) gives

axe"T 5 7974 (7.7)
and a lower bound for Au results by substituting (7.5) into (7.7):
-ont 4.-d
Au > (u,/Xq) e 10. (7.8)

This process is continued through time Ty resulting in enough measured data to
construct N (local) linear models for ecosystem structure and function. Since
only positive inflows of toxic material are possible

Ndu <u, (7.9)
must hold true. Hence, an upper bound on Au is obtained
Au < u /N. (7.10)
In response to the staircase diagnostic input of Figure 10, at least n

samples of the microcosm state are measured with a sampling interval of t. As
in Figure 11, the time difference Ti - 1) during which the transient state

response is observed must be greater than six (6) time constants, i.e.,

Ti Ty 2 6/8 (7.11)

i+]

47


I I¢l
pe tithe

X(ID), X(T 1), X (ITN)
i

® X2x7

Figure 13. Flow chart for subroutine EIVECT.

53


(2) eer

impulse transfer output
al
input °| = %4 %
(vo) Ax, XQ
Xm
x,
° '
:
H
—» ¢ >
T t T t

Figure 1. An example two compartment system with precursor.

time and amplitude are not sequential processes; they are, in fact, related to
one another. This gives rise to the commonly used laboratory procedure for
single dose experiments of taking microcosm samples very rapidly at first
until Xo reaches its maximum and then more slowly to accommodate the compart-

mental dynamics inherent in the free response. The asynchronous sampling
procedure described is difficult to analyze, particularly when the propagation
of measurement errors through the identification algorithm is required. A
periodic sampling scheme would be more desirable.

The basic experiments (d) and (e) of Table 1 eliminate the precursor
problem. This follows from the fact that transitions between steady-states in
compartmental system are monotonic and bounded by the initial and final states
(Thron, 1972). Therefore, the amplitude quantization process is well defined,
and a periodic sampling strategy is possible.

SAMPLING THEORY

The collection of time series data for each compartment x, (t) in the bio-

logical system takes place in two distinct steps in which first time and then
amplitude are quantized. The samples of the state are obtained only at dis-
crete times. These samples are periodically spaced apart by t time units
giving a data base defined by {x(0), x(t), ... , x[(kt+1)t]}. It has been
shown that the minimum data base for solution of the identification problem is
(n+1) state samples or k = n. Next, the finite precision with which the
Measurements of state are made must be accounted for in terms of a mathemat-
ical relationship. The assumption is that x, (t) and its samples are fixed

point numbers which are less than or equal 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 1074, where d is a
positive integer, be the smallest non-zero number in the data set.


NEW SAMPLING THEORY FOR
MEASURING ECOSYSTEM STRUCTURE

by

Robert J. Mulholland
School of Electrical Engineering
Oklahoma State University
Stillwater, Oklahoma 74074

Grant No. R-805564010

Project Officer

dames Hill, IV
Environmental Systems Branch
Environmental Research Laboratory
Athens, Georgia 30605

ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATHENS, GEORGIA 30605

EPA-600/3-78-104
December 1978


test for controllability, Kalman (1963) has shown that a linear, time-invar-
iant system defined by (5.1) is observable if and only if the n x np matrix
Py o= tcls abel; (ay? chs 22.5 catye'ely (5.9)
is of rank n. If there is only one output (p = 1) then a necessary and
sufficient condition for observability is that Py be nonsingular. The test

for observability may be applied to compartmental models by utilizing tech-
niques that were employed when discussing controllability. Upon inspection of
the definition of observability, it is clear that observability does not
depend upon the controlled input matrix, but rather on the output matrix. A
linear compartment model, in the structural sense, is completely observable if
at least one output is reachable from every compartment. It is also necessary
that no outputs be associated with only the state components which form one or
more sink, all influencing the same outputs (Cobelli and Romanin-Jacur, 1976).

Definition 5.5

A subsystem of compartments which receives input from the remainder of
the system, but has no transfers out of the subsystem from any of its compart-
ments, is called a trap. A trap may also be called a closed subsystem.

It is important to mention that the conditions given for controllability
and observability are not sufficient, that is, the preceding conditions
neither guarantee the identifiability of A, nor are the elements of A in the
uncontrollable and unobservable part of the system necessarily unidentifiable.
Counter-examples which clarify these points are considered by Distefano
(1977). Bellman and Astrom (1970) emphasize the difficulty of the problem of
identifiability and since that time, only necessary conditions have been re-
lated in the literature. Cobelli and Romanin-Jacur (1976) give criteria for
sufficient conditions, but a close study (Delforge, 1977) shows by counter-
example that the proposed criteria are necessary, but not sufficient. The
problem reduces to a recurring mathematical question which asks how many
solutions a set of nonlinear equations has.

Let A be unidentifiable. Then it may be possible to formulate an experi-
ment to show that a subset of compartmental levels can be controllable even
when the system fails to meet controllability and observability requirements
(Johnson, 1976). For example, consider a system in which certain states are,
for biological reasons, not accessable and their evaluation, according to the
structural identifiability tests, is not possible. If a feasible experimental
condition can be obtained by adding a new input, thereby yielding necessary
conditions, the system may be considered identifiable (Milanese, 1976). Also,
it has been suggested that the simplest way to deal with a multicompartment
system is to ignore closed subsystems if possible because any or all of the
traps can be deleted without altering the behavior of the rest of the system
(Thron, 1972). In general, this is possible with regard to measurements on
compartments of completely open subsystems; and most systems utilizing the
direct method for identification of the A matrix would fall into this cate-
gory. However, certain observations, such as the kinetics of urinary excre-
tion, require the use of an indirect measurement and hence the accumulation
pools may not be excluded. Necessary conditions for identifiability are
therefore derived for the indirect algorithm,

28


9L

-14.0

-15.0

fi, “160 —-————

-170

“T

-18.0

piri!

T

TT TTT

Koje)|

003

0!

Figure 3. Spectral graph of the estimated mean value for ay] of Smith's model.



RESEARCH REPORTING SERIES

Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:

Environmental Health Effects Research

Environmental Protection Technology

Ecological Research

Environmental Monitoring

Socioeconomic Environmental Studies

Scientific and Technical Assessment Reports (STAR)
Interagency Energy-Environment Research and Development
“Special” Reports

Miscellaneous Reports

This report has been assigned to the ECOLOGICAL RESEARCH series. This series
describes research on the effects of pollution on humans, plant and animal spe-
cies, and materials. Problems are assessed for their long- and short-term influ-
ences. Investigations include formation, transport, and pathway studies to deter-
mine the fate of pollutants and their effects. This work provides the technical basis
for setting standards to minimize undesirable changes in living organisms in the
aquatic, terrestrial, and atmospheric environments.

O©MBNOAAONS

This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.


491 apg «0 (6.18)

Figure 9. Hypothetical compartment system.

Assume the two exogenous inputs uy and Uy drive the system to a steady-state,
resulting in the following equation of the form of (6.12):

° iD a EO 1 0] fu,
O} = Jay ao 0 xy {*+ }O 1 uy |. (6.19)
0 0 Azo 39 3 0 60

Performing m=2 experiments with two linearly independent constant inputs gen-
erates the ensemble of steady-states given by (6.5) as

ap 8 813% 2 Uy “42
a, gp 9 Xap Xan} =~] Yor Yoo}. (6.20)
O 3g ag34L%31 X39 o 9
Taking the transpose of each side yields
X11 %o1 X31) [41 ar 9% 41 421 °
X12 X22 %32| {9 822 832 -| 2 "32 9 (6.21)
a3 9 33

40


For systems which are not bandlimited the frequency domain concept of
bandwidth is inexact. Thus, the bandwidth of an ecosystem is not clearly
defined. However, based upon the calculations presented in this paper, B =
qT) >» taken as an approximation, appears to be useful in expressing the eco-
system bandwidth as the maximum compartmental turnover.

The Mu Thol land-Gowdy theory for the periodic sampling of ecosystem
structure gives

TMTTM <T< , (4.10)

where the upper (1*) and lower (t ) bounds on the sample period t depend on
the system turnovers, the arithmetic precision of the data, and the number of
system compartments. Processing of the transition state data, which are of
bandpass character, requires a window in the frequency domain given by

fps fs fy (4.11)

where fy = 1/t” and f, = Ve" describe, respectively, the highest and lowest

frequency (f) present in the data. If an assumption is made that the deter-
ministic system is perturbed by a white noise Gaussian random process, which
contains all frequencies with equally probable magnitudes, then the high fre-
quency noise can be filtered from the experimental measurements by bandlimit-
ing data plus noise to an upper frequency bound given by f = 1/t~. Although
the transition state data is bandpass, the steady-state data is low pass re-
quiring measurements at essentially zero frequency, hence the low frequency
noise cannot be filtered from the data as easily as the high frequencies.
Least squares methods of data analysis remain the best solution to minimizing
low frequency noise. It is important to study the effect of filtering on the
error analysis of the data plus noise which come from the microcosm experi-
ments.
For Smith's hypothetical ecosystem the value for t is 5 x i0~4 assuming
6 =2andd=4. This gives a value of fy = 2000 per day. for which the Bode

plots for all three compartments are well below -60 db. Thus, a low pass
filter with an upper cutoff frequency of 2000 per day would surely pass the
information content of the microcosm experiments while rejecting the high
frequency noise present in the data.

24


test matrix

n-2
P, os KTS ATK; wos aT?
Pio = = (5.27)
Pa U50; 05.45 0

is of rank n + 2 where I' is the 2-th order identity matrix. Note that if the
submatrix Py is of rank n, then

P,£ to (5.28)

where I is the n-th order identity matrix and hence,

P3] e fo tle fr o
Pig = = ; = , (5.29)
Pay fatof jot

giving an equivalent (n+2) x (n+) identity matrix. Thus, for Pio to be of
the desired rank, P, must be of rank n. By considering the states y as con-

tinuous observations of the internal states x, the necessary conditions pre-
sented by Cobelli and Romanin-Jacur (1976) for observability for the direct
algorithm may be applied to the indirect scheme. In the Structural sense, the
indirect model is completely observable if at least one flow to an accumula-
tion compartment is reachable from every internal compartment. In addition,
no accumulative flows may be associated with only the internal state variables
which form one or more traps. all influencing the same pool.

COMBINED DIRECT-INDIRECT

Applications of the direct and indirect algorithms may tend to overlap in
certain systems. An introduced (labeled) substance may metabolize and accumu-
late, along with its metabolites, in one or more pools and also be directly
Measureable as an output from compartments within the system itself. In other
systems, the outputs available for data collection might be the sum of obser-
vations from a compartment and its accumulation pool.

The method in which the observations are made is not involved in control-
lability which is concerned with the flow of the inputs, so the indirect con-
trollability conditions suffice for the combined measurement technique.

Observability of the two types of combined systems will be discussed in
this section. Requirements for observability are straightforward. First con-
sider the system of (5.10) where the measurement vector is of the form

z= . (5.30)
oO ll ty

Measurements are to be collected from both system compartments and accumula~
tion pools. The system is observable if and only if the test matrix

33


6L

3.0 ~~ T 7 TT T TTT T T T T TOT TT T T
2.0 7
a2 L
om |
|.OF 1
8 1 L po rit i 1 n
.OO! .003 Ol 3 J 3
T

Figure 6. Spectral graph of the estimated variance for 12 of Smith's model.


SECTION 2
CONCLUSIONS AND RECOMMENDATIONS

The results of this research have produced an experimental design for the
laboratory analysis and evaluation of perturbed microcosms. The complete de-
sign has not yet been tested against actual experiments. A continuation of
this line of research would allow such a test to take place. Additional re-
search is also needed for the development of a stochastic version of the ex-
perimental design, thus combining statistics with systems analysis which could
lead to a powerful technique for error analysis. The development of a statis-
tical error analysis should give rise to techniques for filtering raw data,
and hence for minimizing measurement errors in the presence of noise. Testing
of mathematical models against microcosm derived data, based upon the proposed
development of a stochastic experimental design, could lead to a rigorous
methodology for model validation.

With respect to analysis in the frequency domain this research has been
unable to reconcile the Mulholland-Gowdy sampling theory with the Nyquist
theory used in engineering systems analysis. The main problem here, yet to be
resolved, concerns the development of an effective technique for bandlimiting
ecosystems. Unlike engineering systems, it may be impossible to define a use-
ful measure for the equivalent bandwidth of an ecosystem. The Mulholland-
Gowdy sampling theory provides bounds on the sample period used for data
collection in the laboratory, which in turn describe a window in the frequency
domain where the state response spectrum resides. However, it is unclear
whether or not this window, having been developed under prescribed input con-
ditions, gives a universal representation for the frequency spectrum of an
ecological system.

The study of two special forms of compartmental structure has produced a
generalization of the concept of identifiability. First, for systems, in
which the state is not directly observable, but measurements of only accumu-
lation pools are possible, the conditions for identifiability, controllability
and observability have been defined. Second, the conditions under which
systems receiving a constant infusion of labeled material are identifiable
have been produced. Although these conditions on the identifiability for the
constant infusion method are somewhat limited by system structure, it should
be noted that most system scientists believe that systems cannot be identified
from only steady-state data and that excitation of transient states are always
required. In spite of this misconception, it has been shown that complete
sets of model parameters can be obtained when systems are in steady-state
under specific conditions which are often met in practice. The study of
identifiability conditions is necessarily theoretical in nature. So, in order
to be of some use in experimental design, the conditions which have resulted
from this research must be tested in the laboratory against the data derived
from microcosm experiments .


REFERENCES

Bellman, R. and K. J. Astrom. 1970. On Structural Identifiability. Math.
Biosci. 7:329.

Bowden: Re 1973. The Theory of Parametric Identification. Econometrics.
2152.

Brockway, D. L., J. Hil] IV, J. R. Maudsley and R. R. Lassiter. 1978. Devel-
opment, Replicability, and Modeling of Naturally Derived Microcosms.
J. Environ. Studies, to appear.

Carlson, A. B. 1975. Communication Systems. Second edition. New York:
McGraw-Hill.

Chen, C. 1970. Introduction to Linear System Theory. New York: Holt, Rine-
hart and Winston.

Cobelli, C. and G. Romanin-Jacur. 1976. On the Structural Identifiability of
Biological Compartment Systems in a General Input-Output Configuration.
Math. Biosci. 30:139.

Delforge, J. 1977. The Problem of Structural Identifiability of a Linear
Compartmental System: Solved or Not? Math. Biosci. 36:119.

Distefano, J. J. 1977. On the Relationships Between Structural Identifiabil-
ity and the Controllability and Observability Properties. IEEE Trans.
Autom. Control. AC-22:652.

Draggan, S. 1976. Role of Microcosms in Ecological Research. Bioscience.
26:402.

Emanuel, W. R., D.C. West and H. H. Shugart. 1978. Spectral Analysis of
Forest Model Time Series. Ecol. Modelling. 4:313.

Glover, K. and J. C. Willems. 1974. Parametrizations of Linear Dynamical
Systems: Canonical Forms and Identifiability. IEEE Trans. Autom.
Control. AC-19:640.

Graupe, D. 1975. Identification of Systems. New York: Robert E. Krieger.

Hearon, J. Z. 1963. Theorems on Linear Systems. Ann. N.Y. Acad. Sci.
108 :36.

55


which is the matrix form of the following three sets of equations in the six
unknown rate coefficients:

1 9 * %31 813 F wt (6.22)
M2 847 + % 32 443 = Ye
for the first column of A,
X11 821 * X21 B22 = a (6.23)
. 6.23
X12 221 * X22 822 = ~Y29
for the second column, and
Xo 832 + X31 833 = 0 (6.24)
a + a = 0

Xoo 932 * X32 933
for the third column.

The system of Figure 9 may be tested for identifiability by inspecting
the necessary conditions provided in this section. There is a path from some
input to every compartment and one outlet to the environment is reachable from
every compartment, thus satisfying Theorems 6.1 and 6.2. The system has two
exogenous inputs and at most one (m-1) endogenous input to a single compart-
ment. Notice that the sets of equations (6.22) and (6.23) have a unique solu-
tion if and only if

x x

x x
VW 31 WW 21 #0. (6.25)

X12 %32 X12 X22

These conditions are required for identifiability and are obtained when the X-
matrix of steady-states has row rank of m=2 and the individual state vectors

can lie anywhere in the m-dimensional subspace of the state space Rr". A suf-
ficient condition for these constraints on the X-matrix of steady-states is
provided by the controllability of (A, B).

Equation (6.22), (6.23) and (6.24) can provide five independent equations
in the six unknown entries in the A-matrix. Thus, with this example, it is
clear that knowledge of the turnover 33 of compartment 3, or the exogenous

flow (a33)X3 from compartment 3 provides the identification of the systems
model. It is also clear that closure of compartment 3 gives

41


ient matrix [see (6.24)] of (6.30) is of rank (m-1). So, equations (6.29) and
(6.30) give

we + c(m-1) = mn-c (6.31)

independent equations in nm unknowns. However, c compartments are closed,
yielding c independent equations for the colums of A. Thus, the solution of
the resultant set of mn linearly independent equations gives the desired nm
rate coefficients comprising the A-matrix.

Now consider the case in which some states cannot be measured directly.
It will be shown by example that it is impossible to identify a1] unknown
parameters unless there is access to every compartment. Assume that steady-
state measurements can be taken only from compartments 1 and 2 of Figure 9.
The ensemble measurement matrix Y may be given by

Y= CX (6.32)

where C is defined by equation (5.2), resulting in

YN. Me 1 0 0 [ a1 Xe
Yoy Yop) = [9 1 0 | Kyo (6.33)
X31 *32

substituting the measured values Y into the system of equation (6.20) gives
the system

My 41% 430 FH
My 821 * Yay 822 Hy
Yo, 894 + X31 833 0 (6.34)
Ne 1 * %2 430 = he
M2 927 * Yq 822 = “Uae

with the additional equation 413 + a33 = 0 provided by equation (6.24). No-

tice that the result is a set of six linearly independent equations with eight
unknowns since the values of X31 and X30 cannot be measured.

Although the parameters associated with compartments which are immeasur-
able cannot be estimated, again, a maximum sub-model can be identified. The
identifiable sub-model flow matrix is composed of those parameters associated
only with those compartments which can be measured directly and satisfy con-
ditions given when all states are measurable. Thus, for Figure 9, the maximum
identifiable sub-model flow matrix Anax is that of (6.25).

Controllability of the system (6.1) allows the exogenous inputs to pro-
duce arbitrary steady-state values anywhere in state space resulting in

43


TABLE 1. BASIC MICROCOSM EXPERIMENTS

Input Response
4 4
Xol---
(a) NL
Ce
T t T t
7 4
Uo Xs
(b)
od a
T t T t
4 4
(c) Xs py
Ll
T t T t
4 4
u Xs1
(a) a
EERE
T a T t
A 4
(ey) Baron Xs2fp——
ee ooo
T t T t