A STRATEGY OF APPROACHING THE FIT OF A COMPARTMENTAL MODEL APPLIED TO THE BEHAVIOUR OF ASPIRINS IN HUMANS

C. REWLLARD,M. GNESSEN,C. V. Pmmm and J. R. SCHEEUUR

et DWion dr ?hammc&gie Clintqw. UniversityHospital,Genew (Switzerland)

Uaitk d’b$ormatiqne

(Received: 15 January, 1975)

SUMMARY A compartmental model has been set up with five parts. Formulation and testrng of the model have been performed by simultaneous experiments involving measurement of blood levels of ASA and SA for four diflerent preparations of aspirin investigated in four dtrerent patients. The hypotheses of the present study had been based upon a four-compartment model. However, such a model did not appear to be adequate and a new model with jive compartments has been substituted for the former one. The five-compartment model has led to the formulation and testing of new physiological hypotheses: it has also provided the development of a fitting method using powerful convergence algorithms. A general analytical formulation has been set up which has made it possible to circumscribe the multiple roots of the characteristic polynomial when they have to be taken into consideration.

SOMMAIRE Un motile compartimental. h cinq compartiments a Cte formule et test& sur ies observations synchrones d’ASA et de SA du sang pour quatre types d’aspirine et chez quatre sujets diflkrents. Ce nouveau modsie, en lieu et place de celui, inapproprit d quatre compartiments, qui correspondaient aux hypotheses de depart, a et& non seuiement l’occasion de formulation et test d’hypotht%es physiologiques nouvelies. mais aussi ie point de depart du d&loppement dune mtthode d’ajustement utilisant &s algorithmes de convergences puissants h partir d’une formulation analytique gt!nPrale et complete notamment en ce qui concerne la detection puis l’utilisation de racines multiples du polyn6me caracteristique lorsqu’il y a lieu de les considerer. 131 Int. J. Bio-MedicalComputing(6) (1975)-O Applied Science Publishers Ltd, Enghd, 1975 Priitcd in Great Britain

132

C. REVILLARD,

M. GRIESSEN,

C. V. PERRIER,

J. R. SCHERRER

INTRODUCTION

A physical model is used in order to study the behaviour of aspirin in humans. Experimental data are collected by measuring the levels of aspirin and salicylate in the blood. A fitting method had to be found for compartmental models completely formulated on the basis of data referring only to a limited number of compartments with respect to the model considered (Rowland et al., 1972; Brownell et al., 1967; Berman and Schoenfeld, 1956; Lewallen et al., 1959; Robertson, 1957; Rescigno and Segre, 1966; Jacquez, 1972). The method commonly used is the SAAM method (Berman et al., 1962; Berman and Weiss, 1966). Another method has been developed in the present work in order to facilitate the investigation of aspirin in humans; it derives from the fitting methods mainly used in experimental physics (James and ROOS, 1971). The estimated values of the exponential parameters of the model are very close. Consequently the number of exponential terms required for the fitting of the observations was smaller than the number of compartments of the proposed physical model. Therefore it was not advisable to keep the number of the model’s compartments to the number of exponential terms. Similarly, it was not possible to use the parameters of the sum of exponential terms fitted in order to correctly determine the possible fields of recruitment of initial estimates satisfying parameters of the matrix of the differential equations describing the model (Lewallen et al., 1959). In addition, when the hypothesis of indistinct exponents (that is to say the case of the root of the characteristic polynomial having a multiplicity greater than one) has to be kept, the formal relationship permitting the step from exponential solutions to differential equations (Berman and Schoenfeld, 1956; Hearon, 1969) becomes analytically too complicated to investigate the possible fields of good initial parameter estimation. In other respects the use of the SAAM method needs good initial estimates of the parameters: convergence of the fitting is realised only when the initial parameters are close to their final value. Indeed, the algorithm used takes into account only the first term of Taylor’s series development of the functions representing the behaviour of the substance in a given compartment, at every moment (Berman and Weiss, 1966). Assuming that the initial estimates are random ones, a new chain of programs can be developed from a set of programmed algorithms of convergence that vary in performance (James and Roos, 1971). In spite of bad initial values of the parameters of the differential equations, such a set of algorithms has been able to vary these parameters in order to obtain a satisfying fit of the collected data. When such a set of programs was in use, the test of a four-compartment model, as suggested by physiological data (Rowland et al., 1972), was not decisive, since it was impossible to obtain a satisfying fit of the observations. Fitting of the results became systematically satisfying only when five compartments were involved

