ELSEVIER
Computer Methods and Prog:rams in Biomedicine 51 (1996) 183192
CXTMAIN: a software package for determination of analytical form of the pharmacokinetic system weighting function Ladislav Dedik”, Mzhia lhri,50vBb~* “Faculty of Mechanical blnstitute of Experimental
Engineering, Pharmacology,
Slovak Technical University, 812 31 Bratislava, Slovak Rqpublic Slovak Academy of Sciences, 842 16 Bratislava, Slovak Republic
Received 14 February 1996; revised 4 June 1996; accepted 26 June 1996
Abstract A new procedure specific for the determination of the analytical form of the model weighting function of a complex multicomponent pharmacokinetic system with or without a shunt and time delays is described. The procedure is based on the theory of linear dynamic systems and on a circulatory pharmacokinetic model of the living body. The model transfer function of the system under study was obtained by the frequency response method in the form of the ratio of two frequency dependent polynomials. Subsequently, the technique of the partial fraction inversion was employed to determine the analytical form of the model weighting function. Two examples from bioavailability studies in pharmacokinetics are given. The first example presents two estimates of the model weighting function of a pharmacokinetic system obtained by the new procedure and by a polyexponential deconvolution method. To compare these results, two models of the measured system output were determined using the two estimates of the model weighting function, the actual system input and a convolution method. The model weighting function obtained by the new pnocedure yielded a better model approximation of the output data than that obtained by the polyexponential deconvolution method. The second example, usin ‘5 the new procedure, presents the determination of the model weighting function of such a system that the deconvolution methods, commonly used in pharmacolcinetics, cannot be applied to. Keywords:
Linear
system; Weighting
function;
Pharmacokinetics;
Circulatory
* Corresponding author. Tel./fax: + 42 7 375928; email:
[email protected] 01692607/96/$1.5.00 0 1996 Elsevier Science Ireland Ltd. All rights reserved PII 01692607(96)01770l
model; Drug bioavailability
184
L. Dedjk,
hf. duui.fovri
/ Computer
Methods
1. Lntroduction Each dynamic system has its inherent modelindependent properties. The identifiable system properties can be determined by measurements on the system. Based on these measurements, properties of the dynamic system can be described by static and dynamic characteristic functions. The static characteristic function expresses the relationship between the system input and output in steady state. For a linear dynamic system, it is a straight line with a slope called the system gain. The three dynamic characteristic functions of a linear dynamic system, i.e. the impulse characteristic function, the transient characteristic function and the frequency characteristic function are the system responses to three specific system inputs: the Dirac delta pulse, the unit step function and the harmonic signal with unit amplitude, respectively. If the system input has the same form yet not the same magnitude as does either the Dirac delta pulse, the unit step function, or the harmonic signal with unit amplitude, the system response is not the system characteristic function but it only reflects the system impulse function, the system transient function or the system frequency response. The corresponding system charact,eristic Sunctions can be obtained on the basis of these functions by dividing them by the magnitude of the input signals. The system impulse characteristic function, frequently called the system weighting function, is regarded the most important characteristic function of the system. It is a unique characteristic function of the system. since its form in the Laplace sdomain is equal to that of the system transfer function [l]. In practice, this function is usually determined by direct measurement on the system or calculated by the deconvolution method in the tim’e domain. The literature on the deconvolution methods, numerical and analytical, is extensive. The analytical deconvolution methods are mostly devoted to relatively simple models [2,3] which do not contain shunt and time delays and whose transfer functions do not contain complex poles El]. The purpose of this communication is to demom&ate a procedure for the determination of
and Programs
i72
Biomedicine
51 (1996)
183192
the analytical form of the model weighting function of a complex multicomponent pharmacokinetic system, with or without shunt and time delays, i.e. a system whose model transfer function contains real and/or complex poles.
2. Theoretical E,q. (1) is the general linear tdomain model of all linear or linearised complex multicomponent systems, with or without shunt [4] and time delays
151, dkX( t)
x(t) + f bkT k=l
=
(1)
where t is time, X(t) and Y(t) are the output and input of the system, respectively, G is the system gain and the symbols a,,. . ., n,,, b,, . . ., b,,, are the model parameters. If the transfer function [I] of the system is defined by Eq. (2)
H(s)=X(s) __ m 1’ in the Laplace sdomain, rewritten into Eq. (3) H(s)
= G
(2) then Eq. (1) can be
n,+a,s+n,s’+... +n,s” 1 + b,s + b,s2 + . + b,,s’“’
(3)
where s is the Laplace operator and a, C=Z1 [1,6]. Eq. (3) is the general linear sdomain model of the transfer functions of all linear or linearised complex multicomponent dynamic systems, with or without shunt and time delays [4,6&B]. For the system input in the form of a Dirac delta pulse s(t), where m 6(r)dt= 1, (4) s0 and L{&(t)}
=
‘= d(t) e“‘dt = 1, (5) s0 the corresponding system output is the system impulse characteristic function itself, or in other words, the system weighting function. It follows the:n that in the general case the system weighting function W(t) is given by Eq. (6)
L. Dedik, M. hiSovb
W(t) = L ~ ‘{H(s))
=&
/ Computer Methods ana’ Programs in Biomedicine 51 (1996) 183192
:Ix H(s)e”’ ds. sc Iso
The integration takes place in the Laplace main. c > 0 and i is the imaginary unit [9]. Eq. (3) can be rewritten into Eq. (7)
(6)
do
where A(s) and B(s) are norder and morder polynomials, respectively [6,7]. For the specific case YU= n, indicating a system shunt [4], Eq. (3) can be rearranged into the form of Eq. (8)
where C(s) ==A,(s) f+(s), II A,(s) = A(s)  u,f,
(10)
B,(s) = B(s)  ~5,s”.
(11)
The general analytical form of the weighting function of the complex multicomponent system, with or without shunt and time delays, can be expressed in the form of Eq. (12) using the technique of the partialfraction inversion [9]
?Pk,j
sin(@kt>l
,
1112)
where  y = 1 if  pi and complex nomial, ~
rn = Iz, or y = 0 if n < 172, p2 are numbers of real roots &
PZ Yl,/‘+ C r2,j k=l
j=l
=
n1S
 yl;,j and clk,j k ip,c,j are real and complex coeffi
