Med. & Biol. Eng. & Comput., 1979, 17, 216-222
Adequacy of measurements in compartmental modelling of metabolic systems R. F. B r o w n *
E.R. Carson
L. Finkelstein
Department of Systems Science, The City University, London EC1V 0HB, England
K.R. Godfrey
R . P . Jones
Inter-University Institute of Engineering Control, Department of Engineering, University of Warwick, Coventry CV4 7AL, England
A b s t r a c t - - T h e adequacy of present-day measurement techniques for the compartmental modelling of metabolic systems is investigated using numerical examples and analysis of experimentallyobtained p/asma clearance curves. It is concluded that the model parameters obtained are often of questionable accuracy. The situation can be improved by careful monitoring of experimental conditions and judicious spacing of data points on the response curve, but the work shows a clear need for a continuous (or semicontinuous) method of monitoring plasma concentration, To resolve ambiguities between models equal/)/ plausible on physiological grounds, it is necessary to monitor the dynamics of internal variables, for example, of the concentration in the fiver (which is nowadays possible noninvasively). K e y w o r d s - - M e t a b o l i c systems, Modelling 1 Introduction
LINEAR multicompartmental analysis of tracer clearance curves is widely used to yield information about the rates of reaction of chemical substances of interest and their transport into relevant body tissues and fluids (JACQUEZ, 1972; G1BALDI and PERRIER, 1975). Clearly, use of a multicompartmental model to reveal impairment of metabolic pathways in diseased states (JONES et al., 1977) or to formulate drug dosage regimens (RANJBARAN and MCLAREN, 1976) implies acceptance of the validity of the model. This paper concerns itself with the acceptance stage--with the adequacy of model structure and of model parameter values for a given structure (CARSON and FINKELSTE1N, 1977; DISTEFANOand CAMPFIELD, 1977). F o r validation to be of practical worth, it must be with respect to the intended use of the model (CARSON and FINKELSTEZN, 1977). Here, the purpose of the model is assumed to be met if the quantitative representation of the important metabolic pathways afforded by the model is in accord with a priori knowledge, in the following ways. First, the model structure (topology) should be compatible with a priori knowledge of the various internal couplings in the physiological system and their significance. Secondly, the internal states of the model (the compartmental concentrations) should be compatible with a priori knowledge of their dynamic responses (often such a priori knowledge comes from * Currently on sabbaticalleave at the University of Warwick from the University of New South Wales, Kensington 2033, NSW, Australia First received 7th Apri/ and in final form 20th AprU 1978
animal studies and its applicability in the human situation is uncertain). Thirdly, some parameters (such as plasma volume) can be estimated from other experiments without the need to postulate a mathematical model; their values are then tabulated in the literature. Predictions of these parameters from the compartmental model are relatively independent of uncertainties in model structure; e.g. plasma volume can be estimated from the zero time value of a plasma clearance curve, and this is relatively independent of the model postulated. The values of these model-independent parameters should be compatible with the a priori knowledge of their values (see Section 2). Inextricably bound up with the problem of model structure validation is the complementary problem of the reliance that may be placed on the parameter estimates. In a typical measurement of a plasma clearance curve, a dose of tracer material is administered intravenously or orally at time zero, and a succession of plasma samples is taken at times beginning several minutes after the administration and ending several hours (but not usually days) after the administration. The concentration xl(t) of tracer in the plasma is then determined for each sample. The form of the compartmental chosen will yield a series of exponentials to be fitted to the data xl(t) = ~
i=1
Ai exp ( - e ~ t ) . . . . .
(1)
the coefficients then being fitted by standard curve. fitting procedures (BEKEY et al., 1974). Over the time scale of the response, it is not
0140-0118/79/020216 -t-07 $01 -50/0 9 IFMBE: 1979
216
Medical & Biological Engineering & Computing
March 1979
mathematically meaningful to attempt to fit more than three exponentials; if a higher-order model is tried, degeneracy will occur in that either some of the coefficients (As) will be negligible or rate constants (es) will be so close as to be mergeable. Thus, with measurement times from several minutes to several hours, the most complex form of concentration curve that can be fitted is x l ( t ) = A~ exp ( - ~ t t ) + A 2
exp (-cz2 t) + A3 exp ( - e3 t) (2)
parameters of the assumed compartmental structure, are also directly measurable. Such parameters are referred to as model-independent parameters, and they serve a particularly important role in assessing the adequacy of any postulated compartmental model. The best-known model-independent parameter is the plasma volume (or more correctly the initial volume of distribution) that may be inferred from the zero-time intercept of the plasma clearance curve in the following manner. Consider an experiment in which labelled tracer is injected into the plasma, and the subsequent concentration x l ( t ) of tracer in the plasma is measured. Curve fitting of the form of eqn. 1 will give the model-predicted initial concentration xl (0). From a knowledge of the quantity Q o! labelled material injected, the plasma volume q is given in normalised form by
Arbitrarily ordering the terms so that el > c(2 > e3, problems of fitting the first ter m (the fastest exponential) often arise because there are insufficient data points before the term has died away to a small value. Problems of fitting the third term (the slowest exponential) arise because the coefficient A3 is usually small; e.g. BERK et al. (1969) found a value of only 3 % (taking A~ + A2 + A3 = 100) in bilirubin q Q clearance tests on a group of 19 normal subjects. m - mxl(O) . . . . . . . . . (3) (It is not always small, and a value of 12% was found in similar tests on a group of thirteen subjects where m is the weight of the subject. If the tracer is with Gilbert's syndrome (BERN et al., 1970)). radioactive, for example, the usual units are Q in The problems discussed are compounded by the counts per minute, xl(0) in counts per minute per presence of various kinds of disturbances that couple millilitre of plasma, q in millilitres and m in kiIointo the measurements. These include various uncon- grammes. trolled factors that may affect the subject's measured The other main model-independent parameters response behaviour, such as bodily movements, are based on the normalised area under the response aural and visual stimuli (which may affect emotional curve when the input is an impulse-function change and hence physical condition) and whether the in infusion rate (where normalisation is with respect concentration of the material of interest (the to the area under the impulse), or the steady-state unlabelled material) is in a steady state at the outset gain from input to output when the input is a stepof the experiment. function change in infusion rate (these are identical Short-term environmental disturbances often for a linear system). If the experiment is concerned manifest themselves in the measurements by values with the kinetics of bilirubin, for example, the commonly referred to as 'outliers', so cal/ed since model-independent parameter of concern is the they lie outside the region of acceptability. Naturally, unconjugated bilirubin production rate (u.b.p.), every effort needs to be made to minimise environ- which varies inversely as the area under the plasma mental disturbances at their source where this is clearance curve. possible, and by choosing the most favourable time In a typical measurement of a plasma clearance and location for the experiment. Long-duration curve, as noted previously, the first plasma sample is environmental disturbances give rise to trends in the data; these trends are often modelled mathematically 100 by a low-order polynomial, usually the constant to term and linear term only. The constant term (bias) is particularly important when establishing the zero reference as, for example, when subtracting the background count in radioactive tracer experiments. 8 10 9 \.\ Set against the background of the problems "'.'x.__ = data time r a n g e t2 previously described, the remainder of this paper will ID be aimed at determining the adequacy of currentlyTx aD available measurements on which compartmental modelling is based, and at indicating how measurements should be improved so as to arrive at models t._ on which greater reliance can be placed. O C 2 Model-independent parameters Certain parameters of tracer response curves, as well as being inferrable from knowledge of the M e d i c a l & Biological Engineering & C o m p u t i n g
0
06
i
i
J
10
20
30
time F,'g. I Logarithmic plot of eqns. 4 and 5
M a r c h 1979
217
taken several minutes after injection of the radioactive tracer and the last sample several hours (but not usually days) after injection. A compartmental model fitted to this data is then used to extrapolate the data to zero time to deduce the steady-state gain of the system. Implicit in this procedure is the assumption that the model structure is valid outside the time range of data used in estimating the model parameters. An important function of experiment design is to perform studies to test the validity of such assumptions. Some numerical examples will now be presented to illustrate likely problems.
is incorrectly estimated to be 9.09 units (true value 18.18 units). E x a m p l e 2: effects o f neglecting a small time delay
Consider again the idealised plasma clearance curve of the preceding example, but with a small time delay inserted into the response so that the curve is defined by x~(t) = O, O ~< t < 0 . 1 = 9.09e-r176
-a~176
t >! 0.1
(6)
If measurement does not begin until after t = 0.1, then eqn. 6 may validly be written
E x a m p l e 1 : effects o f neglecting a fast transient
Fig. 1 shows a logarithmic plot of an analytic idealisation of a plasma clearance curve, defined in normalised units by
x~(t) = l O . 0 5 e - t + 2 4 7 . 1 2 e - l ~
t > 0.1
(7)
If this equation is fitted to the data
(4)
(a) the plasma volume is incorrectly estimated to be 0.39 units (true value 1 unit)
[with xl(0) = 100]. A possible compartmental model which would represent eqn. 4 is shown in Fig. 2a, where the volume of compartment 1 (the plasma compartment) is defined to be ql = 1 unit. Assume that measurements of plasma concentration are made over the time range 0.6 to 3 units and that the curve
(b) the area under the curve is incorrectly estimated to be 34.76 units (true value 18.18 units).
x~(t) = 9 . 0 9 e - t + 9 0 . 9 1 e -a~
x~(t) = 9.09e -t .
.
.
.
. . . .
.
.
.
.
(5)
is fitted to the data. The compartmental model corresponding to eqn. 5 is shown in Fig. 2b.
"(~1.09 0 '73
q1=1
a
9 - 0 9 6 ( t ) ~
1"00
ql =11 b Fig. 2 a Compartmental eqn. 4
model
corresponding
to
b Compartmental eqn. 5
model
corresponding
to
In this example, neglecting the fast transient term 90.91e - l ~ has dramatic effects: (a) the zero-time intercept is 9.09 (eqn. 5) instead of 100 (eqn. 4), so that the plasma volume is, from eqn. 3, incorrectly estimated to be 11 units (true value 1 unit) (b) the area under the curve, which from eqn. 1 is given by ~A~ i=1
218
Q{i
BERK et al. (1969), in a study of bilirubin metabolism in 13 normal subjects, obtained the following mean plasma clearance curve for unconjugated bilirubin: xz(t) = 66e-~176 + 31e-~176176 + 3 e - ~ 1 7 6 1 7 6 (8)
(normalised so that xl(0) = 100). The small-valued slow transient term 3e -~176176 was essential in their modelling in order to obtain a correct area under the curve, and hence to obtain correct values for unconjugated bilirubin production rate and redblood-cell lifespan.
9.18
1006(t) --C)=
Example 3: effeets o f neglecting a small slow transient
COBELLI et al. (1975), in a later investigation, neglected the slow transient term on the grounds that their experiment lasted only 4h. They argued that 'under this condition, only plasma, liver and bile are involved; tissues are not taken into account because they are only involved at a later time than that accounted for in the experiment.' The model of COBELLI et al. was therefore a 2-compartment one. On the hypothesis that the true plasma clearance curve is given by eqn. 8, the model of COBELLIet al. would be the first two terms of this equation, rescaled to give xl(0) = 100; i.e. xl (t) = 68e- o.o44t + 32e- o.o1ot
(9)
As noted before, u.b.p, varies inversely as the area under the plasma clearance curve. Therefore, if eqn. 8 corresponds to a mean u.b.p, value of 263 rag/day (BERK et al. 1969), then eqn. 9 corresponds to a mean u.b.p, value of 394 mg/day. This demonstrates clearly that the 2-compartment model, derived from data collected over 4h, cannot validly be extrapolated over a period of 48h, the
Medical & Biological Engineering & Computing
March 1979
time needed to account for the effect of the neglected slow transient term. It is interesting to note that COBELLIet aI., using a 2-compartment model, obtained a mean u.b.p, value of 228 mg/day for a group of eight normal subjects. BERK et al. (1969) obtained a mean value of 263 rag/day with a standard deviation of 60 mg/day for 19 studies on normal subjects, so that the mean of the results of CO~ELLt et al. lies within the one standard deviation limit of the results of BERK et al. This was not obtained at the expense of incorrect plasma volume since the mean plasma volume of 2678 ml obtained by COBELLI et al. is within 1% of the mean obtained by BOWDLER (1969) for ten normal subjects. No explanation can be offered for this. 3 Choice of number of compartments The preceding Section has illustrated the need for incorporating the correct number of exponentials in a curve fit to obtain agreement with independent measurements of plasma volume, but, in practice, it frequently proves difficult to do this with the measurements available. To illustrate this, we will lO0
consider data for the plasma clearance of unconjugated 14 C-labelled bilirubin in a cirrhotic patient (cirrhotic patient number 1 in OWENS et al., 1977). The scaled data up to t = 720 min are shown as the circled points in Fig. 3a; the remaining scaled data are tabulated in Table 1. The scaling factor is derived from a least-squares fit of two exponentials to the data x x ( t ) = A e -~t + ( 1 0 0 - A) e -at
(10)
(i.e. with x~(0) = 100). Fitting such a curve to the complete set of data gave xt(t) = 74-47e-~176176176176176
(11)
and this is shown as the solid line curve in Fig. 3a. Sampled values of eqn. 11 for t > 720 rain are tabulated in Table 1. A least-squares fit of three exponentials to the data caused no significant improvement in the Table 1, Data values beyond the time range of Fig. 3a time, min. data
1380
1590
1720
1830
2790
10.1
2'5
3"0
2.5
0"6
eqn. 11
6-4
5"2
4"6
4.1
1 "6
eqn. 13
4"8
3"6
3"0
2"6
0"7
8 o
8
curve fit compared with two exponentials. The 3-exponential family 5O
x l ( t ) = (74.47- C) e-~176
Ce -(~176 +e)t
7D
+25-53e - ~ 1 7 6 1 7 6(12) 176
-6 E
%
;
'
250
5 0
'
750
where 0 < C < 74-47 and e < 0.0161 fits the data equally well. I n practical terms, a 3-compartment model is not identifiable from the data of Fig. 3a.
time,min ~1oo
OWENS et al. (1977) have in fact fitted three exponentials to the data for cirrhotic patient number 1 and have obtained parameter values for the 3compartment model consistent with corresponding parameter values for a further seven cirrhotic patients. Using the data in OWENS et al. (1977), their least-squares fit is given by
O
x
.(2_
E O
u 50
"D
Xl(t) = 27.4e-~176
o
t3
E
+31.5e -~176176(13)
c
00
' 250
, 500
, 750
time, rain b Fig. 3
16.2e -~176176
a Plasma clearance curve of 14C labelled bi/irubin in cirrhotic patient number I of OWENS et aL (1977) o o o scaled data eqn. 11 to the data 9 fitting eqn. 13 to the data b fitting eqn, 16 to the same data
--fitting
M e d i c a l & Biological Engineering & C o m p u t i n g
(with xl(0) = 75.1), the same scaling factor as for eqn. 11 having been used. The dotted line in Fig. 3a is a plot of eqn. 13 up to t = 720 rain; sampled values for t > 720 min are again tabulated in Table 1. Even allowing for the better fit at the tail end of the data (Table 1), the 3-exponential fit is a worse approximation to the data than the 2-exponential fit, the r.m.s, errors being 4.8 and 2.7 normalised units of plasma concentration, respectively. What appears to have happened is that the curvefitting procedure used by OWENS et al. (1977) M a r c h 1979
219
(graphical peeling of exponentials to give initial values of coefficients and exponentials, followed by use of a standard gradient algorithm) has locked into a local minimum. The relatively poor fit explains the lack of any obvious degeneracy (no coefficients of negligible value and no exponents so close together as to be mergeable) in eqn. 13. Applying the Laplace transform to eqns. 11 and 13 gives, respectively
Xl(s) =
100 (s+ 0.0049) (s+ 0.0161)(s+ 0.0010)
(14) 4 Differentiation between models that are inputoutput equivalent
and
XKs)
oscillation. This would be a simple test to implement in future investigations. Naturally, it would be preferable to detect the presence of any oscillation in the concentration of unlabelled bilirubin before the tracer is injected so that the experiment is performed only under conditions where the oscillation is not present. Conversely, any oscillations in existing data should be treated as disturbance effects and no attempt made to model them.
=
75.1 (s+O.O117)(s+O.O056) (s + 0.0145)(s + 0.0093)(s + 0.0014)
(15)
It is interesting to note that, except for the gain constant, eqn. 14 is approximately derivable from eqn. 15 by approximate cancellation of the pole at s = - 0.0093 with zero at s = - 0.0117. Inserting the values of x~(0) obtained from the non-normalised forms of eqn. 11 and 13 in eqn. 3 gives respective plasma volumes of 59.6 ml/kg for the 2-exponential fit and 79.3 ml/kg for the 3exponential fit. The plasma volume for cirrhotic patient number 1 was independently measured by the b.s.p, method to be 47-6 ml/kg (OwENs, 1977) and by the radioactive-chrome red-cell method (STERLINGand GRAY 1950) to be 54-9 ml/kg (OWENS, 1977). The 2-exponential fit gives a lower plasma volume (closer to the independently-measured figures) simply because its zero-time intercept is larger. Inability to detect three exponentials in the data stems from the masking influence of the pronounced oscillation on the data [noticeable oscillation also occurred in the responses of cirrhotic patients numbers 5 and 7 (OWENS et al., 1977)]. A crude approximation to the oscillatory behaviour is given by the expression
Consider a linear n-compartment model containing (p + 1) unknown parameters, namely, p rate constants specifying material transfer rates from compartments to other compartments and to the environment, and an extra parameter specifying the size of the input compartment (the compartment to which the input is applied). The model transfer function contains m coefficients that are algebraic functions of the p rate constants only, and an extra coefficient, the gain constant, which is a function also of the size of input compartment. Physical
path from mouth through g,i. membrane
plasma pool
liver cell pool
Xl(t) = 74.47e -o.ox61t + (25-53-- 4.56 cos 0.0330 e - ~ 1 7 6 1 7 (16) 6176 This is plotted in Fig. 3b and is seen to fit the data reasonably well up to t = 270 min. Eqn. 16 is the impulse response of a 4th-order system, but it cannot be synthesised by a physically-realisable 4-compartment model (i.e. one having positive real rate constants). It is probable that the oscillatory behaviour is not representative of the true tracer kinetics but is indicative of an underlying nonlinear mechanism whereby oscillations in the concentrations of unlabelled unconjugated bilirubin are coupled in the tracer kinetics, Plasma unconjugated bilirubin concentration was in fact estimated during the course of each experiment, but unfortunately not sufficiently frequently to detect unequivocally the presence of 220
path from mouth through g.i.membrane
plasma pool
a31(
liver cell pool b Fig. 4 Alternative models describing metabolic pathways of orally-administored methionine that are input-output equivalent
Medical & Biological Engineering & Computing
March 1979
realisability of the model dictates that the number (rn + 1 - n) of numerator coefficients of the model transfer function satisfies (m + 1 - n) ~< n. The ( r e + l ) coefficients can be identified from inputoutput measurement by standard curve-fitting procedures (BE~EY Ct al., 1974). Clearly, a necessary condition for identifiability of the p rate constants is that p = m. Because the m algebraic relations specified by the p = m rate constants are nonlinear, the parameter estimates are multivalued in general. Ambiguities in parameter values can sometimes be resolved on physiological grounds, but usually require further experimentation involving other input-output ports (DISTEFANO,1976). The above aspects of parameter identifiability have been designated structural identifiability (BELLMAN and /~STR6M, 1970), because algebraic conditions can be formulated in terms of structure only (COBELLI, and ROMANIN-JAcuR, 1976; ZAZWORSKY and KNUDSON, 1977) independent of the (nonzero) parameter values. Two alternative 3-compartment models that have been proposed on physiological grounds to describe the metabolic pathways of orally-administered methionine (BROWN et aI., 1979) are shown ill Figs. 4a and b. In the modelling, an orallyadministered dose of methionine is treated as an impulsive input (Q/ql)5(t) to compartment 1, and sample measurements of plasma concentration are treated as sampled values of the concentration x2(t) in compartment 2. The two models are, in fact, indistinguishable from their plasma responses. For both models, the Laplace transform of the output variable is of the form
K(s+O) X2(s)
= s(s+2,)(s+22)
.
.
.
.
.
(15)
For subject APL (BRowN et aL, 1979), K = 3688 mmoles/l/h -~, 0 = 0.130 h -1, 2a -= 8.74 h -a and 22 = 1-75 h -a, giving respective rate constants for the two models as shown in Table 2.
liver for dynamic measurement is rarely possible; it is however feasible to use, for example, a radioactive tracer and measure its activity within the liver noninvasively. 5 Discussion This paper illustrates, through numerical examples and through the analysis of plasma clearance responses, how measurements currently available can give rise to parameters for compartmental models of metabolic systems that are of questionable accuracy. In terms of model-independent parameters, probably the most important checks of validity of a compartmental model, inaccuracies arise due to having insufficient data points during the decay of the fast exponential, and to the noise level on the data masking the slow exponential, which is usually of small amplitude. These two factors combine to make the system appear to be of lower order than it really is. If fitting of a model of the correct order is attempted, degeneracy occurs, whereby either the coefficient of one or more of the exponentials is very small or two exponents are so close together as to be mergeable. With a noisy response, it is very easy for the gradient algorithm to lock into a local, rather than the global, minimum. These problems could be resolved to some extent by careful control of the environmental conditions of the experiment and ensuring that the unlabelled material concentration is reasonably constant before the start of the experiment. Even then, data points are needed frequently in the early stages of the response to fit the fast exponential adequately, while occasional measurements are needed over a longer time scale than that generally used at present to fit the slow exponential adequately. This points to the obvious need for some continuous (or semicon100
Table 2, Rate constants for the two models
2
M o d e l a21, h -1
~_ ~5
a31, h -1
a12, h -1
823, h -1
a32,h -1
Fig. 4a 2 . 2 6
2'66
5"52
0"06
0
Fig. 4b 2 " 2 6
6"49
0
0'03
1571
The corresponding model predictions of the concentration x3 (t) in compartment 3 are shown in Fig. 5, where x3(t) has been scaled to make the final value 100 units in the case of the model of Fig. 4a. If compartment 3 is interpreted as hepatic, it is seen from Fig. 5 that measurement access to the liver could resolve which of the two models is more representative of the physiological system. Measurement access does not necessarily imply surgical access. In a patient, continuing surgical access to the Medical & Biological Engineering & Computing
mode[ 4b APL
is0 u
-0
I
I
6
12
norma[ised time
Fig. 5 Concentrations in Compartment 3 of the two models shown in Fig. 4, following impulsive input to compartment I
March 1979
221
tinuous, every minute, say) reading instrument to monitor plasma concentration. Even this highly desirable instrument could not resolve the type of ambiguity presented in Section 4, wherein two equal-order models of equal physiological plausibility but giving rise to different liver concentration response curves are proposed. To distinguish the two models in this case, the concentration in the liver pool needs to be monitored. As indicated previously, this may be done noninvasively. Acknowledgment--The work described in this paper is supported by a grant from the UK Science Research Council. References
BEKEY, G. A., UNG, M. T. and KARUZA, S. (1974) Observations on some commonly-used methods for identification of parameters in linear systems. Simulation, Sept. 1974, 69-75. BELLMAN, R. and ASXR6r~, K. J. (1970) On structural identifiability. Math. Jliosciences 7, 329-339. BERK, P. D., HOWE, R. B., BLOOMER,J. R. and BERLIN, N. I. (1969) Studies of bilirubin kinetics in normal adults. J. Clin. Invest. 48, 2176-2190. BERK, P. D., BLOOMER,J. R., HOWE,R. B. and BERLIN, N.I. (1970) Constitutional hepatic dysfunction (Gilbert's Syndrome). Amer. J. Med. 49, 296-305. BROWN, R. F., GODFREY,K. R. and KNELL,A. (1979) Compartmental modelling based on methionine tolerance test data: a case study. Med. & Biol. Eng. 17, 223-229 BOWDLER,A. J. (1969) Regional variations in the proportion of red cells in the blood in man. 2trit. J. Haem. 16, 557-571. CARSON, E. R. and FINKELSTEIN, L. (1977) General problems of validation of control system models in physiology. IFAC symposium on control mechanisms in bio, and ecosystems, Leipzig, Sept. 12-16, 1977.
222
COBELLI, C., FREZZA, M. and TmmELLI, C. (1975) Modelling, identification and parameter estimation of bilirubin kinetics in normal, hemolytic and Gilbert's states. Computers & ~iomed. Res. 8, 522-537. COBELLI, C. and ROMANIN-JACUR,G. (1976) Controllability, observability and structural identifiability of multi-input and multi-output biological compartmental systems. 1EEL Trans. BME-23 93-100. DISTEFANO, J. J. (1976) Tracer experiment design for unique identification of nonlinear physiological systems. Amer. J. Physiol. 230, 476-485. DISTEFANO, J. J. and CAMPFIELD,L. A. (1977) Activity report 1976-77 of the Biocybernetics Laboratory. School of Engineering and Applied Science, UCLA, California. GIBALDI, M. and PERRIER, D. (1975) Pharmacokinetics, Marcel Dekker Publications, New York. JACQUEZ,J. A. (1972) Compartmental analysis in biology and medicine, Elsevier Publishing Company. JONES,E. A., BLOOMER,J. R., BERK,P. D., CARSON,E. R., OWENS, D. and BERLIN,N. I. (1977) The quantitation of hepatic bilirubin synthesis in man. In BERK, P. D. and BERLIN, N. I. (Ed.) The chemistry and physiology of bile pigments, US Department of Health, Education and Welfare, Washington DC, US Government Printing Office, 1977, 189-205. OWENS, D., JONES, E. A. and CARSON, E. R. (1977) Studies on the kinetics of unconjugated x4C bilirubin metabolism in normal subjects and patients with compensated cirrhosis. Clin. Sci. Molec. Med. 52, 555-570. OWENS, D. (1977) Private communication. RANJBARAN,E. and McLAREN, R. W. (1976) Modelling, testing, simulation and optimization of digitalis pharmacokinetics, Proceedings of the IEEE International conference on cybernetics and society, Washington, 1976, 361-367. STERLING, K. and GRAY, S. J. (1950) Determination of the circulating red cell volume in man by radioactive chromium. J. Clin. Invest. 29, 1614-1619. ZAZWORSKY,R. M. and KNUDSEN, H. K. (1977) Comments on COBELLI and ROMANIN-JACUR (1976), IEEE Trans., BME-24, 495-496.
Medical & Biological Engineering & Computing
March 1979
Adequacy of measurements in compartmental modelling of metabolic systems R. F. B r o w n *
E.R. Carson
L. Finkelstein
Department of Systems Science, The City University, London EC1V 0HB, England
K.R. Godfrey
R . P . Jones
Inter-University Institute of Engineering Control, Department of Engineering, University of Warwick, Coventry CV4 7AL, England
A b s t r a c t - - T h e adequacy of present-day measurement techniques for the compartmental modelling of metabolic systems is investigated using numerical examples and analysis of experimentallyobtained p/asma clearance curves. It is concluded that the model parameters obtained are often of questionable accuracy. The situation can be improved by careful monitoring of experimental conditions and judicious spacing of data points on the response curve, but the work shows a clear need for a continuous (or semicontinuous) method of monitoring plasma concentration, To resolve ambiguities between models equal/)/ plausible on physiological grounds, it is necessary to monitor the dynamics of internal variables, for example, of the concentration in the fiver (which is nowadays possible noninvasively). K e y w o r d s - - M e t a b o l i c systems, Modelling 1 Introduction
LINEAR multicompartmental analysis of tracer clearance curves is widely used to yield information about the rates of reaction of chemical substances of interest and their transport into relevant body tissues and fluids (JACQUEZ, 1972; G1BALDI and PERRIER, 1975). Clearly, use of a multicompartmental model to reveal impairment of metabolic pathways in diseased states (JONES et al., 1977) or to formulate drug dosage regimens (RANJBARAN and MCLAREN, 1976) implies acceptance of the validity of the model. This paper concerns itself with the acceptance stage--with the adequacy of model structure and of model parameter values for a given structure (CARSON and FINKELSTE1N, 1977; DISTEFANOand CAMPFIELD, 1977). F o r validation to be of practical worth, it must be with respect to the intended use of the model (CARSON and FINKELSTEZN, 1977). Here, the purpose of the model is assumed to be met if the quantitative representation of the important metabolic pathways afforded by the model is in accord with a priori knowledge, in the following ways. First, the model structure (topology) should be compatible with a priori knowledge of the various internal couplings in the physiological system and their significance. Secondly, the internal states of the model (the compartmental concentrations) should be compatible with a priori knowledge of their dynamic responses (often such a priori knowledge comes from * Currently on sabbaticalleave at the University of Warwick from the University of New South Wales, Kensington 2033, NSW, Australia First received 7th Apri/ and in final form 20th AprU 1978
animal studies and its applicability in the human situation is uncertain). Thirdly, some parameters (such as plasma volume) can be estimated from other experiments without the need to postulate a mathematical model; their values are then tabulated in the literature. Predictions of these parameters from the compartmental model are relatively independent of uncertainties in model structure; e.g. plasma volume can be estimated from the zero time value of a plasma clearance curve, and this is relatively independent of the model postulated. The values of these model-independent parameters should be compatible with the a priori knowledge of their values (see Section 2). Inextricably bound up with the problem of model structure validation is the complementary problem of the reliance that may be placed on the parameter estimates. In a typical measurement of a plasma clearance curve, a dose of tracer material is administered intravenously or orally at time zero, and a succession of plasma samples is taken at times beginning several minutes after the administration and ending several hours (but not usually days) after the administration. The concentration xl(t) of tracer in the plasma is then determined for each sample. The form of the compartmental chosen will yield a series of exponentials to be fitted to the data xl(t) = ~
i=1
Ai exp ( - e ~ t ) . . . . .
(1)
the coefficients then being fitted by standard curve. fitting procedures (BEKEY et al., 1974). Over the time scale of the response, it is not
0140-0118/79/020216 -t-07 $01 -50/0 9 IFMBE: 1979
216
Medical & Biological Engineering & Computing
March 1979
mathematically meaningful to attempt to fit more than three exponentials; if a higher-order model is tried, degeneracy will occur in that either some of the coefficients (As) will be negligible or rate constants (es) will be so close as to be mergeable. Thus, with measurement times from several minutes to several hours, the most complex form of concentration curve that can be fitted is x l ( t ) = A~ exp ( - ~ t t ) + A 2
exp (-cz2 t) + A3 exp ( - e3 t) (2)
parameters of the assumed compartmental structure, are also directly measurable. Such parameters are referred to as model-independent parameters, and they serve a particularly important role in assessing the adequacy of any postulated compartmental model. The best-known model-independent parameter is the plasma volume (or more correctly the initial volume of distribution) that may be inferred from the zero-time intercept of the plasma clearance curve in the following manner. Consider an experiment in which labelled tracer is injected into the plasma, and the subsequent concentration x l ( t ) of tracer in the plasma is measured. Curve fitting of the form of eqn. 1 will give the model-predicted initial concentration xl (0). From a knowledge of the quantity Q o! labelled material injected, the plasma volume q is given in normalised form by
Arbitrarily ordering the terms so that el > c(2 > e3, problems of fitting the first ter m (the fastest exponential) often arise because there are insufficient data points before the term has died away to a small value. Problems of fitting the third term (the slowest exponential) arise because the coefficient A3 is usually small; e.g. BERK et al. (1969) found a value of only 3 % (taking A~ + A2 + A3 = 100) in bilirubin q Q clearance tests on a group of 19 normal subjects. m - mxl(O) . . . . . . . . . (3) (It is not always small, and a value of 12% was found in similar tests on a group of thirteen subjects where m is the weight of the subject. If the tracer is with Gilbert's syndrome (BERN et al., 1970)). radioactive, for example, the usual units are Q in The problems discussed are compounded by the counts per minute, xl(0) in counts per minute per presence of various kinds of disturbances that couple millilitre of plasma, q in millilitres and m in kiIointo the measurements. These include various uncon- grammes. trolled factors that may affect the subject's measured The other main model-independent parameters response behaviour, such as bodily movements, are based on the normalised area under the response aural and visual stimuli (which may affect emotional curve when the input is an impulse-function change and hence physical condition) and whether the in infusion rate (where normalisation is with respect concentration of the material of interest (the to the area under the impulse), or the steady-state unlabelled material) is in a steady state at the outset gain from input to output when the input is a stepof the experiment. function change in infusion rate (these are identical Short-term environmental disturbances often for a linear system). If the experiment is concerned manifest themselves in the measurements by values with the kinetics of bilirubin, for example, the commonly referred to as 'outliers', so cal/ed since model-independent parameter of concern is the they lie outside the region of acceptability. Naturally, unconjugated bilirubin production rate (u.b.p.), every effort needs to be made to minimise environ- which varies inversely as the area under the plasma mental disturbances at their source where this is clearance curve. possible, and by choosing the most favourable time In a typical measurement of a plasma clearance and location for the experiment. Long-duration curve, as noted previously, the first plasma sample is environmental disturbances give rise to trends in the data; these trends are often modelled mathematically 100 by a low-order polynomial, usually the constant to term and linear term only. The constant term (bias) is particularly important when establishing the zero reference as, for example, when subtracting the background count in radioactive tracer experiments. 8 10 9 \.\ Set against the background of the problems "'.'x.__ = data time r a n g e t2 previously described, the remainder of this paper will ID be aimed at determining the adequacy of currentlyTx aD available measurements on which compartmental modelling is based, and at indicating how measurements should be improved so as to arrive at models t._ on which greater reliance can be placed. O C 2 Model-independent parameters Certain parameters of tracer response curves, as well as being inferrable from knowledge of the M e d i c a l & Biological Engineering & C o m p u t i n g
0
06
i
i
J
10
20
30
time F,'g. I Logarithmic plot of eqns. 4 and 5
M a r c h 1979
217
taken several minutes after injection of the radioactive tracer and the last sample several hours (but not usually days) after injection. A compartmental model fitted to this data is then used to extrapolate the data to zero time to deduce the steady-state gain of the system. Implicit in this procedure is the assumption that the model structure is valid outside the time range of data used in estimating the model parameters. An important function of experiment design is to perform studies to test the validity of such assumptions. Some numerical examples will now be presented to illustrate likely problems.
is incorrectly estimated to be 9.09 units (true value 18.18 units). E x a m p l e 2: effects o f neglecting a small time delay
Consider again the idealised plasma clearance curve of the preceding example, but with a small time delay inserted into the response so that the curve is defined by x~(t) = O, O ~< t < 0 . 1 = 9.09e-r176
-a~176
t >! 0.1
(6)
If measurement does not begin until after t = 0.1, then eqn. 6 may validly be written
E x a m p l e 1 : effects o f neglecting a fast transient
Fig. 1 shows a logarithmic plot of an analytic idealisation of a plasma clearance curve, defined in normalised units by
x~(t) = l O . 0 5 e - t + 2 4 7 . 1 2 e - l ~
t > 0.1
(7)
If this equation is fitted to the data
(4)
(a) the plasma volume is incorrectly estimated to be 0.39 units (true value 1 unit)
[with xl(0) = 100]. A possible compartmental model which would represent eqn. 4 is shown in Fig. 2a, where the volume of compartment 1 (the plasma compartment) is defined to be ql = 1 unit. Assume that measurements of plasma concentration are made over the time range 0.6 to 3 units and that the curve
(b) the area under the curve is incorrectly estimated to be 34.76 units (true value 18.18 units).
x~(t) = 9 . 0 9 e - t + 9 0 . 9 1 e -a~
x~(t) = 9.09e -t .
.
.
.
. . . .
.
.
.
.
(5)
is fitted to the data. The compartmental model corresponding to eqn. 5 is shown in Fig. 2b.
"(~1.09 0 '73
q1=1
a
9 - 0 9 6 ( t ) ~
1"00
ql =11 b Fig. 2 a Compartmental eqn. 4
model
corresponding
to
b Compartmental eqn. 5
model
corresponding
to
In this example, neglecting the fast transient term 90.91e - l ~ has dramatic effects: (a) the zero-time intercept is 9.09 (eqn. 5) instead of 100 (eqn. 4), so that the plasma volume is, from eqn. 3, incorrectly estimated to be 11 units (true value 1 unit) (b) the area under the curve, which from eqn. 1 is given by ~A~ i=1
218
Q{i
BERK et al. (1969), in a study of bilirubin metabolism in 13 normal subjects, obtained the following mean plasma clearance curve for unconjugated bilirubin: xz(t) = 66e-~176 + 31e-~176176 + 3 e - ~ 1 7 6 1 7 6 (8)
(normalised so that xl(0) = 100). The small-valued slow transient term 3e -~176176 was essential in their modelling in order to obtain a correct area under the curve, and hence to obtain correct values for unconjugated bilirubin production rate and redblood-cell lifespan.
9.18
1006(t) --C)=
Example 3: effeets o f neglecting a small slow transient
COBELLI et al. (1975), in a later investigation, neglected the slow transient term on the grounds that their experiment lasted only 4h. They argued that 'under this condition, only plasma, liver and bile are involved; tissues are not taken into account because they are only involved at a later time than that accounted for in the experiment.' The model of COBELLI et al. was therefore a 2-compartment one. On the hypothesis that the true plasma clearance curve is given by eqn. 8, the model of COBELLIet al. would be the first two terms of this equation, rescaled to give xl(0) = 100; i.e. xl (t) = 68e- o.o44t + 32e- o.o1ot
(9)
As noted before, u.b.p, varies inversely as the area under the plasma clearance curve. Therefore, if eqn. 8 corresponds to a mean u.b.p, value of 263 rag/day (BERK et al. 1969), then eqn. 9 corresponds to a mean u.b.p, value of 394 mg/day. This demonstrates clearly that the 2-compartment model, derived from data collected over 4h, cannot validly be extrapolated over a period of 48h, the
Medical & Biological Engineering & Computing
March 1979
time needed to account for the effect of the neglected slow transient term. It is interesting to note that COBELLIet aI., using a 2-compartment model, obtained a mean u.b.p, value of 228 mg/day for a group of eight normal subjects. BERK et al. (1969) obtained a mean value of 263 rag/day with a standard deviation of 60 mg/day for 19 studies on normal subjects, so that the mean of the results of CO~ELLt et al. lies within the one standard deviation limit of the results of BERK et al. This was not obtained at the expense of incorrect plasma volume since the mean plasma volume of 2678 ml obtained by COBELLI et al. is within 1% of the mean obtained by BOWDLER (1969) for ten normal subjects. No explanation can be offered for this. 3 Choice of number of compartments The preceding Section has illustrated the need for incorporating the correct number of exponentials in a curve fit to obtain agreement with independent measurements of plasma volume, but, in practice, it frequently proves difficult to do this with the measurements available. To illustrate this, we will lO0
consider data for the plasma clearance of unconjugated 14 C-labelled bilirubin in a cirrhotic patient (cirrhotic patient number 1 in OWENS et al., 1977). The scaled data up to t = 720 min are shown as the circled points in Fig. 3a; the remaining scaled data are tabulated in Table 1. The scaling factor is derived from a least-squares fit of two exponentials to the data x x ( t ) = A e -~t + ( 1 0 0 - A) e -at
(10)
(i.e. with x~(0) = 100). Fitting such a curve to the complete set of data gave xt(t) = 74-47e-~176176176176176
(11)
and this is shown as the solid line curve in Fig. 3a. Sampled values of eqn. 11 for t > 720 rain are tabulated in Table 1. A least-squares fit of three exponentials to the data caused no significant improvement in the Table 1, Data values beyond the time range of Fig. 3a time, min. data
1380
1590
1720
1830
2790
10.1
2'5
3"0
2.5
0"6
eqn. 11
6-4
5"2
4"6
4.1
1 "6
eqn. 13
4"8
3"6
3"0
2"6
0"7
8 o
8
curve fit compared with two exponentials. The 3-exponential family 5O
x l ( t ) = (74.47- C) e-~176
Ce -(~176 +e)t
7D
+25-53e - ~ 1 7 6 1 7 6(12) 176
-6 E
%
;
'
250
5 0
'
750
where 0 < C < 74-47 and e < 0.0161 fits the data equally well. I n practical terms, a 3-compartment model is not identifiable from the data of Fig. 3a.
time,min ~1oo
OWENS et al. (1977) have in fact fitted three exponentials to the data for cirrhotic patient number 1 and have obtained parameter values for the 3compartment model consistent with corresponding parameter values for a further seven cirrhotic patients. Using the data in OWENS et al. (1977), their least-squares fit is given by
O
x
.(2_
E O
u 50
"D
Xl(t) = 27.4e-~176
o
t3
E
+31.5e -~176176(13)
c
00
' 250
, 500
, 750
time, rain b Fig. 3
16.2e -~176176
a Plasma clearance curve of 14C labelled bi/irubin in cirrhotic patient number I of OWENS et aL (1977) o o o scaled data eqn. 11 to the data 9 fitting eqn. 13 to the data b fitting eqn, 16 to the same data
--fitting
M e d i c a l & Biological Engineering & C o m p u t i n g
(with xl(0) = 75.1), the same scaling factor as for eqn. 11 having been used. The dotted line in Fig. 3a is a plot of eqn. 13 up to t = 720 rain; sampled values for t > 720 min are again tabulated in Table 1. Even allowing for the better fit at the tail end of the data (Table 1), the 3-exponential fit is a worse approximation to the data than the 2-exponential fit, the r.m.s, errors being 4.8 and 2.7 normalised units of plasma concentration, respectively. What appears to have happened is that the curvefitting procedure used by OWENS et al. (1977) M a r c h 1979
219
(graphical peeling of exponentials to give initial values of coefficients and exponentials, followed by use of a standard gradient algorithm) has locked into a local minimum. The relatively poor fit explains the lack of any obvious degeneracy (no coefficients of negligible value and no exponents so close together as to be mergeable) in eqn. 13. Applying the Laplace transform to eqns. 11 and 13 gives, respectively
Xl(s) =
100 (s+ 0.0049) (s+ 0.0161)(s+ 0.0010)
(14) 4 Differentiation between models that are inputoutput equivalent
and
XKs)
oscillation. This would be a simple test to implement in future investigations. Naturally, it would be preferable to detect the presence of any oscillation in the concentration of unlabelled bilirubin before the tracer is injected so that the experiment is performed only under conditions where the oscillation is not present. Conversely, any oscillations in existing data should be treated as disturbance effects and no attempt made to model them.
=
75.1 (s+O.O117)(s+O.O056) (s + 0.0145)(s + 0.0093)(s + 0.0014)
(15)
It is interesting to note that, except for the gain constant, eqn. 14 is approximately derivable from eqn. 15 by approximate cancellation of the pole at s = - 0.0093 with zero at s = - 0.0117. Inserting the values of x~(0) obtained from the non-normalised forms of eqn. 11 and 13 in eqn. 3 gives respective plasma volumes of 59.6 ml/kg for the 2-exponential fit and 79.3 ml/kg for the 3exponential fit. The plasma volume for cirrhotic patient number 1 was independently measured by the b.s.p, method to be 47-6 ml/kg (OwENs, 1977) and by the radioactive-chrome red-cell method (STERLINGand GRAY 1950) to be 54-9 ml/kg (OWENS, 1977). The 2-exponential fit gives a lower plasma volume (closer to the independently-measured figures) simply because its zero-time intercept is larger. Inability to detect three exponentials in the data stems from the masking influence of the pronounced oscillation on the data [noticeable oscillation also occurred in the responses of cirrhotic patients numbers 5 and 7 (OWENS et al., 1977)]. A crude approximation to the oscillatory behaviour is given by the expression
Consider a linear n-compartment model containing (p + 1) unknown parameters, namely, p rate constants specifying material transfer rates from compartments to other compartments and to the environment, and an extra parameter specifying the size of the input compartment (the compartment to which the input is applied). The model transfer function contains m coefficients that are algebraic functions of the p rate constants only, and an extra coefficient, the gain constant, which is a function also of the size of input compartment. Physical
path from mouth through g,i. membrane
plasma pool
liver cell pool
Xl(t) = 74.47e -o.ox61t + (25-53-- 4.56 cos 0.0330 e - ~ 1 7 6 1 7 (16) 6176 This is plotted in Fig. 3b and is seen to fit the data reasonably well up to t = 270 min. Eqn. 16 is the impulse response of a 4th-order system, but it cannot be synthesised by a physically-realisable 4-compartment model (i.e. one having positive real rate constants). It is probable that the oscillatory behaviour is not representative of the true tracer kinetics but is indicative of an underlying nonlinear mechanism whereby oscillations in the concentrations of unlabelled unconjugated bilirubin are coupled in the tracer kinetics, Plasma unconjugated bilirubin concentration was in fact estimated during the course of each experiment, but unfortunately not sufficiently frequently to detect unequivocally the presence of 220
path from mouth through g.i.membrane
plasma pool
a31(
liver cell pool b Fig. 4 Alternative models describing metabolic pathways of orally-administored methionine that are input-output equivalent
Medical & Biological Engineering & Computing
March 1979
realisability of the model dictates that the number (rn + 1 - n) of numerator coefficients of the model transfer function satisfies (m + 1 - n) ~< n. The ( r e + l ) coefficients can be identified from inputoutput measurement by standard curve-fitting procedures (BE~EY Ct al., 1974). Clearly, a necessary condition for identifiability of the p rate constants is that p = m. Because the m algebraic relations specified by the p = m rate constants are nonlinear, the parameter estimates are multivalued in general. Ambiguities in parameter values can sometimes be resolved on physiological grounds, but usually require further experimentation involving other input-output ports (DISTEFANO,1976). The above aspects of parameter identifiability have been designated structural identifiability (BELLMAN and /~STR6M, 1970), because algebraic conditions can be formulated in terms of structure only (COBELLI, and ROMANIN-JAcuR, 1976; ZAZWORSKY and KNUDSON, 1977) independent of the (nonzero) parameter values. Two alternative 3-compartment models that have been proposed on physiological grounds to describe the metabolic pathways of orally-administered methionine (BROWN et aI., 1979) are shown ill Figs. 4a and b. In the modelling, an orallyadministered dose of methionine is treated as an impulsive input (Q/ql)5(t) to compartment 1, and sample measurements of plasma concentration are treated as sampled values of the concentration x2(t) in compartment 2. The two models are, in fact, indistinguishable from their plasma responses. For both models, the Laplace transform of the output variable is of the form
K(s+O) X2(s)
= s(s+2,)(s+22)
.
.
.
.
.
(15)
For subject APL (BRowN et aL, 1979), K = 3688 mmoles/l/h -~, 0 = 0.130 h -1, 2a -= 8.74 h -a and 22 = 1-75 h -a, giving respective rate constants for the two models as shown in Table 2.
liver for dynamic measurement is rarely possible; it is however feasible to use, for example, a radioactive tracer and measure its activity within the liver noninvasively. 5 Discussion This paper illustrates, through numerical examples and through the analysis of plasma clearance responses, how measurements currently available can give rise to parameters for compartmental models of metabolic systems that are of questionable accuracy. In terms of model-independent parameters, probably the most important checks of validity of a compartmental model, inaccuracies arise due to having insufficient data points during the decay of the fast exponential, and to the noise level on the data masking the slow exponential, which is usually of small amplitude. These two factors combine to make the system appear to be of lower order than it really is. If fitting of a model of the correct order is attempted, degeneracy occurs, whereby either the coefficient of one or more of the exponentials is very small or two exponents are so close together as to be mergeable. With a noisy response, it is very easy for the gradient algorithm to lock into a local, rather than the global, minimum. These problems could be resolved to some extent by careful control of the environmental conditions of the experiment and ensuring that the unlabelled material concentration is reasonably constant before the start of the experiment. Even then, data points are needed frequently in the early stages of the response to fit the fast exponential adequately, while occasional measurements are needed over a longer time scale than that generally used at present to fit the slow exponential adequately. This points to the obvious need for some continuous (or semicon100
Table 2, Rate constants for the two models
2
M o d e l a21, h -1
~_ ~5
a31, h -1
a12, h -1
823, h -1
a32,h -1
Fig. 4a 2 . 2 6
2'66
5"52
0"06
0
Fig. 4b 2 " 2 6
6"49
0
0'03
1571
The corresponding model predictions of the concentration x3 (t) in compartment 3 are shown in Fig. 5, where x3(t) has been scaled to make the final value 100 units in the case of the model of Fig. 4a. If compartment 3 is interpreted as hepatic, it is seen from Fig. 5 that measurement access to the liver could resolve which of the two models is more representative of the physiological system. Measurement access does not necessarily imply surgical access. In a patient, continuing surgical access to the Medical & Biological Engineering & Computing
mode[ 4b APL
is0 u
-0
I
I
6
12
norma[ised time
Fig. 5 Concentrations in Compartment 3 of the two models shown in Fig. 4, following impulsive input to compartment I
March 1979
221
tinuous, every minute, say) reading instrument to monitor plasma concentration. Even this highly desirable instrument could not resolve the type of ambiguity presented in Section 4, wherein two equal-order models of equal physiological plausibility but giving rise to different liver concentration response curves are proposed. To distinguish the two models in this case, the concentration in the liver pool needs to be monitored. As indicated previously, this may be done noninvasively. Acknowledgment--The work described in this paper is supported by a grant from the UK Science Research Council. References
BEKEY, G. A., UNG, M. T. and KARUZA, S. (1974) Observations on some commonly-used methods for identification of parameters in linear systems. Simulation, Sept. 1974, 69-75. BELLMAN, R. and ASXR6r~, K. J. (1970) On structural identifiability. Math. Jliosciences 7, 329-339. BERK, P. D., HOWE, R. B., BLOOMER,J. R. and BERLIN, N. I. (1969) Studies of bilirubin kinetics in normal adults. J. Clin. Invest. 48, 2176-2190. BERK, P. D., BLOOMER,J. R., HOWE,R. B. and BERLIN, N.I. (1970) Constitutional hepatic dysfunction (Gilbert's Syndrome). Amer. J. Med. 49, 296-305. BROWN, R. F., GODFREY,K. R. and KNELL,A. (1979) Compartmental modelling based on methionine tolerance test data: a case study. Med. & Biol. Eng. 17, 223-229 BOWDLER,A. J. (1969) Regional variations in the proportion of red cells in the blood in man. 2trit. J. Haem. 16, 557-571. CARSON, E. R. and FINKELSTEIN, L. (1977) General problems of validation of control system models in physiology. IFAC symposium on control mechanisms in bio, and ecosystems, Leipzig, Sept. 12-16, 1977.
222
COBELLI, C., FREZZA, M. and TmmELLI, C. (1975) Modelling, identification and parameter estimation of bilirubin kinetics in normal, hemolytic and Gilbert's states. Computers & ~iomed. Res. 8, 522-537. COBELLI, C. and ROMANIN-JACUR,G. (1976) Controllability, observability and structural identifiability of multi-input and multi-output biological compartmental systems. 1EEL Trans. BME-23 93-100. DISTEFANO, J. J. (1976) Tracer experiment design for unique identification of nonlinear physiological systems. Amer. J. Physiol. 230, 476-485. DISTEFANO, J. J. and CAMPFIELD,L. A. (1977) Activity report 1976-77 of the Biocybernetics Laboratory. School of Engineering and Applied Science, UCLA, California. GIBALDI, M. and PERRIER, D. (1975) Pharmacokinetics, Marcel Dekker Publications, New York. JACQUEZ,J. A. (1972) Compartmental analysis in biology and medicine, Elsevier Publishing Company. JONES,E. A., BLOOMER,J. R., BERK,P. D., CARSON,E. R., OWENS, D. and BERLIN,N. I. (1977) The quantitation of hepatic bilirubin synthesis in man. In BERK, P. D. and BERLIN, N. I. (Ed.) The chemistry and physiology of bile pigments, US Department of Health, Education and Welfare, Washington DC, US Government Printing Office, 1977, 189-205. OWENS, D., JONES, E. A. and CARSON, E. R. (1977) Studies on the kinetics of unconjugated x4C bilirubin metabolism in normal subjects and patients with compensated cirrhosis. Clin. Sci. Molec. Med. 52, 555-570. OWENS, D. (1977) Private communication. RANJBARAN,E. and McLAREN, R. W. (1976) Modelling, testing, simulation and optimization of digitalis pharmacokinetics, Proceedings of the IEEE International conference on cybernetics and society, Washington, 1976, 361-367. STERLING, K. and GRAY, S. J. (1950) Determination of the circulating red cell volume in man by radioactive chromium. J. Clin. Invest. 29, 1614-1619. ZAZWORSKY,R. M. and KNUDSEN, H. K. (1977) Comments on COBELLI and ROMANIN-JACUR (1976), IEEE Trans., BME-24, 495-496.
Medical & Biological Engineering & Computing
March 1979