COMPARTMENTAL

MODEL

133

in the model. The goodness of the fit of the data was comparable not only between the tested patients but also between the different preparations of aspirin. Material

A group of normal patients was investigated and four commercially-available preparations of acetylsalicylic acid tested. Blood samples were collected at regular intervals in each patient in order to measure the blood level of acetylosalicylic acid and salicylate. Chemical analysis of plasma was performed in duplicate on blood samples collected according to the following schedule: interval 0 to 30 min: blood sampling every 5 min; interval 30 to 60 min: blood sampling every 10 min; later samples were collected at 90, 120, 180, 300, 360 and 570 mm. The clinical and chemical methods used in the present work have been described in detail by Griessen et al. (in press). Acetylsalicylic (ASA) and salicylate (SA) levels have been measured in moles (M). Sampling errors during the experimental procedure have been reestimated on the following basis: the calculated value of chi-square should approximately correspond to the numbers of degrees of freedom. According to the clinical and chemical methods used in the present work, the results are considered to be stochastically independant. MATHEMATICAL

FORMULATION

OF A COMPARTMENTAL

MODEL

Nomenclature

The nomenclature used in the present work is in agreement with the conventions set up by Brownell et al (1967).

I: Pij

k, = K,

k Oi t n cij

Ai

V

=

p,ijlsi

1

kji

quantity of substance considered in compartment i at time t. variation of the quantity of substance in compartment i, as a function of time. quantity of substance crossing from compartment j to compartment i per unit time. rate constant for transfer of compartment j to compartment i. sum of the rate constants leaving compartment i.

j-0 j#'I

rate constant bound to the substance leaving the system from compartment i. time. number of compartments of the model considered. coefficient of the exponential terms of the solution of the differential equation system. exponent of the same exponential terms. volume in which the substance to be measured might be found.

C. REVILLARD,

134

hf. GRIIBSEN, C. V. PERRIER, J. R. SCHERRJ3t

The di$erential equations Using the nomenclature given above, a compartmental model can be completely described by a system of differential equations which are linear, homogeneous, of first order and with constant coefficients: # 8, =

c $5:

i=

kvSj - KiSi,

(1)

1,n

Such a system can be written as a matrix form:

(2)

3 = KS where S and 2! are vectors and K is an n by n matrix semi-definite negative: -K&i&,3 - - +kin k 21 - K2kz3.. . k,, K=

:

(3)

k.ik,,2knJ . .,. -K,, Using a matrix nomenclature, the solution of eqn. (1) is as follows: S = exp (Kt) S,,

(4)

where S, is the vector of the initial values. The practical resolution of this problem has been widely described. The method chosen in the present study has been used by Gantmacher (1959a, b). It consists in expressing analytically the matrix exp (Kt), and is summarised in the Appendix. Knowing the matrix K, defined above, and noticing that the characteristic values of Kare the 1, defined above, according to the formalism described in the Appendix, qn. (4) can be written: (5)

S = I(K)& and being derived from Davidon (1968), withf(K)

= exp (Kt), it can also be written:

,

S = 2 [exp (&t)&. i + t exp (&t)Z,, z + - - - + tme- ’ exp (A,t)Z,,,,JS, c-l By multiplying vector SO of the initial values by matrices &, obtained and eqn. (6) is now written: S =

2

exp (A.t)f$ I + t exp (&t)C,,,

(6)

vectors C, can be

+ - - - + t’“c- 1 exp (&t)C #.lnC

(7)

e-l This equation can be written in matrix form: s=cv

03)

