Introduction

Children differ from adults in their response to drugs. These differences may be caused by changes in the pharmacokinetics (PK) and/or pharmacodynamics (PD) between children and adults and may also vary among children of different ages. The PK of a drug includes the processes of absorption, distribution, metabolism and elimination of a drug, whereas the PD comprises the physiological and biological response to the administered drug and therefore may represent both efficacy and safety measures. While a child grows, enzyme pathways (involved in the PK), function and expression of receptors and proteins (involved in the PD) mature, which can be referred to as “developmental changes” in childhood. The maturation rates of these developmental changes vary, however, between the pathways and receptors and often do not correlate solely with the increase in body weight of the child. The question is therefore how to obtain data in children that allow for the study of these developmental changes, ultimately resulting in evidence-based dosing regimens for drugs in children.

To date, only a small number of drugs used in children are licensed for use in this specific group. Up to 70% of the drugs in paediatric intensive care, and 90% of the drugs in neonatal intensive care, are prescribed in an off-label or unlicensed manner [1–4]. Paediatric dosing regimens are usually empirically derived from adult regimens using linear extrapolations based on body weight. Since these developmental changes are non-linear dynamic processes, this dosing paradigm may result in under- or over-dosing, particularly in specific age groups. This may cause therapeutic failure, occurrence of severe adverse effects or even fatalities, such as fatalities occurring after long-term sedation with high doses of propofol [5, 6] and occurrence of the grey baby syndrome in neonates after treatment with chloramphenicol [7, 8]. As a result, dose adjustments in the younger age groups are often proposed. For vancomycin, for example, lower doses are administered in neonates younger than 1 week (20 mg/kg/day) compared with 1- to 4-week-old neonates (30 mg/kg/day) and children between 1 month and 18 years (40 mg/kg/day) [9].

Instead of the a priori use of body weight for dosing guidelines in children, detailed information on PK and potentially also the PD needs to be considered in order to define effective and safe dosing regimens throughout the paediatric age range. The lack of PK and PD information on drugs in children has led to the European Regulation, which came into force in 2007. This law imposes on pharmaceutical companies to perform research over the whole paediatric age range for all drugs that are developed for the European market, by requiring the submission of a paediatric investigational plan (PIP) in the early stages of the development of a new drug. In this PIP, a full description has to be given of the studies and of drug formulation in the paediatric population. If little information is available about the efficacy and safety of a drug, studies in children are only performed after more information has been obtained in the adult population to increase the safety of the paediatric study [10–12]. The main aims of introducing the paediatric regulation were to facilitate development and availability of medicines in children aged between 0 and 17 years, to improve the availability of information about medicines used in children, to ensure that the medicines are of high quality, can be administered in a safe and effective way and that paediatric studies are performed in an ethically correct way [10]. The reward for this effort is a 6-month supplementary production certificate for the pharmaceutical company.

Both for industry and for academic researchers, performing (PK–PD) studies in children in order to develop rational dosing schemes is very challenging because of ethical and practical issues. Unlike studies in healthy adults, research in healthy children is considered to be unethical, so all paediatric studies are performed in the vulnerable group of children suffering from a disease. In all clinical trials, informed consent has to be signed by the patient before he or she can be enrolled into a trial. In paediatric trials, this informed consent cannot be obtained by the patient that participates in the trial, and is therefore replaced by the consent of the parents or guardians. In older age groups, in addition to this consent, an assent is used in which the aim of the study is explained in age-appropriate language so that children can understand. [1, 4, 13].

Apart from ethical issues, practical challenges also occur when performing studies in children. There are limitations to the number and volume of samples that can be obtained, resulting in infrequent sampling possibilities and the need for advanced drug assay techniques with improved sensitivity. Another complicating factor is the limited available number of subjects that suffer from the same disease. Finally, pharmacodynamic endpoints that measure the efficacy of the drug and that are validated for children may be lacking. All these factors call for highly advanced study designs and analysis techniques so that the burden for each child can be kept to a minimum while still addressing all the study objectives.

