Method to evaluate the systemic exposure in toxicological and pharmacological studies

Abstract

A method of estimating the average exposure within a population of subjects to a pharmacological substance, administered according to a given protocol, by estimating the area under the concentration curve (AUC), characterised in that the area under the population concentration curve is estimated by the steps of obtaining measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken; building a hierarchical stochastic model composed of the population and of the individual levels; determining the posterior probability distribution of the average AUC from the sample data and hence the average population AUC. The exposure of an individual and its precision within the sample of individuals can be determined by determining the posterior probability distribution of the individual AUC from the sample data. The posterior probability model may be obtained by using a Markov Chain Monte Carlo algorithm, such as Gibbs Sampling.

Description

[0001] The present invention relates to pharmacokinetic studies which are used to evaluate the systemic exposure of animals or subjects to a pharmacological substance. It also relates to the correlation of the descriptor of such exposure to therapeutic and toxic effects (pharmacodynamics) of the drug.

[0002] One of the most important measures of the exposure of animals or subjects to a pharmacological substance is the area under the plasma concentration-time curve (AUC). The area under the curve is the integral of the concentration of the drug in plasma with respect to time. In practice it is rarely possible to take continuous measurements of the plasma concentration of the drug and instead samples are collected from the subjects at times distributed over the exposure time. In cases were it is possible to collect a sufficient number of samples from the subject the AUC and hence the exposure of the drug can easily be determined by numerical integration techniques such as the linear trapezoidal rule. In the following text, the term AUC will indicate the area under the plasma concentration-time curve calculated up to some fixed time (AUC0-t). However, all the procedures can be extended to the total exposure, i.e., the AUC calculated up to infinite time, by means of a suitable mathematical model describing the decline of the plasma concentration-time curve.

[0003] Furthermore, because population values, such as population AUC and population mean, are obtained via sample values, that which is estimated is a sample AUC and sample mean, respectively, which, in turn, are unbiased estimates of the population AUC and population mean. In the following text the terms population AUC and population mean will be used for simplicity.

[0004] There are two aspects of pharmacokinetic and pharmacokinetic/pharmacodynamic studies. From one side, the aim is to determine the exposure of a specific individual to the drug to specify the dosing regimen required in that individual to achieve a specific therapeutic goal. From the other side, the aim is to determine the average exposure to the drug that would be experienced by the individuals in a population. One method of deriving the average exposure to the drug is to determine the AUC for several individuals and then use conventional descriptive statistics to have an estimate of the population AUC. Alternatively, the measure of the population AUC can be determined directly from all the individual measurements of plasma concentrations considered together.

[0005] An important further consideration in pharmacodynamic studies is the practical limitations of the sampling protocols. In many cases the number of samples that can be collected from an individual will be limited by technical, ethical and cost reasons. For example in toxicokinetic studies only a few samples can be taken in order not to perturb the subjects under observation. Furthermore in experiments which involve destructive measurements (in which the subject is sacrificed) only one sample can be collected from each subject. An example of such a situation is the study of metabolite levels in an internal organ. In these cases the only possible option is to consider all of the individuals' samples together.

[0006] Solutions have been proposed to determine the central tendency (e.g. mean) and the dispersion (e.g. SE) of the AUC of the population directly from the individuals' samples. For example the Bailer method (“Testing for the Equality of Area Under the Curves When Using Destructive Measurement Techniques”, Journal of Pharmacokinetics and Biopharmaceutics, Vol 16, No. 3, 1988), which can be used to estimate the descriptive statistics of the AUC when destructive sampling is required and the Yeh method (“Estimation and Significant Tests of Area Under the Curve Derived from Incomplete Blood Sampling”, Proceedings of the American Statistical Association, 1990, p. 7481) which can be used to estimate the average AUC and the corresponding SE when the population is subdivided in different groups and a few samples are collected from each subject at different times.

[0007] Considering a study with the following parameters:

[0008] N subjects are utilized in the study,

[0009] samples are taken at m time points (t1, . . . , tm),

[0010] at each sampling time nj subjects are randomly selected to form m sampling groups,

[0011] such that 1$\sum _{j=1}^m{n}_j=N,$