COWARTMENTAL

MODEL

135

Matrix C with the size of n by m (m I n) has for columns the m vectors:

which are defined above. Vector V with the size m, exp being the exponential base, is composed of: exp (Lit), texp (A,t), . . ., fml-l P-

exp (A,t), exp (A,f), t exp (A,?), . . ., ’ exp (A,t), . . ., exp (A,t), . . ., fmrT1 exp (St).

If i is derived from eqn. (8). it can be written: LCAV

(9)

Matrix A, with the size m by m, is composed of r blocks with the size m,, each corresponding to a particular root of the minimal polynomial of matrix K. If the root Ai has an order of multiplicity mi, the corresponding block will be: -

mi .

Ai

0

. . .

0

I 0

Ai 2

0 Izi

0 0

0

.

.

... mi

.

‘-

I -W?i

.

I ;ti i

If the values of 3 and of S obtained from eqns. (8) and (9) are now replaced in eqn. (2), it follows: CAY

= KCV

(10)

This relationship holds true for any value off; thus it follows: CA = KC

(11)

Finally, if the characteristic and minimal polynomials of K have the same degree n, itfollows: K = CAC-(12) In the particular case where the Ai are all different, eqn. (7) can be simplified as follows : n

S =

2

exp (A$) Ci

(13)

i=l

In this case only the matrix C is the matrix of the characteristic vectors of K; and the matrix A is reduced to a diagonal matrix. Its components are the characteristic values Ai of the matrix K. Equation (l3), which expresses a particular case, cannot be applied to any supposed compartmental system.

136

C. REVILLARD, M. GRIESSEN, C. V. PERRIER, J. R. SCHERRER

Thus it is possible to calculate the elements kij of the matrix K by fitting the measures initially by using functions of the type of eqn. (7) rather than functions of the type of eqn. (13). Such a procedure is advised when available measures are referring to each of the compartments of the supposed model. Such fits allow the matrices C and A to be set up. Equation (12) is particularly appropriate for calculating the elements of matrix K. Furthermore, it is not possible to define the order of matrix K or the number, n, of compartments corresponding to the data when starting from fits according to eqn. (17) since data are not available for every compartment. Indeed such a procedure consists in defining the degree, m, of the minimal polynomial of K. Knowing that m I n, it can only be said that the number of compartments used is a minimum and that it may not correspond to the true number. A model could be set up on the basis of a previously existing model. The data could then be directly fitted from the k, determining the modelunderconsideration. It is quite probable that the number of compartments enclosed in this model shows up to be too small. In such a case a function of the type of eqn. (7) may be useful to test the correctness of new hypotheses. However that may be, the elements kij of the supposed model have to be fitted to the observations. Such a fit can be directly drawn from the k,,; the initial values of the kij can be very distant from their final values provided that the iterative method of fitting is assumed by a powerful convergence algorithm. On the other hand, when the convergence algorithm is less efficient, the best method for the initial kij is that described by Berman and Schoenfeld (1956). This method consists in recruiting initial values as close as possible to their final values by performing linear programming starting from the known constraints of the matrices K, C, and A. However, the method of Berman and Schoenfeld is not always easy to handle. Furthermore, Berman himself points out that his method needs great care. Indeed, within a limited time interval many differential functions may fit data satisfactorily (Julius, 1972). However, such a method has to be considered when the algorithm used to fit the kij to the experimental data does not converge to a satisfying solution, the kij being too far from their final values. Fortunately such a step could be avoided in the present study: a series of powerful convergence algorithms have provided good results even when starting from random kij.

DESCRIPTION OF THE CONVERGENCE ALGORITHMS

The MINUIT program (James and ROOS,1971) has been used to fit the components kij of the present model. This program has been developed by the CERN (Centre Europeen de Recherches Nucleaires). It permits the minimisation of a function of