This article aims to inform clinical pharmacologists, paediatricians and pharmacists about population PK–PD modelling in paediatric drug research. Advanced statistical tools are discussed that can be used to develop rational dosing schemes based on the PK and PD of a drug in children, despite practical and ethical restrictions. Using these tools, covariates can be identified in order to define appropriate doses and dosing intervals based on the individual characteristics of each child with minimal burden to each patient. The article also describes how to evaluate the predictive performance of the models by different validation methods including a prospective clinical trial. Ultimately, the efforts result in an individualised dosing regimen based on the PK–PD relation throughout the paediatric age range.

PK–PD in children

Developmental changes in childhood can affect all PK processes from absorption until elimination as well as the pharmacodynamic effects. For example, in neonates intra-gastric pH is elevated (>4), which may increase the bioavailability of acid-labile compounds (penicillin G) and decrease the bioavailability of weak acids (phenobarbital) when given orally [14]. Additionally, gastric emptying in neonates is delayed, which means that also the absorption of drugs, e.g. paracetamol is slower in neonates [15, 16]. Other examples are changes in the metabolising enzyme capacity in children. Although most uridine 5′-diphosphate (UDP)-glucuronosyltransferases (UGTs) and P-450 cytochromes (CYPs) are expressed during the first week of life, the activity at birth in comparison with adults is often low, e.g. UGT2B7 activity at birth is around 10% of the adult level and different enzyme systems are known to mature at different rates [14, 17–20].

In addition, renal function and liver flow are influenced by physiological changes depending on age, e.g. the glomerular filtration rate (GFR) in mL/min/70 kg in full-term neonates is 35% of the adult value, while mL/min/70 kg adult values are reached at approximately 1 year old [21]). When using units of mL/min/70 kg, however, it should be realised that actual values of GFR in children are still very low compared with adult values because of the correction for differences in total body weight between adults and infants.

Furthermore, the body composition of children changes continuously, resulting in an age-dependent proportion of body water and fat, which influences the distribution of drugs. For example, the total amount of body water (80–90% of the body weight) is higher in neonates than in adults (55–60%). Hydrophilic drugs like aminoglycosides have a larger volume of distribution in neonates, which can be explained by the larger amount of extra-cellular fluid (45% of the body weight) compared with adults (20%) [14, 22].

In order to characterise the specific influence of developmental changes in childhood on the PK of a drug, concentration–time profiles are necessary, which require measurements of drug concentrations. For ethical reasons, in paediatric studies, discomfort, like pain and anxiety associated with venipuncture, must be restricted and practical issues limit the volume and amount of blood samples that can be obtained. Therefore, sensitive analysis techniques requiring only small blood samples should be used. While HPLC methods have been reported to require only 50 μL of blood [23], more recently LC-MS methods can measure up to ten different drugs in volumes as low as 50–100 μL [24]. Additionally, also alternative matrices such as saliva should be explored as a non-invasive, more child-friendly alternative to measure a drug concentration. An example in this respect is an LC-MS/MS method that was developed and validated for the measurement of busulphan in saliva [24]. Also, the use of a dried blood spot method, e.g. for tacrolimus, can facilitate the measurement of drugs in children [25]. Another method is capillary electrophoresis, which requires only a low sample volume for the quantification of drugs in biological fluids [26].

Changes between children and adults may also result from differences in the pharmacodynamics of a drug in children, e.g. by changes in the relative number and function of receptors. These age-related PD differences are until present rarely reported in the literature, but one of the few examples is the increased sensitivity to d-tubocurarine, an antagonist of nicotinic neuromuscular acetylcholine receptors, in neonates and infants compared with children and adults [27]. Other examples are the observed lower minimum alveolar concentration (MAC) of isoflurane in pre-term neonates compared with full-term neonates and older children [28, 29] and the different sensitivity to bronchodilators because of the lack of smooth muscles in the airways in neonates [30] .

