CXT-MAIN: a software package for determination of the analytical form of the pharmacokinetic system weighting function

CXT-MAIN: a software package for determination of the analytical form of the pharmacokinetic system weighting function

ELSEVIER Computer Methods and Prog:rams in Biomedicine 51 (1996) 183-192 CXT-MAIN: a software package for determination of analytical form of the ph...

785KB Sizes 1 Downloads 70 Views

ELSEVIER

Computer Methods and Prog:rams in Biomedicine 51 (1996) 183-192

CXT-MAIN: a software package for determination of analytical form of the pharmacokinetic system weighting function Ladislav Dedik”, Mzh-ia 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; e-mail: [email protected] 0169-2607/96/$1.5.00 0 1996 Elsevier Science Ireland Ltd. All rights reserved PII 0169-2607(96)01770-l

model; Drug bioavailability

184

L. Dedjk,

hf. duui.fovri

/ Computer

Methods

1. Lntroduction Each dynamic system has its inherent model-independent 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 s-domain 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 de-mom&ate a procedure for the determination of

and Programs

i72

Biomedicine

51 (1996)

183-192

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 t-domain 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 s-domain, 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 s-domain 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) 183-192

:-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 n-order and m-order 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 partial-fraction 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 CXT-MAIN. 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)

183-192

(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 s-domain 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 cardio-pulmonary 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)

183-192

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.70x10-3t

+ (1.18 x

10-4)

e-4.09x

+ (2.89 x

10-3),-8.lOxlOP3f~

lo-2/

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 45Ca-labelled 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)

183-192

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).

b-r 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-“) e-1.76x lo-% .+ (3.74 x 10-3) e-6.52x 10-3t -+ e-7.l9

x IO-32

[(7.27 x 10-4) cos(1.84 x 10-2t)

- (3.21 x 10s4) sin(1.84 x 10-2t)])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) 183-192 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 10-l 2.97 2.16 1.65 8.86 x 10-l 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 + ei-0.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) 183-192

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 double-site absorption model in the study [16], used to explain the double-peak 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

183-192

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)

Qcp-H,,(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. Veng-Pedersen, An algorithm and computer program for deconvolution in linear pharmacokinetics, .J.Pharmacokin. Biopharm. 8 (1980) 463-481. [3] W.R. Gillespie and P. Veng-Pedersen, 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) 289-307. and M. DuriSovC, CXT-A programme for 141 L. Dedik analysis of linear dynamic systems in the frequency domain Int. J. Bio-Med. Comput. 39 (1995) 231-241. and L. Bousquet, Time-delay for two [51 Y. Plusquellec compartment models used for study of enterohepatic circulation of drugs, IEEE Trans. Biomed. Eng. 31 (1984) 469--4X. 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) t-93. 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) 504-510. 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)

183-19,

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. AC-4 (1959) 37-43. 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) 531-536. Y. Plusquellec, G. Campistron, S. Staveris et al., A double-peak phenomenon in the pharmacokinetics of veralipride after oral administration: a double site model for drug absorption, J. Pharmacokin. Biopharm. 15 (1987) 225-239. 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. 27-96 (Academic Press. New York, NY, 1976).