COMPARTMENTAL

MODEL

137

n variables by one of four methods. These four methods can be briefly described in the following way: Calculation of a minimum by a Monte-Carlo method (SEEK) (James and ROOS, 1971). This method can be used when no reasonable starting point is available. Minimisation on the basis of the method developed by Rosenbrock and Sheppey (TAUROS) (Rosenbrock, 1960; Sheppey, 1968). This algorithm is an extension of the method based upon the variation of coordinates. Algorithm based upon the simplex method (SIMPLX) (Nedler and Mead, 1967). Method providing a simultaneous convergence towards a minimum and the true matrix of covariance when close to the minimum (MIGRAD) (Davidon, 1968). Moreover, the MINUIT program (James and Roos, 1971) allows one to calculate the true errors by using a method developed by its authors. For each available sample an initial minimum has been calculated using the algorithm SEEK. Minimisation has been further performed through TAUROS or SIMPLX. Finally, MIGRAD has provided the convergence and estimation of the errors of the parameters kij. Next the kij found have been checked by again using SEEK with the final parameters to make sure that they correspond to a true minimum of the function. Starting from parameters with initial values generally more than ten times greater or less than the final values, the minimum of the function has been obtained (11 variable parameters) after 2000 iterations. Since the function is calculated at least once at each iteration, it is advisable to use an integration method as rapid as possible in order to reduce the time used for calculation.

METHOD OF INTEGRATION

An algorithm of pseudo-analytical integration has been developed in the present work. It is based upon the calculation of a matrix function (Gantmacher, 1959a, 6); and has been compared with several methods of numerical integration (Martin, 1958; Hanning, 1959; Bulirsch and Stoer, 1966). In all cases the algorithm used in the present study has been 4 to 10 times more rapid than the numerical algorithms quoted above. In the present analysis the elements of vector S, of the initial values are all equal to zero, except the first which is equal to 1. Such a condition appreciably reduces the number of operations needed for the integration and may partly account for the rapidity of the procedure. In any case the pseudo-analytical integration gives at any time good information about the structure of the matrix K and the plausibility of the magnitude of kij,

138

C. REVILLARD, M. CRIESSEN, C. V. PERRIER, J. R. SCHERRER