cients, respectively of the Heaviside expansion [9,101.
3. Procedure description
A(s) H(s) := GB(s)’

185
(13)
The procedure designed for the determination of the model weighting function in the analytical form consisted of the following steps: (1) The program CXT (CompleX Tools for Linear Dynamic System Analysis) based on the frequency response method [ 1.4,6,7,11,12] was employed to model the system under study in the frequency domain: (a) The frequency response of the system was calculatecl as the ratio of the Fourier transforms of the system output and input signal. The Fourier transforms were determined as sums of the Fourier transforms of straight lines approximating the measured drug concentrations between two adjacent sampling points and the Fourier transforms of exponential functions approximating the concentration time courses in time intervals from the last sampling time to infinity. (b) Point estimates of the parameters a,,. . ., a,, k,..., b,,, of the model transfer function in the form of Eq. (3) were obtained using the Levy method [13], implemented in the CXT program. To decide on optimal numerator and denominator polynomial degrees n and 171 the Complex Criterion (CC) and the Akaike information criterion (AIC) in the complex and time domain, respectively were applied [7,181. (2) The program
ROOTS, based on the GaussNewton algorithm, was used to determine real and imaginary (nonmultiple and multiple) roots Ak and Qk & i$k of the denominator polynomial B(s) of the model transfer function in the form of Eq. (3). based on the Heav(3) The program ANALYT, iside expansion theorem, was used to determine the analytical form of the model weighting function in the form of Eq. (12).
186
L. Dedik,
M.
duris’ovb
/ Computer
Methods
The coefficients of the Heaviside expansion theorem ~1~.j and ak. j + ia,,/ were determined individmally for the real and imaginary (nonmultiple and multiple) roots of the denominator polynomial B(s) of the model transfer function in the form of Eq. (3). The programs ROOTS and ANALYT can cooperate in the environment of the CXT program, using the control program CXTMAIN. The version ‘96 of the CXT program, including its environment, is available at: http://www.cpb. uokhsc.edu/pkin/pkin.html, option: software. It can be run on IBM compatible 386 and 486based machines under the operating system MSDOS without a math coprocessor,
and Programs
in Biomedicine
51 (1996)
183192
(15) Eq. (16) defines the system between the drug input into) the blood circulation after extravascular administration Y,,(s) and the drug output measured in a peripheral loop of the blood circulation J&.h) (161
Eq. (17) can be derived for the transfer function H,,(s) on the basis of the. model shown in Fig. 1 using the method of the Laplace block algebra (see derivation of Eq. (17) in Appendix A)
B 4
4. Procedure application
‘Yiv In pharmacokinetics, the drug bioavailability in the blood circulation after extravascular administration is defined in two ways: as the extent and the rate of bioavailability, i.e. the amount of the drug and the rate at which the drug reaches the blood circulation [14]. Fig. 1 shows a circulatory model of the living body. The equations summarized below can be derived on the basis of this model. The system describing the drug input into the blood circulation after extravascular administration can be defined in the Laplace sdomain by Eq. (14) (14) where Ha,,(s), X,,(s) and Y,,(s) are the system transfer function, output and input, respectively. X,,(s) simultaneously represents the drug input into the blood circulation after extravascular administration. This input is not available for measurement. Eq. (15) defines the system between the drug input into the blood circulation after intravascular administradon YJs) and the drug output measured in a peripheral loop of the blood circulation &,ds)
I
X PN X p,ex
FJ
Fig. 1. A circulatory model of the living body. Y, X, H and Q are the inputs, outputs, transfer functions and blood flows, respectievly; Indexes cp and ab correspond to the cardiopulmonary and absorption system, respectively. Indexes r and p correspond to the system of remaining body tissues and the systnm of a peripheral loop of the blood circulation, respectively. Indexes iv and ex correspond to the intravascular and extravascular inputs, respectively.
L. Dedik,
Hex(s) ‘= Hiv(s) ’ Hab(S).
M. duris”owi
1 Computer
Methods
and Programs
in Biomedicine
51 (1996)
183192
187
(17)
The transfer function H&S) can be than expressed in the form of Eq. (18), using the definitions given by Eqs. (15) and (16)
(18) For the special inputs riv(s> = K’ yex(s>,
(19)
where K is a constant K > 0, Eq. (18) can be rewritten in the simple form of Eq. (20) H&(S) = 3c.Ax
($1
(20)
xp,iv(s)’
Eqs. (18) and (20) enable us to estimate the transfer function H&S) on the basis of the measured functions I,,,,, Y,, Xp,iv and Y,,. Using the known transfer function Ha&), the drug input into the circulatory system after extravascular ,administration X,,(s) can be determined according to Eq. (21) x&)
= K,(s)
Yex(s).
(21)
The gain of the model transfer function approaches the extent of extravascular bioavailability of the drug studied. The model weighting function corresponding to the model transfer function approaches the rate of the drug extravascular bioavailability. In the specific case given by Eq. (19), the model weighting function can be estimated only on the basis of the two functions, i.e. X,,,,(s) and X,,,(S), as is commonly used in the deconvolution methods in pharmacokinetics. The arguments presented are based on the precondition that the disposition kinetics of the drug is linear.
Fig. 2. Schematic representation The corresponding time course %a (0); (b). (Example 1).
of 45Ca i.v. input in rats (a). of the plasma concentration of
forte 500 mg (open circles in Fig. 3b). The time courses were approximated by polyexponential functions. Since forms of bot,h 45Ca inputs (Dose. s(t)) were equal (Fig. 2a, Fig. 3a), the extent of calcium bioavailability was estimated as the ratio of the areas under the polyexponential functions. The estimation yielded the value of 0.72. The rate of calcium bioavailability was determined by an analytical deconvolution method based on the polyexponential functions. The estimation yielded Eq. (22) (thin line in Fig. 4) Wpolyexp(f ) = (6.04 x 10 ~4),~1.70x103t
+ (1.18 x
104)
e4.09x
+ (2.89 x
103),8.lOxlOP3f~
lo2/
P2>
= a
5 D 8
4.2, Examples 4.2.1. Example
1
In our previous study [15] oral bioavailability of calcium in rats was estimated on the basis of 45Ca plasma concentration time courses after i.v. administration of the solution 45CaC1, (black points in Fig. 2b) and after p.o. administration of 45Calabelled prep aration of Calcium Sand02
Time
Time lninl
Fig. 3. Schematic representation of 45Ca p.o. input in rats (a). The corresponding time course of the plasma concentration of 45Ca (0), (b). The model system output, obtained on the basis of the model weighting function from Eq, (22) (thin line); the model system output, obtained on the basis of the model weighting function from Eq. (23) (thick fine). (Example 1).
18s
L. Dedik,
M. duriiovb
/ Computer
Methods
and Programs
in Biomedicine
51 (1996)
183192
b
Time (h)
:F: ~0.000 ‘5
i
0
60
120
240
180
300
360
3
Time
(mid
Fig. 4. The two estimates of the model weighting function approaching the rate of calcium oral bioavailability in rats. The estimate from Eq. (22) (thin line). The estimate from Eq. (23) (thick line). (Example 1).
br our current study, the frequency response method was used to estimate the model transfer function of the system describing oral calcium bioavaiiability in rats (Eq. (20)), on the basis of the data published in the above mentioned study. The model transfer function was obtained in the form of Eq. (3) using the CXT program. The point estimates of the model parameters are given in Table 1. The estimation of the gain of the model transfer function, approaching the extent of calcium bioavailability, yielded the value of 0.67. The parameters of the model transfer function (Table 1) and the programs ROOTS and ANALYT were used to estimate the model weighting function on the basis of Eq. (12). This model weighting function, approaching the rate of oral calcium bioavailability, is given by Eq. (23)
Time 04
Fig. 5. Veralipride input during 0.5 h i.v. infusion in ‘B’ volunteer, (a). The corresponding time course of the plasma concentration of veralipride (O), (b). (Example 2). Wcxt(t>
= 0.67{(7.58 x IO“) e1.76x lo% .+ (3.74 x 103) e6.52x 103t + e7.l9
x IO32
[(7.27 x 104) cos(1.84 x 102t)
 (3.21 x 10s4) sin(1.84 x 102t)])7
(23)
and shown as the thick line in Fig. 4. 4.2.2. Example 2 F’harmacokinetics of different dosage forms of veralipride in volunteers was studied [16]. To estimate the oral bioavailability of veralipride, the data corresponding to ‘B’ volunteer of that trial were reevaluated in our current study. The inputs of veralipride corresponding to 30 min i.v. infusion and a p.o. dose in the form Dose. S(t) are shown in Fig. 5a and Fig. 6a, respectively. The respective veralipride plasma concentrations are
Table 1 Point estimates of the model parameters of the system given by Eq. (20) and describing calcium oral bioavailability in rats Model paramerer
Estimate
Dimension
a0 al a2
1 4.29 x lo2 1.83 x lo4 l.17x106 7.58 x IO2 1.16x lo5 5.08 x lo6 2.24 x lo*
min min’ min3 min min* min3 mm9
a3 b, b, b, ba,
b

00
Time
6
12
Time
d 18
24
(h)
Fig. ‘6. Schematic representation of veralipride po. input in ‘B’ volunteer. (a) The corresponding time course of the plasma concentration of veralipride (0), (b). (Example 2).
L. Dedik, M. LjuriSovd 1 Computer Methods and Programs in Biomedicine 51 (1996) 183192 Table 2 Point estimates of the model parameters of the system given by Eq. (18) and describing veralipride oral bioavailability in ‘B’ volunte’er Model
a0 aI a2 a3 a4 b, b, b, 4 b,
parameter
Estimate
Dimension
9.91 x 10l 2.97 2.16 1.65 8.86 x 10l 8.83 1.29 1.44 5.04 1.85
h h2 h3 h4 h h2 h3 h4 h5
illustrated as black points and open circles in Fig. 5b and Fig. 6b, respectively. S&e the two veralipride inputs failed to meet $e condition expressed by Eq. (19) (Fig. 5a, Fig. 6a), the transfer function given by Eq. (18) was used to determine the extent and rate of (oral bioavailability of veralipride in ‘B’ volunteer. The point estimates of model parameters are given in Table 2. The estimation of the gain of the model transfer function, approaching the extent of veralipride bioavailability, yielded the value of 0.86. The parameters of the model transfer function (Table 2) and the programs ROOTS and ANALYT were used to estimate the model weighting function, on the Esasis of Eq. (12). This model weighting function approached the rate of oral veralipride bioavailability. It is given by Eq. (24)
189
5. Discussion
The mathematical model of the linear system expressed in the form of the transfer function is frequently called the disposition function in pharmacokinetics [ 171. The system weighting function is the solution of the disposition function for the specific input function, i.e. the Dirac delta pulse s(t). In the great majority of pharmacokinetic studies the disposition, function is used in such a form as having only real roots and thus not modelling the system time delay [2,3,17]. In such cases also the corresponding analytical form of the weighting function contains only real roots and is given by the first term of Eq. (12)
where y usually equals zero. The whole general analytical form of the weighting function described by Eq. (12) allows the estimation of the model weighting function also for complex multicomponent systems which contain shunt and time delays in their structure [8]. If the system under study contains time delays, only highorder A(s) and B(s) polynomials in Eq. (7) can optimally approximate this system [8] and then both the disposition and weighting function of this system contain complex conjugate roots.
= 0.8612(0.1043 ec0.1364t + ei0.45s8’[(4.608 x lo “) cos(0.7541t) + 0.218 1 sin(0.7541 t)] + e 0.*3s4r[0.3687 cos(2.0956t)  0.1199 sin(2.0956t)]), and shown as a solid line in Fig. 7.
Time
(24)
(hl
Fig. 7. The estimate of the model weighting function approaching the rate of veralipride oral bioavaiiability in ‘B’ volunteer. (Example 2).
190
L. Dedik, M. DuriSovci / Computer Methods and Programs in Biomedicine 51 (1996) 183192
Both systems analyzed in our study were complex. II. was thus highly probable that they would contain time delays. In example 1 we estimated calcium oral bioavailability using the different methods and models. The results obtained were very similar. To compare these results, models of the measured system output were calculated using the two model weighting functions and a numerical convolution method. The output models are shown as a thin line (the model weighting function based on Eq. (22)) and thick line (the model weighting function based on Eq. (23)), respectively in Fig. 3b. The values of the AIC criterion corresponding to the descriptions of the measured output data by the model outputs based on the estimates of the weighting function in the form of Eq. (22) and Eq. (23) were  72.07 and  75.55, respectively. These two descriptions along with the AIC values demonstrated that the model weighting function in the form of Eq. (23) yielded a better prediction of the output data than that in the form of Eq. (22). The essential difference between the two estimates of the weighting function is that Eq. (22) does not contain the trigonometric terms while Eq. (23) does. The presence of the trigonometric terms in Eq. (23) confirmed the hypothesis that the system under study was complex and contained time delays [8,12]. In example 2 we presented the determination of the model weighting function of such a system that the deconvolution methods, commonly used in the pharmacokinetic literature, cannot be applied to because of the different inputs (Fig. 5a, Fig. 6a). The mode1 weighting function in Fig. 7 shows two significant peaks at 0 h and 2.5 h. The presence of these peaks is in agreement with the doublesite absorption model in the study [16], used to explain the doublepeak phenomenon in the pharmacokinetics of veralipride after oral administration in human volunteers.
with or without shunt and time delays; ~ the analytical form of the model weighting function for a complex multicomponent system after an oscillating output;  the analytical form of the model weighting function of a multicomponent system even in the case when a deconvolution method cannot be used since the system inputs fail to meet the condition expressed by Eq. (19), e.g. example 2.
Acknowledgements
The authors express their thanks to Prof. Y. Plusquellec (Universite P. Sabatier, Toulouse, France) for providing the data describing veralipride kinetics and to the Journal reviewers for constructive comments on the manuscript. This work was supported in part by Grant 1021/94 from the Slovak Grant Agency.
Appendix A: Derivation of Eq. (17)
Basic equations: X,(4
= K,(s).
X(s),
Qc, = (zab+ Qr + Qp, x,,(s)
6. Concluding remarks
The procedure presented allows the user to calculate:  the analytical form of the model weighting function of a complex multicomponent system
+ Qr . K(s). X,(s) + Qp.f&(s). X,6> (A7)
L. Dedik,
M. buriiovd
1 Computer
Methods
and Programs
in Biomedicine
[email protected]~L,.&,(s)+ Q~.Hr(s>+ Q;HpW (AS) + fc&) Ye,(s)+ Yj”@)l. Intravascular
administration
KY($) f O;
51 (1996)
191
183192
H ($1 zE. (Qab.Ha&)
Q“P
(A15)
of the drug:
yex(s) = 0
H,,(s). fL&) QCP
=
l _ H,,(s)
e(Qad&,,(S)
+
Q, *f&(s) + Q, .H,(s))
1
+
Q, .Hr(s)
+ Qp .N,(s))
CP
>
6416)
649) x,,(s)
_
Y,,(s)
QcpH,,(s).
( Qa,,. Ha,,(s) + Q, . H,(.F) + Qp.H,(s))’
(Al71
(A 10)
xp ex(S) H,,(s) = b us> ’
6418)
Hex(s) = x,,(s)
_
H,(WC,(WC&)
US)
Q,,H,,(s).(Q,,.H,,(s)+Q,.Hr(s)fQ;H,(s)) K,(s)
Qc,
(Qa,,. fL&)
l’%,(s).
+
Qr . H,.(s)
+
(Al
Hiv($>
=
xp y
(Al91
Qp. H,(d)’
iAs)
1)
Dividing
Eq. (A13) by Eq. (A19) yields 6420)
G412)
IV(S)
It follows then that
HLv(S)= ~,(~~~H,p(~)
Qcr,
K,(s).
(QA. Ha&)
+
Qr . H,.(s)
+
Qp.H,(,s))’ (A131
ExtravascuIar YJS)  0;
administration L(s)
of the drug:
#0
.(Q!at,.Ija&)+ Q;K(s) + Q;HpW + Ha&). Y&)1,
(Al4)
References [l] M.G. Singh, Systems and Control Encyclopedia, Theory, Technology and Application (Pergamon Press, Oxford UK, 1987). [2] P. VengPedersen, An algorithm and computer program for deconvolution in linear pharmacokinetics, .J.Pharmacokin. Biopharm. 8 (1980) 463481. [3] W.R. Gillespie and P. VengPedersen, A polyexponential deconvolution method. Evaluation of the “gastrointestinal bioavailability” and mean in viva dissolution time of
192
L. Dedik,
M. &Gowi
/ Cowputer
Methods
some ibiuprofen dosage forms on appropriate constraints on the initial input response when applying deconvolution, J. Pharmacokin. Biopharm. 13 (1985) 289307. and M. DuriSovC, CXTA programme for 141 L. Dedik analysis of linear dynamic systems in the frequency domain Int. J. BioMed. Comput. 39 (1995) 231241. and L. Bousquet, Timedelay for two [51 Y. Plusquellec compartment models used for study of enterohepatic circulation of drugs, IEEE Trans. Biomed. Eng. 31 (1984) 4694X. of Linear 1’51J. Schoukens and R. Pintelon, Identification Systems. A Practical Guide for Accurate Modelling (Pergamon Press, London, UK, 1991). Frequency response method L71 L. Dedik and M. DuriSova, in pharmacokinetics, J. Pharmacokin. Biopharm. 22 (1994) 2933307. approchee dune function PI H. Pad& Sur la representation par des fractions rationelles, Ann. de I’Ecole Normale 9 (1892) t93. de Mathematiques a l’usage des [91 A. Angot, Complements Ingtnieurs de I’Electrotechnique et des Telecommunications (MASSON et Cie, Paris France, 1952). Operators in mathematical physics, Proc. !lOl 0. Heaviside, R. Sot. A52 (1893) 504510. study of human illI M. &iriSova and L. Dedik, Comparative pentacaine pharmacokinetics in time and frequency do
and Progmm
[12]
[13] [14] [15]
[16]
[17] [18]
in Biomedicine
51 (1996)
18319,
main, Meth. Find. Exp. Clin. Pharmacol. 16 (1994) 219232. M. Durisova, L. Dedik and M. Balan, Building a structured model of a complex pharmacokinetic system with time delay, Bull. Math. Biol. 57 (1995) 7677808. EC. Levy, Complex curve fitting, IRE Trans. Automat. Contr. AC4 (1959) 3743. J.G. Wagner, Fundamentals of Clinical Pharmacokinetics (Drug Intelligence Publications, Hamilton IL, 1975). M. burilova, T. Trnovec, S. Bezek, Z. Kallay and M. Zemanek, The rate and extent of calcium bioavailability from two oral dosage forms in rats, Gen. Physiol. Biophys. 4 (1985) 531536. Y. Plusquellec, G. Campistron, S. Staveris et al., A doublepeak phenomenon in the pharmacokinetics of veralipride after oral administration: a double site model for drug absorption, J. Pharmacokin. Biopharm. 15 (1987) 225239. M. Gibaldi and D. Perrier, Pharmacokinetics (Marcel Dekker, New York, NY, 1982). H. Akaike, Canonical correlation analysis of time series and the use of an information criterion, in: System Identification Advances and Case Studies. eds. A.F. Mehra and D.G. Lainiotis, pp. 2796 (Academic Press. New York, NY, 1976).