To study the PD of a drug in children, the use of a PD endpoint that has been validated for use in children is a prerequisite. An illustrative example is the measurement of pain in young children. Since they are not able to report their pain using a visual analogue scale, an observational scale has been developed. This comfort behavioural (COMFORT-B) scale was developed and validated for use in children under the age of 3 years [31]. The scale assesses six behavioural items: alertness, calmness, muscle tone, body movement, facial tension and crying (non-ventilated children) or respiratory response (ventilated children). All items range from 1 (no distress) to 5 (severe distress), resulting in a total score varying from 6 to 30. This validated scale can then be used as a PD endpoint for the development of PD models for pain and/or sedation in children of different ages [32–34].

The influence of covariates such as the developmental changes, disease status and genetics on the PK and PD of drugs in children is depicted in Fig. 1.

Fig. 1
figure 1

Schematic representation of the relationship between dose and concentration (pharmacokinetics, PK) and between concentration and a pharmacological (side) effect (pharmacodynamics, PD). Important covariates that may affect both the PK and/or PD are body weight, age, disease status (e.g. critically ill versus healthy children) and genetics

Full size image

When both the PK and PD of a drug in children are characterised, the models developed can be used to derive rational dosing regimens with predictable efficacy and concentration profiles. An example of such a PK–PD model with a derived dosing regimen is an article published by Peeters et al. In this article both the PK and the PD were characterised in children, the latter with the use of the COMFORT-B scale as a pharmacodynamic endpoint [33]. Based on the model it was found that propofol clearance is two times higher in non-ventilated children than in ventilated children and adults. For the PD, a model was derived in which an effect of propofol was characterised within a naturally occurring sleep pattern of children in the ICU. Both models (PK as well as PD) were used to simulate concentrations as well as the effects that could be expected using different dosing schemes (Fig. 2). As a result, based on this PK–PD model, a propofol dose of 30 mg/h was recommended for a child of 10 kg, which will result in adequate COMFORT-B scales in the night following craniofacial surgery.

Fig. 2
figure 2

Simulation of propofol concentrations and response using COMFORT-B score versus time based on developed PK and PD models, after administration of different doses of propofol (0, 18, 30, and 36 mg/h) for a 10-kg and a 5-kg non-ventilated infant in the first night at the Intensive Care Unit following craniofacial surgery. Target COMFORT-B scores are between 12 and 14 preferably. Reproduced from [Peeters MY, Prins SA, Knibbe CA, DeJongh J, van Schaik RH, van Dijk M, et al. Propofol pharmacokinetics and pharmacodynamics for depth of sedation in nonventilated infants after major craniofacial surgery. Anesthesiology 2006 Mar;104(3):466-74.]

Full size image

Methods of analysing data: standard two-stage or population approach

When concentration–time and concentration–effect datasets obtained in children are considered for analysis, two different methods can be applied: the standard two-stage approach and the population approach using non-linear mixed effect models [35–38]. When using the standard two-stage approach or classical approach, in a first step parameters are estimated in each individual based on individual concentration–time profiles (Fig. 3a). In a second step, these parameters are summarised by calculating the mean or median of the parameters and the variability between subjects (SE or IQR). A major drawback of this methodology is that it requires a relatively high number of samples in each individual patient (Fig. 3a), while each patient has to contribute roughly the same number of samples. Moreover, it is very difficult to distinguish between inter-individual (variability between subjects) and intra-individual or residual variability (variability within one subject, measurement error and model misspecification) and as a result inter-individual variability is often overestimated [39].

Fig. 3
figure 3

Concentration–time profiles of the same study using two different approaches. In a the standard two-stage approach is applied to a rich dataset. b shows the population approach with mixed effect modelling applied to the same dataset using only two data points for each individual so that a sparse dataset is created. In a, in each of the six individuals 10 samples are available. The different symbols correspond to different individuals. Each black line corresponds to a separate fit to the 10 data points of each individual. In b, which uses the mixed effect modelling approach, two samples of the 10 per subject in a are used. The different symbols correspond to the six different individuals. The black line illustrates the concentration–time plot based on the population mean values of the PK parameters (PRED). The grey lines show the plots of the individual patients, which are based on the population mean values together with the measured concentrations of the specific individual (IPRED)