(7 Compute Xl

no

1

(iy!q

*

cr> T

Compute matrlces C

~2

calculation

0

Fig. I.

Flow chart of the calculation program of the kQ through iterations. The box ‘kv variation by MINUIT’ represents the various convergence algorithms used.

since the elements

of the matrices

not the case when a numerical The algorithm (a) Calculation 1961/62). (b) Calculation

of integration

of multiple

at every iteration.

iLi, of the matrix

the following

steps (see Fig. I):

K by the method

QR (Francis,

roots. Two roots are equal when the ratio

than 0.02. When all the roots are distinct

of the coefficients (c) Otherwise,

is directly calculation

from the distinct

roots the matricial

Such is

is used.

used here performs

of the eigenvalues

Max (/I.,[, liLjl) is smaller

C and A are defined

method

calculated. of the minimal product

polynomial

of matrix

is first calculated:

j/ii -

the matrix

~jl/ C

K. Starting

COMPARTMENTAL

P = n

MODEL

139

(K - liE)

(14)

i=l

where E means the matrix unit. Each term of the matrix P has to be, in terms of absolute value, smaller than lo- ‘. If such is the case, the minimal polynomial can be expressed. Otherwise a matrix P has to be calculated with regard to the multiplicity of the li until the test of the elements of P comes to a satisfying value. In any case, when the number of the terms of the product reaches the level of the matrix K the test is satisfactory. (d) Calculation of the elements of the matrices Zij starting from the minimal polynomial, according to a method described in Gantmacher (1959~). As regards the structure of the vector S, of the initial values, only the elements of the first column of each matrix Zij are useful for the next calculations. The products ZijS, are thus implicitly calculated. They supply the vectors Cij which are used to set up the matrix C of the coefficients. However, in case the roots should all be distinct the theorem of Sylvester is used as follows: n

f(K) =

2 f(li)xi i=l

with : n

r-r

(K - ljE)

Xi

=

j#i ‘3’

r-I i+i j=l

(15) CAi

-

nj)

Such a procedure cuts down one-third of the calculation time needed for the integration. Once the elements of the matrix C have been calculated, the elements of the vector S can be calculated from eqn. (8); comparing these elements with the measurements, chi-square can be formulated and then minimised.

ANALYTICAL

FORMULATION

AND

TESTS OF THE MODEL

Physiological hypotheses

The hypotheses concerning the physiological behaviour of aspirin can be completely described by Fig. 2 which is a deterministic compartmental model. First, a model consisting of four compartments has been set up (Fig. 2(A)). Second, on the basis of the results presented further in the example of Table I, a five compartment

140

C. REVILLARD,

M. GRIESSEN,

C. V. PERRIER,

J. R. SCHERRER

B

Figs. 2(A) and 2(B). Model 1A contained all the starting hypotheses of the present study. It did not provide a simultaneous satisfying fit of ASA and SA. On the other hand, model lB, by admitting an additional compartment, has evidenced a satisfying fit of all the observations concerning the four different preparations and the four patients studied. In Figs. 2(A) and 2(B) compartment 1 represents the stomach where the substance ASA has been directly introduced. It has been assumed that the total dose of 1 g of ASA has reached the stomach at time zero. Compartments 2 and 3 constitute plasma where the levels of ASA and SA have been simultaneously measured. Rate constants of conversion Kss and Ko2 show the definite rate of distribution within the tissues as well as the combined output of ASA and of SA. Compartments 4 and 5 of Fig. 2(B) cannot be identified except by conjecture. As a matter of fact compartment 5 appears like a barrier lying between the stomach (compartment 1) and the plasma (compartments 2 and 3) for the transfer of ASA. Compartment 4 (see Fig. 2(A)) could be regarded as the digestive tract, the stomach not being included. TABLE 1 TABLE

1 SHOWS

MODEL OF FOUR WHEN

THE

AND

OTHER

IS INVOLVED.

OF THE

EXPONENTS

AND A THREEFOLD

SA CURVES

ALONE

(2) (3)

SA ASA SA ASA SA ASA

UP A MODEL

ARE

CALCULATED

OF FIVE COMPARTMENTS

THE RESULTS

COMPATIBLE, OF THREE

-00554 -0.0833 -0.0723 -0.0819 -0.0743 -0.0816

ON THE

ROOT Is STUDIED, WHEN

ROOTS ARE NOT COMPATIBLE.

AI

(1)

OF SETTING

RATHER

THAN

A

COMPARTMENTS.

VALUE

COMPARTMENTS THE TWO

THE OPPORTUNITY

THE

IN ANY

CHARACTERISTIC 11

-0~0554 -0.0833 -0.0723 -0.0819 -0.0743 -0.0816

BASIS

IT APPEARS ERRORS

OF A RELATION

THAT ARE

TAKEN

CASE THE PRESENCE SAMPLES

ARE

~-0.0554 -0.0833 -00723 -0.0819 -0.0743 -0.0816

INTO

A2

WHEREAS

COMPARTMENT _____-

SA

~_____-- oaO202 -000284 -0+-)0262

________

FOUR

ROOTSOF ASA

ACCOUNT,

OF A FIFTH

SHOWN. _---

11

INCLUDING

THE THREEFOLD

A2ASA -0.0209 -0.2167 -0.2349 ___~__

COMPARTMENTAL

141

MODEL

model (Fig. 2(B)) has been devised; it is the only one which has been able to fit the observations (Griessen et al., in press). Fitting test of the model

Starting from the physiological model described in Fig. 2(A), an attempt has been made to fit the data from the eight kij involved in this model and from the two volumes defining the size of the pools of ASA and of SA (the only compartments where measurements have been carried out), that is ten variables in all. The results obtained show that the first model (Fig. 2(A)) does not quite account for the experiments. Indeed, even if the data of the SA can be correctly fitted, the theoretical curve corresponding to the ASA shows an initial slope that is too steep; in addition to its peak being too far from the peak of the measured curve. According to the example shown in Table 1, it has to be determined whether it might be possible to account for the observations by adding extra kij or if a fifth compartment has to be set up. The fitting of the data has been further tested, not from the kij, but from functions such as eqn. (7), still with regard to a mode1 consisting of four compartments. In these conditions only four relationships are in keeping with such a model : (a) Four distinct roots: f(t)

= cl exp @,t) + c2 exp (A2t) + c3 exp (L,t)