[0012] , ie. every subject is part of one of the sampling groups, and

[0013] the measurement taken in the i-th subject at the time tj is denoted yij.

[0014] Then the Bailer method calculates the estimate of the population AUC and its standard error (SEAUC) as: 2$\begin{array}{cc}\mathrm{AUC}=\sum _{j=1}^m{a}_j{c}_j& \left(1\right)\\ {\mathrm{SE}}_{\mathrm{AUC}}=\sqrt{\sum _{j=1}^m{a}_j^2{s}_j^2/n_j}\mathrm{where}\text{:}& \left(2\right)\\ c_j=\sum _{i=1}^{n_j}\frac{y_{\mathrm{ij}}}{n_j}& \left(3\right)\\ s_j^2=\sum _{i=1}^{n_j}\frac{{\left(y_{\mathrm{ij}}-c_j\right)}^2}{\left(n_j-1\right)}& \left(4\right)\\ a_j=\{\begin{array}{cc}\left(t_2-t_1\right)/2& \mathrm{for}j=1\\ \left(t_{j+1}-t_{j-1}\right)/2& \mathrm{for}j>1\mathrm{and}j<m\\ \left(t_m-t_{m-1}\right)/2& \mathrm{for}j=m\end{array}& \left(5\right)\end{array}$

[0015] In order to find the standard error, an estimate, within each sampling group, of the average variation sj of the individual measurements yij from the average of the group cj is calculated according to (4). The standard error SEAUC is then calculated according to (2) by applying the estimate of the average variation sj over the width of each section.

[0016] One problem with the Bailer method is that in order to calculate the standard error it is assumed that the measured points of the concentration curve are independent of each other. Furthermore the Bailer method makes no distinction between measurement error and inter-individual variability.

[0017] The Yeh method is a slight generalization of the Bailer method. Consider a further study in which the parameters of the study are:

[0018] N subjects are utilized in the study,

[0019] the subjects are divided randomly into S groups (hereafter called sessions),

[0020] samples are taken at m time points,

[0021] each session contains nk subjects,

[0022] such that 3$\sum _{k=1}^S{n}_k=N$

[0023] ie. each subject is allocated to one session,

[0024] each of the nk subjects in each session is sampled more than once,

[0025] each of the nk subjects in each session is sampled at mk time points,

[0026] the sampling times of the sessions are different, and

[0027] the measurement taken in the i-th subject of the group k at time j is yijk.

[0028] Then the Yeh method calculates the estimate of the population AUC and its standard error (SEAUC) as: 4$\begin{array}{cc}\mathrm{AUC}=\sum _{k=1}^S{f}_k& \left(6\right)\\ {\mathrm{SE}}_{\mathrm{AUC}}=\sqrt{\sum _{k=1}^S{s}_k^2/n_k}\mathrm{where}\text{:}& \left(7\right)\\ f_k=\frac{1}{n_k}\sum _{i=1}^{n_k}{f}_{\mathrm{ik}}& \left(8\right)\\ f_{\mathrm{ik}}=\sum _{jϵm_k}{a}_j{y}_{\mathrm{ijk}}& \left(9\right)\\ S_k^2=\sum _{i=1}^{n_k}\frac{{\left(f_{\mathrm{ik}}-f_k\right)}^2}{n_k-1}& \left(10\right)\end{array}$

[0029] The conditions for aj are the same as for the Bailer method as shown in (5). Equation (9) shows that the first step is to calculate the estimate of the AUC for each individual in each session by using a method similar to that which was used in the Bailer method to calculate the estimate of the population AUC from the average concentrations at each of the sampling times. An estimate for the average AUC for each session is then calculated according to (8) and finally the estimate for the population AUC is calculated according to (6) by summing the contributions of each of the sessions.

[0030] The standard error is again obtained by assuming that the sampled points on the concentration curve related to different sessions are independent of each other. Also the Yeh method, as in the Bailer method, does not make a distinction between measurement error and inter-individual variability. Moreover the sampling protocol is restrictive as all the subjects in each session are sampled at the same time points and the subjects of different sessions are sampled at different time points.

[0031] In addition, with these methods it is only possible to calculate descriptive statistics (e.g., mean, SE) of the exposure to the drug in a sample of individuals. This means that it is not possible to calculate the individual values of the exposure to the drug, which, in turn, does not allow one to apply the pharmacokinetic/pharmacodynamic relationships to identify the individual dosing regimen able to realize the optimal therapeutic outcome.

[0032] It is therefore the aim of the present invention to provide an improved method of deducing the average exposure to a drug with fewer limitations on the sampling pattern and which can more accurately assess the standard error, taking account of the effects of inter-individual variability and measurement error. It is a further aim of the invention to provide a method which evaluates the individuals' exposure to a drug.

[0033] According to a first aspect of the present invention there is provided a method of determining the average exposure within a population of subjects to a pharmacological substance administered according to a given protocol by determining the area under the concentration curve, characterised in that the area under the population concentration curve is estimated by the steps of:

[0034] obtaining measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken;

[0035] building a hierarchical stochastic model composed of the population and of the individual levels;

[0036] determining the posterior probability distribution of the average AUC from the sample data and hence the average population AUC.

[0037] According to a further aspect of the present invention there is provided a method of determining the exposure of an individual and its precision within the sample of individuals given a pharmacological substance administered according to a given protocol by determining the area under the concentration-time curve, characterised in that the area under the individual concentration-time curve is estimated by the steps of:

[0038] obtaining measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken;

[0039] building a hierarchical stochastic model composed of the population and of the individual levels;

[0040] determining the posterior probability distribution of the individual AUC from the sample data.

[0041] According to a further aspect of the present invention there is provided a computer program comprised of computer code means for estimating the average total exposure within a population of subjects to a pharmacological substance, administered according to a given protocol, by estimating the area under the concentration curve (AUC); the computer code means being comprised of means to:

[0042] record measurements of the drug concentration level in each of the subjects at any time within the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken;

[0043] determine the posterior probability distribution of the population AUC from the recorded sample data using a hierarchical stochastic model of the population; and

[0044] hence determine the average population AUC.

[0045] According to a further aspect of the present invention there is provided a computer program comprised of computer code means for determining the exposure of an individual and its precision within the sample of individuals given a pharmacological substance administered according to a given protocol by determining the area under the concentration-time curve, the computer code means being comprised of means to:

[0046] record measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken;

[0047] build a hierarchical stochastic model composed of the population and of the individual levels; and

[0048] determine the posterior probability distribution of the individual AUC from the sample data.

[0049] The method of the present invention is advantageous because it can be applied to any arbitrary sampling scheme. Furthermore it allows the rigorous derivation of the posterior distribution of each parameter which means that the precision of the estimates can be made as well as the point estimates. This method is also advantageous because it provides estimates of the intra- and inter-subject variability.

[0050] The invention will now be described by way of a non-limiting example with reference to the drawing, in which:

[0051] FIG. 1 is a representation of the stochastic population model according to the invention;

[0052] The present invention is based upon a stochastic model of both the individual concentration curves and the population curve. A hierarchical model is used to describe the population model, whereas random walks describe the individual plasma concentration curves and the population curve. Model estimation is performed according to a Bayesian approach based on a Markov Chain Monte Carlo stochastic simulation algorithm. In a preferred embodiment, described below, Gibbs Sampling is used to derive the population AUC, the individuals' AUC and the other parameters describing the inter-individual variability.

[0053] In a study carried out according to the present invention there are no restrictions to the sampling protocol. The parameters of the study are therefore:

[0054] N subjects are utilized in the study,

[0055] ni measurements are taken for each subject,

[0056] the set T {t1, . . . , tm} represents all of the times when a sample is taken from any subject,

[0057] the measurement of the i-th subject at time tj is yij,

[0058] By deriving the posterior probability distribution of the population concentration curve {overscore (x)}j (j=1, . . . , m) it will be possible, by for example the trapezoidal rule, to compute the probability distribution of the population AUC. A posterior probability distribution is the probability distribution that a parameter lies in each of the various regions of the set of all its possible values after (posterior to) the observation of any values from a distribution which depends on the parameter. Thus in this case the posterior distribution will be the distribution of the probability that the concentration will be at a given level at a certain time.

[0059] In the present invention it is assumed that the curve {overscore (x)}j is, to a certain degree, regular. This is a reasonable assumption given that the signal that the invention relates to is a biological signal. A simple stochastic model of a regular signal is given by a random walk and therefore {overscore (x)}j is modelled by:

{overscore (x)} j ={overscore (x)} j −1+{overscore (ω)}j,{overscore (x)}0=0  (11)

[0060] where the variables {overscore (ω)}j are independent and normally distributed with zero mean and variance {overscore (λ)}2. The sampling times t1, . . . tm are assumed to have been chosen properly to ensure that when the signal is less regular the sampling is more frequent. Equation (11) is usefully be rewritten in the form:

P{overscore (x)} ={overscore (w)}   (12)

[0061] where:

{overscore (x)} =[{overscore (x)} 1 , . . . ,{overscore (x)} m ]T  (13)

{overscore (w)} =[{overscore (w)} 1 , . . . ,{overscore (w)} m ]T  (14)

[0062] 5 $$ P _ = [ 1 0 … 0 - 1 1 … 0 ⋮ ⋰ ⋰ ⋮ 0 ⋯ - 1 1 ] ( 15 )

{overscore (w)} ˜ .N( 0,{overscore (λ)}2I)  (16)

[0063] I is the identity matrix

[0064] A similar modelization is also made for the individual concentration curves (xij), assuming that the assumption of regularity made above is also true at the individual level. Thus for the i-th subject:

x ij =x ij −1+wij,xi0=0  (17)

[0065] where the variables wij are independent and normally distributed with mean {overscore (w)}j and variance λ2. In other words, the individual curves are random walks with expected value equal to the population curve and inter-individual variability modelled by λ2. This too is usefully rewritten as:

Px i =wi  (18)

[0066] where:

x i =[xi1, . . . ,xim]T  (19)

w i =[wi1, . . . ,wim]T  (20)

w i ˜ N ({overscore (w)},λ2I)  (21)

[0067] The measurement error is also modelled, assuming that each measurement is affected by an additive error:

y ij =x ij +v ij   (22)

[0068] where the variables vij are independent and normally distributed random variables with zero mean and variance σij2. If the measurement errors are considered as having a constant coefficient of variation (CV):

σij=CVxij  (23)

[0069] It is possible to re-write (22) as:

y i =Lixi+vi  (24)

[0070] where:

y i =[yi1, . . . ,yini]T  (25)

v i =[vi1, . . . ,vini]T  (26)

v i ˜ N (0,Σi)  (27)

[0071] 6 $$ ∑ i _ ⁢ = [ σ i1 2 0 ⋯ 0 0 σ i2 2 ⋯ 0 ⋮ ⋰ ⋰ ⋮ 0 ⋯ 0 σ ini 2 ] ( 28 )

[0072] and Li is a suitable (ni x m) matrix. Li is required because yi contains the set of all the measurements taken from the i-th individual (ie. ni values) whilst xi contains the set of all the actual values of the concentration in the i-th individual at the times T when any individual has a sample taken (ie. m values). Therefore Li is chosen such that Lixi contains the actual values of the concentration in the i-th individual at the times when a measurement is taken from the i-th individual. It is therefore the same size as yi and each value of Lixi is the actual value of the concentration in the i-th individual at the time when the measurement for the value of yi in the corresponding position was made.

[0073] By replacing xi according to (18), (24), and hence (22), is rewritten as:

y i =LiP−1wi+vi  (29)

[0074] FIG. 1 is a representation of the stochastic population model. As shown, standard prior distributions are chosen for the hyper-parameters {overscore (λ)}2 and λ2. In the case of FIG. 1 the distributions are gamma distributions. (In a gamma distribution Γ(a, b), the mean is a/b and the variance is a/b2.)

[0075] The aim of the estimation procedure is to derive the joint (and then the marginal) probability distributions of {overscore (λ)}2, λ2, {overscore (w)}, w1, . . . , wN given the data. From the joint distribution it is possible to compute probability distributions of concentration curves and AUCs (both for the population and the individuals). Point estimates and their reliability can easily be derived.

[0076] Since the computation of the joint distribution cannot be derived in closed form, a Markov Chain Monte Carlo sampling scheme is used. In this case Gibbs Sampling is chosen, and then the samples are extracted, iteratively, from the following distributions: 7$\begin{array}{c}p\left(\frac{1}{{\stackrel{\_}{\lambda }}^2}|\stackrel{\_}{\underset{\_}{w}}\right)=\Gamma \left(\frac{m}{2}+{\gamma }_1,\frac{1}{2}{\stackrel{\_}{\underset{\_}{w}}}^T\stackrel{\_}{\underset{\_}{w}}+{\gamma }_2\right)\\ p\left(\frac{1}{{\lambda }^2}|\stackrel{\_}{\underset{\_}{w}},{\underset{\_}{w}}_1,\dots ,{\underset{\_}{w}}_N\right)=\Gamma \left(\frac{N}{2}+{\gamma }_3,\frac{1}{2}\sum _{i=1}^N{\left({\underset{\_}{w}}_i-\stackrel{\_}{\underset{\_}{w}}\right)}^T\left({\underset{\_}{w}}_i-\stackrel{\_}{\underset{\_}{w}}\right)+{\gamma }_4\right)\end{array}$ p({overscore (w)}|{overscore (λ)}2,λ2,w1, . . . ,wN)=N(A−1B,A−1)

p (wi|λ2{overscore (w)},yi)=N(C−1D,C−1)

[0077] where:

A =({overscore (λ)}2+N*λ−2)I 8 $$ B _ = ∑ i = 1 N ⁢ λ - 2 ⁢ w _ _ i  C=LiTΣi−1Li+λ−2I

D =LiTΣi−1yi+λ−2{overscore (w)}

[0078] The population AUC thus determined provides a measure of the average exposure to a drug. It is further possible to determine the distribution of the individual AUCs within the population.

[0079] The method according to the present invention has been evaluated by comparing it to the Bailer and Yeh methods. The information relating to the comparison of the methods is shown in the tables below. Table 1 summarises the “experiments” which were performed. The analysed data set included plasma concentration data after a xenobiotic administration in 27 human subjects, each sampled at nine times (Rocchetti & Poggesi 1997).

[0080] In experiment 1 the methodology of the present invention was applied to the complete data set. After 1500 runs of the Gibbs Sampler (convergence is verified by using the Raftery Criterion; parameters of the prior distribution are fixed in the following way: γ1=0.6, γ2=3, γ3=0.6, γ4=10 CV is fixed to 0.1) the mean of the population AUC posterior distribution was found to be 805 and the standard error was found to be 104.

[0081] Experiments 2 and 3 were carried out to study the performance of the methodology of the present invention when the number of subjects involved in the study is reduced. In each experiment 100 randomly generated “synthetic” experiments were simulated by using a proportion (randomly selected) of the complete data set of experiment 1. A comparison of the results of experiments 2 and 3 with the results of experiment 1 is shown in table 2. It can be seen that although the mean of the population AUC posterior distribution means remains close the value obtained from the complete data set, the variance of the means increases markedly.

[0082] Experiment 4 was performed to compare methodology of the present invention with that of the Bailer method. Table 3 compares the results of the methodology of the invention, using the sampling scheme of experiment 4 (as given in Table 1), and the results of using the Bailer method, using the sampling scheme of experiment 4, with the results from experiment 1, in which the complete data set was used. The final column of the table is the difference in the results obtained from the two different methods. Again 100 “synthetic” experiments were conducted in which the samples were randomly drawn from the complete data set. It can be seen that two methods obtain almost exactly the same result if you consider the point estimate. However it can also be seen that the Bailer method significantly underestimates the Standard Error.

[0083] Experiments 5 and 6 were carried out to perform similar comparisons of results of the methodology of the present invention with the Yeh method when using the same sampling scheme. Again 100 “synthetic” experiments were performed. The results are shown in Tables 4 and 5 respectively. It can be seen from the results of this experiment that the Yeh method, as with the Bailer method, significantly underestimates the Standard Error.

[0084] Experiment 7 was performed to test the methodology of the present invention on a scheme that cannot be managed by the Bailer or Yeh methods (set out in Table 1). The scheme provides for sampling subjects in a completely “free” manner, in which each of the subjects can be sampled at any time. Again 100 “synthetic” experiments were performed by selecting random portions of the complete data set provided by experiment 1. Table 6 shows the mean and standard deviations of the mean and standard error of the population AUC posterior distribution derived from each of the “synthetic” experiments. In this case the estimates obtained are only slightly worse than the estimates obtained from the complete data set, however there is a big loss in terms of precision.

[0085] Whilst we have described above a specific embodiment of the invention, it will be appreciated that the invention may be practised otherwise than described. The description is not intended to limit the invention.

TABLE 1
Simulated Sampling Protocols
			Samples per
	Experiment	Subjects	subject	Samples	Note
	1	27	9	243	Complete
					Data Set
	2	18	9	162
	3	9	9	81
	4	27	1	27	Bailer
	5	9	3	27	Yeh 1
	6	9	3	27	Yeh 2
	7	10	—	27	Free
	Simulated sampling protocols:
	Yeh 1: subjects in session 1 are sampled at time points [t1, t4, t7], those of session 2 are sampled at [t2, t5, t8], and those of session 3 at [t3, t6, t9];
	Yeh 2: subjects in session 1 are sampled at time points [t1, t3, t5], those of session 2 are sampled at [t2, t4, t6], and those of session 3 at [t7, t8, t9];
	Free: one subject is fully sampled, one at time points [t1, t3, t5, t9], one at [t2, t3], one at [t4, t6], one at [t5, t9], one at [t6], one at [t7], one at [t8], one at [t7, t8, t9], and the last one at [t2, t1].

[0086]

TABLE 2
Population AUC Estimates Using the Method of the Present Invention
	Experiment
	1	2	3
	Mean of Population AUC Posterior	805	806	801
	Distribution Means
	Standard Deviation of Population AUC	N/A	30	71
	Posterior Distribution Means
	Mean of Standard Errors	104	127	179
	Standard Deviation of Standard Errors	N/A	15	44

[0087]

TABLE 3
Comparison of Bailer Method and the Method of the Present Invention
	Exp	Present	Bailer
	1	Invention	Method	Difference
Mean of Population AUC	805	825	825	0.057
Posterior Distribution Means
Standard Deviation of		49	48	19
Population AUC Posterior
Distribution Means
Mean of Standard Errors	104	123	70	53
Standard Deviation of Standard		31	22	34
Errors

[0088]

TABLE 4
Comparison of Yeh Method and the Method of the Present Invention
from Experiment 5
	Exp	Present	Yeh
	1	Invention	Method	Difference
Mean of Population AUC	805	786	808	−22
Posterior Distribution Means
Standard Deviation of		72	76	20
Population AUC Posterior
Distribution Means
Mean of Standard Errors	104	149	79	69
Standard Deviation of Standard		56	29	46
Errors

[0089]

TABLE 5
Comparison of Yeh Method and the Method of the Present Invention
from Experiment 6
	Exp	Present	Yeh
	1	Invention	Method	Difference
Mean of Population AUC	805	844	825	19
Posterior Distribution Means
Standard Deviation of		94	95	43
Population AUC Posterior
Distribution Means
Mean of Standard Errors	104	201	92	109
Standard Deviation of Standard		71	50	87
Errors

[0090]

TABLE 6
Population AUC Estimates Using the Method of the Present Invention
on the “Free” Sampling Scheme of Experiment 7
	Experiment
	1	2
	Mean of Population	805	796
	AUC Posterior
	Distribution Means
	Standard Deviation of		87
	Population AUC
	Posterior Distribution
	Means
	Mean of Standard Errors	104	176
	Standard Deviation of		59
	Standard Errors

Claims

1. A method of estimating the average exposure within a population of subjects to a pharmacological substance, administered according to a given protocol, by estimating the area under the concentration curve (AUC), characterised in that the area under the population concentration curve is estimated by the steps of: obtaining measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken; building a hierarchical stochastic model composed of the population and of the individual levels; determining the posterior probability distribution of the average AUC from the sample data and hence the average population AUC.

2. A method of determining the exposure of an individual and its precision within the sample of individuals given a pharmacological substance administered according to a given protocol by determining the area under the concentration-time curve, characterised in that the area under the individual concentration-time curve is estimated by the steps of: obtaining measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken; building a hierarchical stochastic model composed of the population and of the individual levels; determining the posterior probability distribution of the individual AUC from the sample data.

3. A method according to claim 1 or 2, wherein the posterior probability distribution is determined by using a Markov Chain Monte Carlo algorithm.

4. A method according to claim 1, 2 or 3, wherein the hierarchical model is obtained by: modelling the variation of the population concentration level by a random walk; modelling the individual's concentration level as a random walk with an expected value equal to the population curve; and modelling the measurement error of each sample as an additive error.

5. A method according to any one of claims 1 to 4, wherein the Markov Chain Monte Carlo algorithm used is the Gibbs Sampling method.

6. A method of determining the therapeutic and toxic effects related to the administration of a pharmacological substance comprising the steps of: administering the pharmacological substance to a population of subjects according to a given protocol; determining the average of a measure of the therapeutic or toxic effect of the pharmacological substance within the population; comparing this to the average total exposure to the pharmacological substance within the population which is determined by a method according to any one of claims 1 to 5.

7. A method of determining a safe exposure level of a subject to a pharmacological substance comprised of the steps of: administering the pharmacological substance to a plurality of populations of subjects according a plurality of different protocols; determining the average of a measure of toxic effect of the pharmacological substance within each population of subjects; comparing the average measures of the toxic effects to the average total exposure to the pharmacological substance within each population which is determined by a method according to any one of claims 1 to 5.

8. A method of determining a beneficial exposure level of a subject to a pharmacological substance comprised of the steps of: administering the pharmacological substance to a plurality of populations of subjects according a plurality of different protocols; determining the average of a measure of therapeutic effect of the pharmacological substance within each population of subjects; comparing the average measures of the therapeutic effects to the average total exposure to the pharmacological substance within each population which is determined by a method according to any one of claims 1 to 5.

9. A method of determining a dose protocol of a pharmacological substance comprising at least one of the steps of: determining the therapeutic and toxic effects related to the administration of the pharmacological substance according to claim 6;determining a safe exposure level of a subject to the pharmacological substance according to claim 7; and determining a beneficial exposure level of a subject to the pharmacological substance according to claim 8.

10. A method of treatment comprised of preparing at least one dose of a pharmacological substance according to a dose protocol determined according to claim 9; and administering said dose to the subject.

11. A computer program comprised of computer code means for estimating the average exposure within a population of subjects to a pharmacological substance, administered according to a given protocol, by estimating the area under the concentration curve (AUC), the computer code means being comprised of means to: record measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken; build a hierarchical stochastic model composed of the population and of the individual levels; determine the posterior probability distribution of the average AUC from the sample data and hence the average population AUC.

12. A computer program comprised of computer code means for determining the exposure of an individual and its precision within the sample of individuals given a pharmacological substance administered according to a given protocol by determining the area under the concentration-time curve, the computer code means being comprised of means to: record measurements of the drug concentrations in each of the subjects at any time during the study, independently from the time when measurements of the drug concentration level in the other subjects are being taken; build a hierarchical stochastic model composed of the population and of the individual levels; and determine the posterior probability distribution of the individual AUC from the sample data.

13. A computer program according to claim 11 or 12, wherein the computer code means to determine the posterior probability from the recorded sample data uses Markov Chain Monte Carlo algorithm.

14. A computer program according to claim 11, 12 or 13, wherein the computer code means obtains the hierarchical model, used to determine the posterior probability distribution from the recorded sample data, by: modelling the variation of the population concentration level by a random walk; modelling the individual's concentration level as a random walk with an expected value equal to the population curve; and modelling the measurement error of each sample as an additive error.

15. A computer program according to any one of claims 11 to 14, wherein the computer code means to determine the posterior probability from the recorded sample data uses the Gibbs Sampling method.