Full size image

Since usually only a limited number of observations can be obtained in paediatric subjects, the population approach using non-linear mixed effect modelling to obtain PK and PD parameters is the preferred approach [37]. The population approach differs from the standard two-stage approach in the fact that the analysis is based on simultaneous analysis of all data of the entire population, while still taking into account that different observations come from different patients (Fig. 3b). Additionally, the population approach allows not only for the analysis of dense data, but also for sparse (limited number of observations per individual) and unbalanced data (unequal distribution of observations in various parts of the concentration–time profile in the individuals) or a combination of both. Finally, both the inter-individual and intra-individual variability are separately estimated in the dataset using this approach.

As a result of this methodology, when designing a paediatric study of which the data will be analysed using the population approach, it is advisable to collect samples at different times (or time windows) or to set alternating sampling schemes in subgroups of patients. This also means that (part of the) samples can be collected during routine clinical sampling. Consequently, the burden for the child who participates in the trial is reduced and the statistical power to develop a model describing the concentration–time or concentration–effect profile is not affected or improved.

The term “mixed” in non-linear mixed effects modelling represents a mixture of fixed and random effects. For the fixed effects, a structural model describing the PK or PD is chosen (e.g. a two-compartment model for PK or an Emax model for PD). The random effects quantify the variability that is not explained by the fixed effects. These random effects include inter-subject and intra-subject random variability (Fig. 4), which are both simultaneously and separately estimated. It is often assumed that the variability between subjects follows a normal distribution with a mean of zero and variance ω2. Equation 1 is used to describe the relationship between individual and population parameter estimates.

$$ \mathop \theta \nolimits_i = \mathop \theta \nolimits_{mean} \bullet \mathop e\nolimits^{\eta I} $$
(1)