+ c4 exp (A,?)

(16)

with regard to the fact that the vector S, = (1 , 0, 0,O)the above ci are compelled to the following constraints, only the second and the third compartment being fitted : c, +c,

+c,

+c,=o

(b) One double root: f(t)

= (cl f c,t) exp (AIt)

+ c3 exp (l,t)

+ c4 exp (A,t)

(17)

with c, + c3 + c4 = 0, c2 being at random since for t = 0 the expresiion with t is zero. (c) One threefold root: f(t)

withc,

= (cl + c2t + c3t2) exp (l,t)

+ c4 exp (J,t)

(18)

+c,=O.

(d) Two double roots: f(t)

= (cl + c,t) exp V,t)

+ (c3 + c4t) exp (n2t)

(19)

C. REVILLARD,

142

M. GRIESSEN, C. V. PERRIER,

J. R. SCHERRER

with cl + cJ = 0. A function comprising a fourfold root is not to be considered since cI would be equal to zero. Such a condition would be realised by a model with three compartments. With respect to the constraints, the number of variable parameters or degrees of freedom for each type of function is as follows: 7, for 4 distinct roots, 6, for I double root, 5, for 1 threefold root, 5, for 2 double roots. The curves drawn for ASA and SA have been separately fitted for three samples. The four above relationships have been used, and the x2 formulated have been compared in each case.

TABLE MULTIPLICITY ARE FOUND

OF THE

ROOTS

IN THE

FOR ASA AND SA, AS WELL TAKES

Type of relation

4 distinct roots 1 double root 1 threefold root 2 double roots

MINIMAL

AS FOR THE THREE

INTO ACCOUNT

Number of variable parameters (degrees of freedom) : :

2

POLYNOMIAL

A MULTIPLICITY

OF MATRIX

SAMPLES

11 10 6 10

K.

THE

BEST

THE TYPE

OF THE ROOTS OF

Case 1 x* SA

WHEN

2

OR

x* SA

OF X2

3

Case 3

Case 2

x*ASA

VALUES

OF RELATIONSHIP

x* ASA

x*

SA

x*ASA

27 :i 23 47

it: 12 13

E 25 31

:f ::, 23

:;:

Paradoxically, the best x2 (Table 2) correspond to the fits for which the number of degrees of freedom was smallest (I double root or 1 threefold root). Multiple roots have to be taken into account to analyse such data. On the other hand, the analysis of the exponents li connected with the separate fits of the ASA and SA data has shown that when three of them are compatible, the fourth is not (Table 2). Thus, in the first analysis it appears that it was the number of compartments which was insufficient and not the number of kij. That is why the present model had to take into account a fifth compartment which is suggested in Fig. 2(B) and which is physiologically plausible. The number of kij has not exceeded nine, and the number of compartments of the present mode1 has been set to five. Indeed, for any aspirin studied the fit of the data has been systematically satisfying, not only from one individual to another but also from one preparation of aspirin to another (Fig. 3).

143

COMPARTMENTALMODEL

ASA

SA mmoles/l

mmoles/l

0.02 0.01 0 , 0

200 400 600

40

80

80

160 240

0061 0.02

0.01

C 0

100 200 300 400 500 600 Minutes