where θ i represents the parameter of the ith subject, θ mean the population mean, and ηi the variability between subjects. The residual error is generally described using a proportional error (error is dependent on the concentration, which means a higher absolute error at higher concentrations (Eq. 2) or an additive error (constant for all observations (Eq. 3) or a combination of both. This means for the jth observed concentration of the Ith individual the relation (Yij):

$$ \mathop Y\nolimits_{ij} = \mathop c\nolimits_{pred,ij} \bullet \left( {1 + \mathop \varepsilon \nolimits_{ij} } \right) $$
(2)
$$ \mathop Y\nolimits_{ij} = \mathop c\nolimits_{pred,ij} + \mathop \varepsilon \nolimits_{ij} $$
(3)

where cpred is predicted concentration and εij is a random variable with a mean of zero and a variance of σ 2.

Fig. 4
figure 4

In a, the inter-individual variability among three individuals who received the same dose is shown. b presents the intra-individual or residual variability by showing the concentration–time profile after repeated administration. Both these random variables are assumed to be normally distributed with a mean of zero and a variance of ω2 or σ2 respectively

Full size image

In general, model building requires three different steps. First, a structural model (fixed effects) has to be designed, then a statistical sub-model (random effects) has to be developed and in the final step a covariate sub-model is identified.

The structural model describes the overall trend in the data. The choice of structural model (e.g. one-, two- or three-compartment model for PK and an Emax model for PD) is to be based upon the best a priori information about the drug to be studied [40]. The structural model uses fixed effects parameters such as clearance and volume of distribution for PK or Emax and EC50 for PD. The population values for these parameters are called typical values (TV).

After selecting the structural model, the statistical sub-model, which accounts for the inter-individual as well as the residual variability, is chosen and tested. Information on the inter- and intra-individual or residual variability is of clinical value, because it describes differences in clinical response between and within patients and may therefore provide guidance for rational dose adjustments. With the population approach, both these random effects are obtained, apart from estimates of both the population values (TV) and the individual values of PK and PD parameters (so-called post hoc parameter estimates).

In the final step the covariate sub-model is determined, which expresses relationships between covariates and parameters of the structural model (e.g. influence of body weight on volume of distribution or clearance). Covariates can be individual-specific (age, body weight, genetic profile, etc.) or time-varying (renal function, haemodynamic parameters, body temperature, etc.). The covariate analysis will be explained in more detail in the following section.

As these three models are inter-related, the choice of the structural (and statistical) model may affect the choice of the covariate model and vice versa. The process of finding a model that adequately describes the data is thus an elaborate task, where model checking/refining is performed in several steps. To assess model fit in relation to the observed concentrations or effect measures, scatter plots or the so-called goodness-of-fit plots are created (see Validation of the PK–PD models). Free software packages (Xpose, PSN, etc.) are available to generate these plots.

The most commonly used software package for model building, which is also supported by the European Medicines Agency (EMEA) is the non-linear mixed-effect modelling program NONMEM (GloboMax/ ICON, Ellicott City, MD, USA) [4, 41–43]. NONMEM estimates parameters (e.g. clearance, volume of distribution or EC50) via a maximum likelihood approach. This means that with the given data, the estimations of the parameters are those that occur with the highest probability. Alternative software packages that can be used are, for example, Monolix, WinNonMix, USC*PAC, which uses non-parametric maximum likelihood methods [44], or ADAPT, using maximum a posteriori (MAP) methods [45].

Covariate analysis

To determine the optimal dose based on the individual characteristics of a patient, a covariate analysis has to be performed [40, 46, 47]. The aim of the covariate analysis is to identify specific predictors (covariates) of PK and PD variability and can typically be studied in population models. Covariate analysis involves the modelling of the distribution of the individual parameter estimates as a function of covariates, which can be of demographic (e.g. age, body weight, gender), patho-physiological (e.g. renal or hepatic function) and genetic/environmental origin, and/or be the result of the concomitant use of other drugs, which may influence the PK and/or PD. The identification of predictive covariates for variability provides the scientific basis for rational and individualised, patient-tailored dosing schemes.

The influence of developmental changes in childhood can be explored primarily by using size and/or age as covariates. Size (body weight) can be incorporated into the model using two different approaches. The first approach, the “allometric size approach”, includes size a priori by using a body weight-based exponential equation with a fixed exponent of 0.75 for clearance and 1 for volume of distribution [48–52]. Once size is incorporated in the model using this fixed method, the influence of age is investigated, being the difference between the actual value of the PK parameter and the 0.75 allometric equation. When incorporating age as a covariate, different age descriptors may be used like postmenstrual age (PMA), gestational age (GA) or postnatal age (PNA) [53]. The choice of any of these age descriptors is based on the results of the systematic covariate analysis as described below [50, 54]. In the second approach, the “systematic covariate analysis”, body weight is regarded as a covariate as any other, which means that the descriptive properties on the PK parameters are evaluated in a systematic covariate analysis as described below [55–57].

In a systematic analysis, when studying the influence of covariates, scatter plots and summary plots of individual parameter estimates and/or weighted residuals versus covariates are used to screen for appropriate covariates to include in the covariate sub-model. Additionally, these plots are used to explore the nature of the influence of the covariate (linear, exponential, allometric, subpopulations, etc.). Likely candidate covariates are then added to the model (forward inclusion). The influence of each covariate on the parameters is examined separately and compared with the simple model (no covariates). To assess whether the model with the covariate statistically improved the fit to the data, the difference between their objective function value, referred to as the log–likelihood ratio, is calculated. This ratio is assumed to be Chi-squared distributed, which means that a reduction in objective function of 3.84 is considered to be significant (P < 0.05) [43, 58]. In addition, the reduction in objective function, goodness-of-fit plots of the simple model and covariate model are explored for diagnostic purposes. Furthermore, the confidence interval of the parameter estimates, the correlation matrix (indicates the relationship between two structural parameters) and visual improvement of the individual plots are used to evaluate the model. Finally, a superior model is expected to reduce the inter-subject variance and/or the residual error terms. This procedure of covariate modelling implies that each covariate is only implemented if it can be fully justified by the data and the results of the statistic evaluations.

When two or more covariates are found to significantly improve the model, the covariate that reduced the objective function most is included in the model after which the other covariates are tested again for their significance. After all covariates that significantly improved the objective function are added to the simple model, a backward deletion is performed, which means that each covariate is removed from the full model, one at a time (the one that causes the smallest increase in objective function first). Retaining or removing the covariate is statistically tested by the use of the objective function (Chi-squared test) until each covariate has been tested.

In datasets containing sparse data, there may not be enough information to accurately estimate inter- and intra-individual variability. This causes the values of these parameters to shrink to 0, resulting in individual parameter estimates that are closer to the population parameter estimates than they really are. This phenomenon is called shrinkage [59]. Shrinkage may cause individual predictions, individual parameter estimates and diagnostics based on them to be less reliable. It can also hide, falsely introduce or distort the shape of covariate relationships.

Shrinkage is the result of the properties of the data and is therefore difficult to avoid. One can only be aware of the presence of shrinkage, realise the influence it may have on the covariate analysis and use diagnostics other than those based on individual predictions or individual weighted residuals in the model building and model evaluation procedures.

Validation of PK–PD models

The objective of a PK or PK–PD modelling exercise is usually not just to describe the dataset of the sample of individuals that were studied. Generally, models are used to simulate which concentrations and/or effects and their variability can be expected when different doses are given to future patients. These simulations may therefore lead to optimised dosing recommendations or to optimisation of new studies for the entire population to which the sample of individuals belongs. It is often said that “all models are wrong, but some are useful” [60]. In order to define whether a model is useful and valid for clinical and trial simulations, thorough evaluation and validation of the model is necessary. Although validations of PK models are only performed in 17% of the published paediatric studies [4] and in 28% of the adult studies [61], proper model validations are an essential step in model building. For this purpose, different evaluation and validation methods are available. As described before [62], a proper validation and evaluation procedure includes an internal model evaluation followed by an external evaluation and a prospective clinical study.

The first evaluation method is the basic internal model validation used to assess whether the model is able to describe the learning dataset (the dataset used to develop the model) accurately and without bias. This evaluation should actually be considered the final stage of the model building procedure. Subsequently, in the external evaluation, it is assessed whether the model is able to describe one or more external datasets (the datasets other than the one used to develop the model) adequately. Alternatively, if a dataset is sufficiently large the original dataset may be split in two so that the model is developed using one part (about two thirds) of the dataset and evaluated externally using the other part (one third) of the dataset. In paediatric studies, it is then especially important to stratify the data correctly and ascertain that all age groups are represented in equal proportions in both datasets.

Various techniques are available for the validation and evaluation of population PK and PK–PD models (both for internal and external validation procedures).

Basic goodness-of-fit plots

  1. 1.

    Individual predicted vs observed concentrations

  2. 2.

    Population predicted vs observed

  3. 3.

    (Conditional) weighted residual vs time

  4. 4.

    (Conditional) weighted residuals vs dependent variable

Weighted residuals (WRES) and conditional weighted residuals (CWRES) are calculated as follows:

$$ WRES = \frac{{{{\vec y}_I} - {E_{FO}}\left( {{{\vec y}_I}} \right)}}{{\sqrt {{Co{v_{FO}}\left( {{{\vec y}_I}} \right)}} }} $$
(4)
$$ CWRES = \frac{{{{\vec y}_I} - {E_{FOCE}}\left( {{{\vec y}_I}} \right)}}{{\sqrt {{Co{v_{FOCE}}\left( {{{\vec y}_I}} \right)}} }} $$
(5)

where \( {\vec y_i} \) is the vector of the measurements, \( E\left( {{{\vec y}_I}} \right) \) is the expectation of the data and \( Cov\left( {{{\vec y}_I}} \right) \) is the covariance matrix of the data [63].

These plots are used in model building, but can also be used to ascertain that there is no trend or bias in the model predictions of the final model. Furthermore, these plots can be used for both the internal and external evaluation of the model.

Bootstrap analysis

In a bootstrap analysis new datasets are generated by resampling from the original dataset and it is therefore an internal validation of the model. The new datasets are subsequently refitted to the original model, yielding mean values and standard errors for every model parameter.

A bootstrap analysis provides information on the stability of the model and its dependence on specific individuals in the learning dataset. With the freely available PSN or Wings for NONMEM software packages an automated bootstrap analysis can be performed.

Visual predictive check

In a visual predictive check (VPC) [64] a PK or PD profile is simulated a 100 to 1,000 times and lines for the median values and their 90% prediction interval are plotted in a graph. The observed values in the internal or external dataset are subsequently plotted on top of this. It can then be visually checked whether 90% of the observations are within the indicated prediction interval and whether there is no bias in the observations compared with the prediction interval. In Fig. 5, two examples of a VPC are given, showing when a model does not work and when a model does work on the same data.

Fig. 5
figure 5

Two examples of a visual predictive check (VPC) are illustrated based on the same dataset (warfarin concentrations and prothrombin complex activity [PCA]) using two different models. In a the VPC of the effect compartment is shown, while in b the VPC of the turnover model is demonstrated. Both real (grey) and simulated (black) data are displayed with mean (thick lines) and 95% intervals (thin lines). Based on both graphics, the turnover model is the most appropriate model since real and simulated lines are in good agreement. Reproduced from [Karlsson, M. O. and N. H. G. Holford (2008). “A Tutorial on Visual Predictive Checks.” PAGE 17 (2008) Abstr 1434 (www.page-meeting.org/?abstract=1434)]

Full size image

The VPC is a simulation-based diagnostic that can be used when the PK or PD profiles for all individuals in the dataset are similar and it allows for easy interpretation of the result. For this diagnostic tool, there are no statistical tests and all evaluations are based on visual assessment. When the individual profiles are expected to deviate largely from one another because there is for instance a large variability in the time and amount of dose administrated, or when there are many covariates, the use of this diagnostic becomes more difficult.

Normalised prediction distribution error

Another simulation-based diagnostic that can be used for both internal and external validations is the normalised prediction distribution error (NPDE) [65]. An example of an NPDE previously published is shown in Fig. 6 [55]. This method yields information on how accurate the model predicts the median value of the observations and the variability within them. The interpretation of this diagnostic is less straightforward than for the VPC, but the advantage of this method is that it can be used when the variability in dosing regimen (both in time, amounts and rates) is high or when there is a large number of covariates in the model. This can for instance be the case for data obtained during routine paediatric clinical practice. Software (e.g. NPDE add-on package for R) [66] to perform this analysis is freely available. For the NPDE, besides visual evaluation of the plots, statistical tests are available. These statistical tests are, however, reported to be highly sensitive and powerful, so that decisions for the model should primarily be based on visual assessments. An example is the statistically significant deviation of zero of the mean value because of the large number of data, while the actual deviation is small (e.g. 0.074) and not of clinical relevance.

Fig. 6
figure 6

Example of a normalised prediction distribution error (NPDE) analysis, which shows the NPDE distributions for morphine. The normal distribution is represented by the solid line. The values for the mean and standard deviation of the observed NPDE distribution are given below the histogram, with an asterisk indicating a significant difference of a mean of 0 and a variance of 1 at the P < 0.05 level, as determined by the Wilcoxon signed rank test and Fisher’s test for variance. The distribution of NPDE vs time after the first dose and NPDE vs the log of the concentrations are also shown. The dotted lines represent the 90% distribution of the NPDE. Reproduced from [Knibbe CA, Krekels EH, van den Anker JN, DeJongh J, Santen GW, van Dijk M, et al. Morphine glucuronidation in preterm neonates, infants and children younger than 3 years. Clin Pharmacokinet 2009;48(6):371-85] with permission from Wolters Kluwer Health | Adis (© Adis Data Information BV [2006]). All rights reserved

Full size image

If the model performed well in both evaluation procedures, the dosing algorithm that results from the PK–PD model needs to be tested and challenged in a prospective (clinical) trial. If the predictive performance of the model is corroborated by the trial it can be used with confidence in clinical practice.

Optimal design of paediatric studies

When new population PK–PD studies are performed, it is important to design these studies in the most efficient manner possible to obtain maximum information about the PK and PD parameters so that they can be determined with the highest precision [51, 67]. When designing PK–PD studies in paediatrics certain factors need to be taken into account, e.g. the age range of the paediatric group, the therapeutic index, the possibility of collecting blood samples, the availability of validated PD endpoints for children and the availability of sensitive analytical methods.

When optimising a PK or PK–PD study design, using literature data from adults or children of different age ranges or possibly in vitro or pre-clinical data, a concentration–time or effect–time profile for a study can be simulated. This can help to identify possible shortcomings in the design or to perform a power analysis. Alternatively, software packages are available (WINPOPT [68], PopED [67] and PFIM [69]) that can help to identify the optimal number and time points of observations in a study based on the prior information on a drug [70]. To determine the appropriate sample size certain factors, which are summarised in Table 1, need to be taken into account. Each of these factors can influence the required number of patients and/or samples in a positive or negative way. In a study by Peeters et al. [32] only 24 patients (aged between 3 and 24 months) were required to determine both the PK and PD, since rich sampling was performed (a median of 11 samples per child) and no covariates were found in the relatively homogeneous population. This is in contrast to a study performed by Knibbe et al. [55], in which 250 children were included. This higher number was required because in addition to the large dispersion in age from (pre-term) neonates up to toddlers of 3 years of age, only 1 to 4 samples were available for each subject. Moreover, infusion rates and additional bolus doses varied for each child during the study to obtain the desired analgesic effect. In another example [71], only 6 patients (aged between 1 and 5 years) were required, in which 7 samples per patient were collected. This lower number of patients (n = 6) compared with the study by Peeters et al. (n = 24) can be explained because there often exists a lower variability in PK than in PD, which results in a lower number of patients required (Table 1).

Table 1 Factors influencing the required number of patients and/or samples per patient
Full size table

Conclusions and perspectives

In view of the European Regulation, which came into force in 2007, it now seems time to use the progress that has been made in the field of integrated PK–PD modelling [72] to develop rational and individualised dosing schemes for children. Because of the possibility of analysing sparse and unbalanced datasets, thereby minimising the burden for each child, population PK–PD modelling and simulation using non-linear mixed effect modelling has become the preferred tool to develop effective and safe dosing regimens for children. Specifically in paediatrics, where developmental changes have to be taken into account, which may influence the PK and/or the PD of the drugs, this advanced statistical tool is of critical value.

Before dosing regimens can be tested in clinical practice, proper validations of the models should be performed, for which adequate tools have recently been developed. Besides internal and external validation, prospective clinical trials, which allow for the evaluation of the model-based dosing regimens, are needed, not only to adjust the proposed dosing regimen, but also to convince paediatricians to use the information that has been generated using these modelling exercises.

Furthermore, one of the future goals may be to explore the possibilities of cross-validation of the models, in which the reported influences of developmental changes on a certain PK or PD parameter of one drug are evaluated for use in another drug that goes through the same metabolic route or shares the same mechanism of action. In this respect, physiologically based pharmacokinetic (PBPK) models are needed. PBPK models consider the physiological and biochemical processes by using in vitro data to describe the PK of drugs [73, 74]. The combination of these two approaches may use the information that is already available in an optimal way in defining effective and safe dosing regimens for each individual patient.

In conclusion, analyses of paediatric data using population PK–PD modelling and covariate analysis will result in individualised dosing regimens for children of different ages, body weights and genetic backgrounds. Thus, population PK–PD modelling constitutes an innovative approach to the study of drug effects in this very special patient population, which is otherwise difficult to study.