0

400 600

Minutes

Fig. 3. The model of Fig. 2(B) has been fitted plausibly on the basis of synchronous observations respectively of ASA and SA: four preparations of aspirin commercially available have been tested successively in four patients. Figure 3 depicts the fits for ASA and SA of the four prepvations of aspirin in one of the patients. The fits shown in Fig. 3 are comparable with those concemm$ the other patients.

144

C. REVILLARD,

M. GRIESSEN, C. V. PERRIER, J. R. SCHERRER

Thus it can be assumed that the formulated model may account for the available data. Nevertheless, considering the efficiency and usefulness of our fitting method and the ill-conditioned matrices in respect of particular sets of data, the proposed model is not necessarily supposed to have multiple roots. Indeed the x2 criterion, as used here, cannot have an absolute meaning. Finally, the convergence has been checked so that every time it reaches the same final values, starting from different initial values of kij. For such a test the same initial value for all the k,, has been used. The test has been repeated three times with values equal to 5 * lo-‘, 5 * 10e2, and 5 * IO-‘. The three fits performed have been converged to the same final values of the k,j. However it is not p’ossible to demonstrate whether the field of the possible kij has been explored in such a way that the convergence reaches a true minimum and that the relationship found for each sample is really unique.

CONCLUFjIONS

The objective of the present work has been to adapt calculation methods using a compartmental model in order to fit experimental observations even in the case of initial estimates very distant from their final value. Such a versatile investigation has resulted in formulating and testing the most likely plausibility of a five compartment model for four preparations of aspirin and also for the patients tested. Such a five compartment model has appeared more plausible than the four compartment model which was the starting point of the physiological hypotheses of the present study. Therefore the hypotheses concerning the behaviour of aspirin within the human body have to be formulated in another way, that is by considering an additional fifth compartment to the model. The complexity of the physiological transfer mechanism also corresponds to such an assumption. The pharmaco-dynamic aspect of aspirin, namely the linearity of the fitted model, has been discussed previously (Griessen, et al., in press). In addition, the occurrence of multiple roots in the fitting of the compartmental model has the following consequence: the number of exponential components necessary to fit the observations has to be dissociated from the most plausible number of compartments contained by the system. As a matter of fact, the number of exponential components may be smaller than the number of compartments. Furthermore, the method used has resulted in a clear-cut reduction of the computer calculation time. The integration method used in the present work differs from the integration methods commonly used in the field of compartmental analysis. Computer calculation time has been all the more clearly reduced as the observations have been more numerous.

COMPARTMENTAL MODEL

145

APPENDIX

Assuming a square matrix M with a rank d. The characteristic polynomial of this matrix is: A(a) = (a - a$‘l(a

- a#9

. . . (a - a,)pr

(20)

with : d=

ZPi

where ai are the roots of this polynomial or the eigenvalues of the matrix M, pi is the multiplicity order of the root ai and r is the number of roots distinct from A(a). A(a) is a polynomial which annuls the matrix M, that is to say: A(M) = 0

(21)

The polynomial annulator Y(a) with the smallest degree is called the minimal polynomial of M. It can be written: ‘Y(a) = (a -

a l)ml(a -

aJml

. . . (a - a,)mr

(22)

and Y(M) = 0

(23)

m, is the multiplicity order of the root ai and m = C;=, m, the degree of Y(a). This polynomial is unique and m < d. The number r of distinct roots of the polynomials A(a) and Y(a) is the same. It follows that a matrix function can be written (Gantmacher, 1959~): f(M) = LJWZ,,

I +f'Wze.2+ **-+fm~-'kJze,me (24)

The matrices Z,,{ are square matrices of order d. They constitute the matrix M.

ACKNOWLEDGEMENT

The authors gratefully acknowledge the assistance and advice of Dr Mones Berman (National Institute of Health, USA).

REFERENCES

BERMAN,M. and SCHOENFELD,R. L., Invariants in experimental data on linear kinetics and the formulation of models, J. Appf. Phys., 27 (1956) pp. 1361-70. BERMAN,M., SHAHN, E. and Wans, M. F., The routine fitting of kinetic data to models: A mathematiall formalism for digital computers. Bfophys.J., 2 (1962) pp. 289-315.

146

C. REVILLARD, M. GRIESSEN, C. V. PERRIER, J. R. SCHERRER

BERMAN,M. an+ WEISS, .M. F., Users’. manual for SAAM (Simulation, Analysis and Modeling), 7$66by Nattonal Instrtute of Arthrttts and Metabolic Diseases, N.I.H., Bethesda, Maryland, BROWNELL,G. L., BERMAN,M. and ROBERTSON J. S Nomenclature for tracer kinetics. A Report to the International Commission on Radiolo&al Units and Measurements by the Task Group on Tracer Kinetics, 20 April, 1967, BNL-11819. BULIRSCH,R. and STOER,J., Numerical treatment of ordinary differential equations by extrapolation methods, Num. Math., 8 (1966) p. 1. DAmoN, W. C., Variance algorithm for minimisation, Comp. J., 10 (1968) p. 406. FRANCIS,J. C. F., The QR transformation: An analogue to the LR transformation, Comput. J.. 4 (1961/62) pp. 265-71 and 332-3. GANTMACHER, F. R., Matrix theory, Chelsea Publishing Co., New York, 1959a. (First American edition.) GANTMACHER, F. R., Application of the theory of matrices. Interscience. Publishing Co., New York, 19596. (First American edition.) GRIESSEN,M., REVILLARD,C., SCHERRER,3. R., SITAVANC,L. and PERRIER,C. V., Absorption kinetics of aspirins in humans, in press. HAMMING, R. W., Stable predictor-corrector methods for ordinary differential equation, J. Assoc. Comp. Mach., 6 (1959) p. 37. HEARON,J. Z., Interpretation of tracer data, Biophys. J., 9 (1969) pp. 1363-70. JACQUEZ,J. A., Compartmental analysis in biology and medicine, Elsevier Publishing Co., Amsterdam, London, New York, 1972. JAMES, F., Monte Carlo for particle physicists, Section 6.1: Proceedings of 1968 Herceg-Novi School (TO be published by Gordon and Breach, London). JAMES,F. and Roes, M., MINUIT, a package of programs to minimise a function of n variables, compute the covariance matrix, and find the true errors. CERN Computer 6000 Series, Program Library, D 506/D 516, 1971. JULIUS, R. S., The sensitivity of exponentials and other curves to their parameters, Computers and Biomed. Res.. 5 (1972) DV. 473-8. LEWALLEN,C. G., BERMAN,M: and RALL, J. E., Studies of iodoalbumin metabolism. I. A mathematical approach to the kinetics, J. Clin. Invest., 38 (1959) pp. 66-87. MARTIN, D. W., Runge-Kutta methods for integrating differential equations on high speed computers, Computer J., l(3) (1958) pp. 119. NEDF~~;_ A. and MEAD, R., A stmplex method for function minimisation, Comp. J., 7 (1967) RE~CIGNO,A. and SEGRE, G., Drug and tracer kinetics, Blaisdell Publishing Co., Waltham, Massachusetts: Toronto, London, 1966. (First American edition.) ROBERTSON,J. S., Theory and use of tracers determining transfer rates in biological systems, Physi 1. Rev., 37 (1957) pp. 133-54. g K, H. H., An automatic method for finding the greatest or least value of a function, ROSENBRO Comp. J., 3 (1960) p. 175. ROWLAND,M., RIEGELMAN,S., HARRIS, P. A. and SHOLKOFF,S., Absorption kinetics of aspirin in man following oral administration of an aqueous solution, J. of Pharmaceutical Sciences, 61-3 (1972) pp. 379-85. SHEPPEY,G. C., Minimisation and curve fitting, CERN Yellow Report, 68-5, 1968.