PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Eur J Pharm Sci. Author manuscript; available in PMC 2017 September 20.
Published in final edited form as:
PMCID: PMC4984529
NIHMSID: NIHMS803234

Pharmacokinetic Models of Morphine and its Metabolites in Neonates

Systematic Comparisons of Models from the Literature, and Development of a New Meta-Model

Abstract

Morphine is commonly used for pain management in preterm neonates. The aims of this study were to compare published models of neonatal pharmacokinetics of morphine and its metabolites with a new dataset, and to combine the characteristics of the best predictive models to design a meta-model for morphine and its metabolites in preterm neonates. Moreover, the concentration-analgesia relationship for morphine in this clinical setting was also investigated.

A population of 30 preterm neonates (gestational age: 23–32 weeks) received a loading dose of morphine (50–100 μg/kg), followed by a continuous infusion (5–10 μg/kg/h) until analgesia was no longer required. Pain was assessed using the Premature Infant Pain Profile. Five published population models were compared using numerical and graphical tests of goodness-of-fit and predictive performance. Population modelling was conducted using NONMEM® and the $PRIOR subroutine to describe the time-course of plasma concentrations of morphine, morphine-3-glucuronide, and morphine-6-glucuronide, and the concentration-analgesia relationship for morphine.

No published model adequately described morphine concentrations in this new dataset. Previously published population pharmacokinetic models of morphine, morphine-3-glucuronide, and morphine-6-glucuronide were combined into a meta-model. The meta-model provided an adequate description of the time-course of morphine and the concentrations of its metabolites in preterm neonates. Allometric weight scaling was applied to all clearance and volume terms. Maturation of morphine clearance was described as a function of postmenstrual age, while maturation of metabolite elimination was described as a function of postnatal age. A clear relationship between morphine concentrations and pain score was not established.

Keywords: Morphine, Morphine-6-glucuronide, Preterm Neonates, Pharmacokinetics/pharmacodynamics, Modelling

Graphical abstract

An external file that holds a picture, illustration, etc.
Object name is nihms803234u1.jpg

1 INTRODUCTION

Preterm neonates of low birth weight experience significant pain and stress during their stay in the neonatal intensive care unit (Canadian Paediatric Society, 2000; Johnston and Stevens, 1996; Simons et al., 2003). The physiological instability and underlying diseases of these infants necessitates advanced perinatal care. Different factors such as mechanical ventilation, invasive procedures, postoperative issues, and acute medical illness due to prematurity (Bhalla et al., 2014) may be responsible for pain in this vulnerable population. Morphine is the main analgesic used for pharmacological pain relief in preterm neonates (Chay et al., 1992). However, the response to morphine is highly variable and difficult to predict, so optimal analgesia remains a challenge (Carbajal et al., 2005; Tibboel et al., 2005).

Neonatal pharmacokinetics is affected by developmental aspects of morphine metabolism and excretion, particularly maturation of organ function, and variability in body size (Alcorn and McNamara, 2002; de Wildt et al., 1999; van den Anker, 1996). The major metabolic pathway of morphine is the formation of morphine-3-glucuronide (M3G) and morphine-6-glucuronide (M6G). Changes in enzyme systems and decreased hepatic function causes a reduced glucuronidation capability in preterm neonates (Choonara et al., 1989). Reduced glomerular filtration rate and impaired renal function also affects morphine pharmacokinetics and especially the excretion of metabolites (Choonara et al., 1989). These rapid changes in body composition and renal function results in extensive inter-individual and intra-individual variability in morphine exposure. This makes the prediction of morphine concentrations after a given dose a challenge, and therefore makes it difficult to achieve a desired target concentration.

The effects observed after morphine administration are primarily due to μ-opioid receptor activation. The active morphine metabolite, M6G, is also a μ-opioid receptor agonist and is believed to contribute to analgesia (Kilpatrick and Smith, 2005). No significant signs of pharmacological activity have been shown for M3G in humans (Penson et al., 2001). The concentration-response relationship for morphine used in acute pain in preterm and term neonates is however still not characterised (Anand et al., 2008; Anderson and van den Anker, 2014; Tibboel et al., 2005).

In recent years, several studies have examined variability and maturational effects to increase our understanding of morphine pharmacokinetics in paediatrics (Anand et al., 2008; Bouwmeester et al., 2004; Knibbe et al., 2009; Wang et al., 2013). However, as far as we know there are currently no models that capture the antenatal maturation of morphine clearance separated from size-related factors, and describe metabolite pharmacokinetics in preterm neonates at the same time. Based on a new dataset from preterm neonates receiving morphine on clinical indication, collected in two neonatal intensive care units (NICU) over 5 years, the aims of the current study were to compare the observed concentrations of morphine and its metabolites with the predictions based on the currently existing published models, to investigate if the new data would result in an update of currently described parameters, and to combine the characteristics of the best models to propose a meta-model that is able to describe the new data. Moreover an attempt to describe the analgesia-concentration relationship for morphine was planned.

2 METHODS

2.1 Study design

The current study used data from a clinical trial conducted at the Children’s National Health System, Washington DC, and the Kosair Children’s Hospital, Louisville, KY, USA. Thirty preterm neonates were enrolled in the study from December 2005 to June 2009. The original proposal was designed with the goal to enroll 60 neonates but despite multiple efforts including adding other clinical site, only 30 neonates were enrolled. Subjects received morphine as part of their routine medical management of pain. The inclusion criteria were neonates with a gestational age (GA) at birth from 22 to less than 32 weeks, a postnatal age (PNA) of less than 30 days, an indwelling arterial line already in place for clinical purposes, and a clinical indication for receiving intravenous (IV) morphine with an anticipated therapy for at least 24 hours. Neonates who received morphine prior to the study were eligible for enrolment. The exclusion criteria were neonates with severe asphyxia, serious intraventricular haemorrhage, major congenital/facial malformations, neurological disorders, and neonates with clinical or biochemical evidence of hepatic and/or renal failure. Neonates who were receiving continuous or intermittent neuromuscular blockers and had received drugs that were UGT2B7 substrates were also excluded.

2.1.1 Study medication

The study drug was morphine sulphate. Neonates with a gestational age of less than 29 weeks received a 30-minute IV loading infusion of 50 μg/kg at time of study enrolment. This loading dose was immediately followed by a continuous IV infusion of 5 μg/kg/h as a maintenance dose. Neonates with a gestational age of 29 weeks or more received a 100 μg/kg loading infusion and a continuous infusion of 10 μg/kg/h as a maintenance dose. Decisions to administer additional breakthrough doses of morphine (30-minute IV bolus dose of 50 μg/kg) or to discontinue the morphine infusion were made at the neonatologist’s discretion. The total duration of the study for each infant was dependent on the duration of the continuous infusion.

2.1.2 Blood samples

Sampling started at time of enrolment, and 200 μL blood was drawn at each predetermined blood sampling time point. The first PK sample was drawn just before the administration of a loading dose, and subsequent PK samples were scheduled at 1, 4, 8, 12, 24, 48, 72, and 96 hours after study enrolment. One additional sample was collected 24 hours after morphine discontinuation. If doses for breakthrough pain were given, an additional PK sample was collected just before the additional morphine administration. Morphine, M3G, and M6G were quantitated in plasma, using a validated HPLC-MS/MS method (Meng et al., 2000). The analysis of the samples was carried out at the Paediatric Pharmacology Research Unit of the University of Utah.

2.1.3 Pain assessment

At each of the scheduled blood sampling time points (t = baseline, 1, 4, 8, 12, 24, 48, 72, and 96 hours), each subject was videotaped by a whole-body camera and a separate camera focused on their face. Neonatologists scored each patient’s pain using the videotapes and the Premature Infant Pain Profile (PIPP; Stevens et al., 1996). Each infant had continuous monitoring of vital signs, such as temperature, respiratory rate, blood pressure, heart rate, and oxygen saturation (SaO2). The PK/PD analysis was performed on the PIPP score, and other components of the study will not be discussed.

2.2 Pharmacokinetic-pharmacodynamic modelling

An overview of the full modelling process can be found in Figure 1. All steps in the modelling analysis are explained in detail in the following sections.

Figure 1
Flowchart for the entire modelling analysis. The top part (red and blue) shows the step for the literature model comparison (as described in section 2.2.2). The green part shows the step for development of a pharmacokinetic meta model (as described in ...

2.2.1 Software

The PK/PD analysis was performed using nonlinear mixed effects modelling. The population modelling process, simulations, and bootstraps were performed with NONMEM version VII, level 3.0 (Beal et al., 2009) using the Wings for NONMEM interface (http://wfn.sourceforge.net/) and Intel Fortran compiler. The first-order conditional estimation with interaction (FOCE-I) method in NONMEM was used for estimation of the population parameters. Individual parameter estimates were obtained using the Bayesian POSTHOC functionality of NONMEM. The ADVAN 13 subroutine in NONMEM was used for the PK/PD modelling. Raw data manipulation, processing of NONMEM output, visual representation of the data, and graphical outputs were performed using the R data analysis language, version 3.1.1 (R Development Core Team, 2014), with the ggplot2 (Wickham, 2009), doBy (Søren Højsgaard et al., 2014), plyr (Wickham, 2011), reshape2 (Wickham, 2007) and the npde (Comets et al., 2008) packages. Morphine sulphate doses and metabolite concentrations were converted to their morphine base equivalent after fitting in NONMEM using a conversion factor of 0.752 for morphine sulphate, and 0.618 for the two morphine glucuronides. Missing concentrations and pain score data were omitted from the analyses.

2.2.2 Literature models

Four previously published reports of population PK analyses of morphine in neonates (Anand et al., 2008; Bouwmeester et al., 2004; Knibbe et al., 2009; Wang et al., 2013) were selected to provide the basis for a systematic model comparison. The models were selected based upon a literature search, and by scanning reference lists of relevant articles for additional studies. The performance of these models and their suitability for simulation purposes was assessed in the current study using the new dataset. The four studies contained five population PK models, two of which were based solely on morphine plasma concentrations, one was based on morphine and M3G concentrations, and two were based on morphine, M3G, and M6G concentrations. The characteristics of the five models are shown in Table 1.

Table 1
Details of the five population pharmacokinetic models of morphine

The evaluation of the literature models was performed in a sequential manner for the parent drug and then the metabolites. First, the morphine PK for all five models was evaluated, including only the morphine plasma concentrations and the morphine parameters. Next, the three models with metabolites were evaluated.

2.2.2.1 Comparing literature models with new data

Individual posthoc parameter values, objective function values (OBJ), predictions, and residuals were computed for each model using the new dataset as the input. The models were coded in NONMEM and run without fitting the models (MAXEVAL=0 in the estimation step) with the values specified for each model as the initial parameter estimates (Owen and Fiedler-Kelly, 2014).

2.2.2.2 Tools used for systematic model comparison
  1. The first metric used to compare the models was OBJ, which was used to rank the models in terms of finding the best fit for the new dataset (Mould and Upton, 2013). To account for the different numbers of parameters in the models, the Akaike Information Criterion (AIC) was also compared.
  2. All models were examined with respect to their ability to describe clearance for individuals in the new dataset. Distribution densities of Individual Empirical Bayesian posthoc ηclearance values were compared to the theoretical η-distribution N(0,ω2), where ω2 was the reported value of ηclearance for the respective model. ηclearance describes the individual’s difference from the population typical value. It was hypothesised that if a model adequately described the new dataset, the posthoc estimates for the ηs would imitate the theoretical η-distribution and η-shrinkage would be low (Savic and Karlsson, 2009). The overlap between the two distributions was visually examined and interpreted in context of the η-shrinkage.
  3. To evaluate the predictive performance of the models, 200 simulations using the current dataset as input were performed to obtain model-based predicted concentrations of the dataset with each model. The simulations were evaluated by Visual Predictive Checks (VPCs) (Holford, 2005; Karlsson and Holford, 2008).
  4. Normalised prediction distribution errors (NPDEs) were used to compare the simulated model-predicted concentrations to the observations in the current dataset (Comets et al., 2008). A Wilcoxon signed rank t-test, a Fischer test for variance, and a Shapiro-Wilks (SW) test was used to test whether the NPDEs followed a normal distribution. Furthermore, a global p-value was reported (Brendel et al., 2006). The models were compared and ranked according to p-values from the various tests.

2.2.3 Development of a pharmacokinetic meta-model

The parent drug model and the metabolite model that performed best in the systematic model comparison stage were carried forward and combined to propose a new meta-model for morphine and metabolite PK. The compartmental structure and covariate relationships for morphine PK, was combined with the compartmental structure and covariate relationships for metabolite PK. Morphine plasma concentrations were initially fitted alone, and subsequently simultaneously with M3G and M6G plasma concentrations to update all structural parameter estimates for the combined model.

2.2.3.1 Updating parameter estimates

Different strategies for updating the parameter estimates in the meta-model were investigated. First, the model was fitted to the concentrations in the new dataset with FOCE-I estimation, using the parameter values from the published model as initial parameter estimates. Secondly, the model was fitted using the $PRIOR subroutine in NONMEM to stabilize estimation, and help the model converge. The prior value was selected from the model that had the best predictive performance in the systematic model comparison, and priors were tested both on the structural parameters (θ) for clearance and volume, and on between-subject variability (BSV) parameters (ω). Incorporating a prior value on a θ required the variance of that prior parameter estimate to reflect uncertainty about the prior. For priors on ω’s, the number of subjects in the prior study was used as degrees of freedom for the prior ω.

The structural parameters in the parent drug model and the metabolite model were updated to describe the new dataset, whilst the covariate parameters were fixed to those values obtained in the previous studies.

2.2.3.2 Stochastic models

The structure of the stochastic components of the models were also investigated when updating parameters and developing the meta-model. Between-subject variability of the structural parameters was described using an exponential model as per the original models. Omega blocks for random effects were tested if graphical inspection of η-scatter plots showed correlation. A combined additive and proportional error model was applied to describe the residual unexplained variability (RUV). This error model was reduced to either additive or proportional if an error term was very small, or the removal did not show significantly impact on the model, or if removing the error term increased stability of the model. Different error terms were tested for morphine, M3G, and M6G concentrations.

2.2.4 Pharmacodynamic data analysis

A sequential PK/PD modelling approach using the PPP&D method was employed (Zhang et al., 2003). The PD model used the final pharmacokinetic meta-model to define input concentrations. The concentration-analgesia relationship for morphine plasma concentrations and the PIPP pain measure was described as baseline effect with either additive or proportional morphine effect.

Baseline effect was tested with or without BSV either normally or log normally distributed. To describe the morphine effect, structural models of linear, log-linear, Emax, and sigmoidal Emax functions of PIPP score versus either plasma concentration (direct effect) or effect compartment concentrations (delayed effect) was tested. The equilibration rate constant between the central compartment and the effect compartment, ke0, was estimated or fixed to a literature value (Bouw et al., 2000; Sverrisdottir et al., 2014). Normally or log normally distributed BSV was tested as appropriate for the different structural model parameters. The drug effect was tested for both morphine and M6G concentrations, either alone or combined. Drug effect was not tested for M3G. When a combined morphine and M6G effect was tested, M6G potency was described as a fraction of morphine response, which was also estimated or fixed to a literature value (Romberg et al., 2003). The influence of covariates was assessed graphically with plots of η-estimates versus relevant covariates after selecting the base model candidate. To ensure this step was unbiased, shrinkage of ηs had to meet the criteria described in section 2.2.5 (Savic and Karlsson, 2009). Covariates were tested in NONMEM if the graphical assessment showed any trends.

2.2.5 Model selection and evaluation

The PK and PD models were evaluated by a range of goodness-of-fit criteria. These criteria included η-shrinkage < 30% to assess the validity of the individual parameter estimates, η p-values > 0.1 to assess that ηs were distributed around zero, and standard error of θs < 30% and ηs < 50%. For addition of a single parameter in nested models, a decrease of 3.84 units of the OBJ (corresponding to a p-value < 0.05 according to the Chi2 distribution) was required, and for non-nested PD models, a decrease of 2 units in AIC was used to select one model over another. The model performance were also assessed with diagnostic plots, and the parameter precision and stability of the models were evaluated through a bootstrap analysis (n=1000). The predictive performance was evaluated by VPCs and NPDEs.

2.2.6 Simulations with the pharmacokinetic meta-model

The pharmacokinetic meta-model was coded in the R programming language (R Development Core Team, 2014), and made available for simulations through an open-source web-application using the Shiny package (RStudio Inc, 2014; Wojciechowski et al., 2015) in Rstudio (RStudio, 2009–2013). The population model was coded with differential equations, using deSolve (Soetaert et al., 2010).

3 RESULTS

3.1 Characteristics of the study population

A total of 30 preterm neonates comprised the study population. Data from 27 subjects were used in the PK analysis, and data from all 30 subjects were used in the PD analysis. Demographic data and a summary of the dataset are given in Table 2. Two subjects received morphine for less than 24 hours, but were still included in the analysis. Five infants died during or after the study. All deaths were due to complications of extreme prematurity and unrelated to the study, and data from these subjects were included in the pharmacometric analysis.

Table 2
Demographic data and summary of dataset

3.2 Systematic comparison of literature models

3.2.1 Comparing objective function value

The OBJ computed by NONMEM, when the five parent drug models were applied to the new dataset without fitting, showed a difference between all the models of more than 3.84 units. Model 1 (Anand et al., 2008) was significantly superior to the other four models in terms of the best fit based on OBJ. When calculating the AIC values for the model, the ranking of the models remained unchanged (Table 3).

Table 3
Objective function values for running the five models with MAXEVAL=0

The OBJ or AIC values for the metabolite models were not comparable, because the models were based on different numbers of observations as model 3 only included M3G concentrations, whereas model 4 and model 5 included both M3G and M6G concentrations.

3.2.2 Posthoc estimates of clearance variability

The distribution densities of Individual Empirical Bayesian posthoc η estimates for morphine clearance were evaluated for the five parent drug models. The mean for the posthoc ηs were less than zero for all five models, but the deviation was less for model 5. For all five models, the posthoc ηs were different from the theoretical η-distribution for the published models, and model 2, 3, 4 and 5 all demonstrated a distorted shape of the distribution. Model 4 had very high η-shrinkage (>30%), suggesting that the individual estimates of clearance shrunk towards the model’s population value of morphine clearance. The distribution density of Individual Empirical Bayesian posthoc η estimates was also compared to the theoretical distribution density (based on the reported omega-value) for metabolite elimination clearance for the three metabolite models. The mean of the posthoc estimates was markedly different from zero for model 3, and model 5 demonstrated a distorted shape of the distribution for both M3G and M6G. Model 4 had the best overlay between the posthoc distribution and the theoretical distribution. (For plots, see Supplementary material)

3.2.3 Predictive performance

VPCs were used to evaluate the models’ abilities to predict the observed concentration-time data. The VPCs for the parent drug models indicated some minor misspecifications, but did not clearly differentiate between the five models. The VPCs also showed misspecifications for all three metabolite models, but it was not possible to rank the models based on this tool. (The VPCs are presented in supplementary material)

The NPDEs were to a greater extent able to differentiate between the parent drug models, and the results of the morphine NPDE analysis are presented in Figure 2. The global p-value was less than 0.05 for all 5 models, but model 1 overall demonstrated the best simulation performance based on these test results.

Figure 2
Distribution of NPDEs (histogram) with theoretical N(0,1) distribution (dashed line) for morphine predictions. (*) indicates that the NPDEs are different from N(0,1) distribution for the specified test on a 95% significance level. A) Model 1 Anand et ...

The NPDEs of the three metabolite models were analysed separately for each of the metabolites in the models. Model 4 predicted the M6G concentrations well, and also showed the best performance of the three models with regards to the M3G concentrations. The visualisation of the NPDE distribution for all three metabolite models is shown in Figure 3.

Figure 3
Distribution of NPDEs (histogram) with theoretical N(0,1) distribution (dashed line) for metabolite predictions. (*) indicates that the NPDEs are different from N(0,1) distribution for the specified test on a 95% significance level. A) Model 3 Wang et ...

Overall, model 1 showed the best predictive performance for morphine concentrations based on OBJ, VPCs, and NPDEs, and this model was chosen as the basis for updating the parameter estimates for the parent drug pharmacokinetics and proposing a new pharmacokinetic meta-model.

Even though none of the three metabolite models predicted the concentrations of M3G and M6G without bias, model 4 showed the overall best predictive performance. The covariates of this model also accounted for maturation of metabolite elimination clearance, so this metabolite model was used for the pharmacokinetic meta-model.

3.3 Proposal of a pharmacokinetic meta model

3.3.1 Structure of the final pharmacokinetic meta-model

The structure of the parent drug model that performed best in the systematic model comparison (model 1) was combined with the structure of the metabolite model that performed best (model 4). The structure of the final model is described in detail here for completeness, although it incorporates elements of the published models. A one-compartment structural model for morphine was parameterised as total morphine clearance and volume of distribution, and combined with a compartment for each of the metabolites. The formation of metabolites was estimated as a fraction of the total morphine clearance. A schematic representation of the model is shown in Figure 4.

Figure 4
Structural pharmacokinetic meta-model of morphine (MOR) and its metabolites (M3G, M6G). The blue colour denotes components from the morphine literature model, and the red colour denotes components from the metabolite literature model. All parameter estimates ...

Covariates were included a priori in the meta-model, from the literature models. All parameter values were scaled for size with an allometric power model, as shown in equation (1).

θi=θstd(WTi70kg)PWR
(1)

where θi is the parameter value of the ith individual, θstd is the parameter value of a standardised adult with a bodyweight of 70 kg and WTi is the bodyweight of the ith individual. The exponent (PWR) was fixed to 0.75 for clearance terms, and 1 for volume terms.

Morphine is excreted through hepatic pathways, and to account for the maturation of morphine clearance, a sigmoidal Hill equation was used to incorporate PMA as a covariate. Furthermore a scaling parameter was applied to morphine clearance if the neonates was preterm, in order to be able to compare with clearance in term neonates (Anand et al., 2008). The relationship between clearance and PMA is shown in equation (2).

CLMORPOP=CLMORSTDPMAHillPMAHill+TM50HillFDEVCL
(2)

where CLMORPOP is clearance in a subpopulation with a given PMA, CLMORSTD is the population estimate (standardised to 70 kg), PMA is the postmenstrual age in weeks (defined as the post conception age at birth plus the PNA of the neonate), TM50 is the PMA at which clearance is 50% that of the mature value and Hill is the coefficient that describes the slope of the sigmoidal curve. FDEVCL is a constant used to scale clearance for premature neonates. A similar scaling parameter, FDEVV, was applied to the volume of the morphine compartment.

To describe age-related changes in the clearance of the metabolites which are excreted renally, an exponential maturation function (Bouwmeester et al., 2004) was applied to the metabolite elimination clearance terms. This relationship is shown in equation (3).

CLEMxGPOP=CLEMxGstd(1βrfexp(PNALn(2)Trf))
(3)

where CLEMxGPOP is the elimination clearance of either M3G or M6G in a subpopulation with a given PNA, CLEMxGstd is the population estimate (standardised to 70 kg), PNA is the postnatal age in days, βrf describe the fraction below the mature value, and Trf is the maturation half-life of the metabolite elimination clearance.

3.3.2 Parameter estimation

For updating the parameter estimates of the parent drug PK, different combinations with and without priors on the structural parameters of the model, and on the BSV terms were tried using the $PRIOR subroutine in NONMEM. The results of this analysis are shown in Table 4, where the models are listed according to their OBJ, and it can be seen how priors affect the parameter estimates.

Table 4
Objective function values and parameter estimates for screened combinations of priors

Having a prior on the volume of distribution for the morphine compartment and the BSV on volume (V), gave the best fit for the parent drug PK. These priors were found necessary because the current dataset was based on a continuous infusion of morphine, with very few PK samples after discontinuation of infusion, and therefore the data did not have enough information to provide a precise estimate of V. The prior value was selected from the model that performed best in the external evaluation, i.e. model 1 (Anand et al., 2008). Prior parameter estimates from the other literature models were not used, because of the poor ability of these models to describe the present data. For the metabolite PK, the elimination clearance and the metabolite formation fraction was estimated using FOCE-I. The metabolite volumes of distribution were not possible to estimate with the present study design and these two parameters were therefore fixed at values from the selected metabolite population PK model (Bouwmeester et al., 2004).

3.3.3 Description and validation of the final meta-model

The proposed model described the concentration time data for morphine and the two glucuronide metabolites, M3G and M6G, accurately. The distribution density of Individual Empirical Bayesian posthoc η estimates for morphine clearance overlaid the model predicted distribution with good agreement, and the mean of the posthoc ηclearance is zero (See supplementary material). The model predicted the individual clearance estimates well, with a η-shrinkage of only 3.9% (shrinkage for other parameters was between 6.8–11.4%). There was a very strong correlation between elimination clearances for the two metabolites, which is explained by their mutual excretion pathways. The goodness-of-fit plots for the meta-model, shown in Figure 5, revealed no significant trends in the residuals, for either morphine, M3G, or M6G. Parameter precision for all structural parameters, all variance terms, and all residual error terms was generally good (%SE < 20%). All parameter estimates for the meta-model are shown in Table 5.

Figure 5
Diagnostic plots for the final pharmacokinetic meta-model. The black solid lines are identity lines with a slope of 1 for plots A) and B) and a slope of 0 for plots C) and D). The black dashed lines are regression lines for the data. The red symbols represent ...
Table 5
Parameter estimates for the pharmacometric models

The PK model showed acceptable predictive performance when assessed by VPCs (Figure 6A). The model predictions overlay the observed concentrations with good agreement. The NPDE analysis of the final pharmacokinetic meta-model generally showed that the model was able to predict concentrations well (Figure 6B).

Figure 6
A: Visual Predictive Check of the final meta-model of morphine and metabolite PK. The observed data are represented by blue symbols, a red solid line (median), and red dashed lines (5th and 95th percentiles). The red shaded area represents the 95% empirical ...

The meta-model is available for simulations as a shiny web application at this URL: https://unicph.shinyapps.io/MorphineNeonates/

3.4 Pharmacodynamic model

The structural model for analgesia-concentration relationship that obtained the lowest AIC in the analysis of the PIPP score was linear and additive to baseline. The effect was delayed, and a model with an effect compartment with the equilibration half time fixed to 0.5 hours (Bouw et al., 2000; Sverrisdottir et al., 2014) was statistically superior to a model with an equilibration half time fixed to 0.28 hours (P = 0.014) (Inturrisi and Colburn, 1986), or a model without an effect compartment (P = 0.026). A model with only morphine concentrations linked to analgesia was superior to models with M6G concentrations (P < 0.001), or a combination of morphine and M6G concentrations (P < 0.001) linked to effect. Parameter estimates for the PK/PD model are given in Table 5, and a VPC for the PD model is shown in Figure 7.

Figure 7
Visual Predictive Check of PD model for PIPP score. The observed data are represented by blue symbols, a red solid line (median), and red dashed lines (5th and 95th percentiles). The red shaded area represents the 95% empirical confidence interval around ...

The model that best described the data had a positive population parameter estimate of the slope, corresponding to an increase in the PIPP score of 0.22 units from baseline (increase in pain), at a clinically relevant concentration of 20 ng/mL (Bouwmeester et al., 2003). However, it should be noted that the additive BSV on this parameter allows the slope to take negative values for some subjects which corresponds to pain relief when dosed with morphine. This range of minor positive and negative slopes (95% CI: −0.030 – 0.052) across the population indicates that there is no clear concentration-effect relationship for these data across the population, with some subjects tending to have increased pain and others experiencing analgesia on average over the study period. The high standard error on the slope and on the variability parameters confirms this, as it was not possible to obtain precise parameter estimates.

4 DISCUSSION

The pharmacokinetics of morphine, M3G, and M6G in preterm neonates was described by developing a meta-model that accounted for the maturational changes in clearance of morphine and its metabolites. The best structural and stochastic features of previously published literature models were combined to propose this meta-model, and parameter estimates were updated to describe the new dataset. The concentration-response relationship was investigated for the analgesic effect of morphine, but it was not possible to establish a concentration-effect relationship using PIPP scores as endpoint in this clinical setting.

4.1 Comparison of literature models

Several population PK models of morphine in paediatrics have been published over the last decade. Five of these models were compared and evaluated using a novel dataset. Physiological processes change rapidly in preterm neonates, and therefore it is important to adequately describe these changes if such a model is to be used as a tool for prediction and simulation purposes. Maturation of the liver has significant impact on glucuronidation processes (Alcorn and McNamara, 2002; de Wildt et al., 1999; Pacifici et al., 1982), thus affecting the hepatic clearance of morphine. Kidney maturation is also a well-known phenomenon (Bueva and Guignard, 1994; Mannan et al., 2012; van den Anker, 1996), which affects the renal clearance of morphine metabolites.

Model 1 (Anand et al., 2008), incorporating body weight, a sigmoidal maturation model as a function of PMA, and a scaling factor to differentiate between term and preterm neonates, had the best predictive performance for morphine. Model 4 (Bouwmeester et al., 2004) showed the poorest predictive performance of morphine concentrations. However, this model was the best in terms of predicting the metabolite concentrations. This might suggest that the exponential function of PNA incorporated in this model is not descriptive for the maturation of morphine clearance in this population, but it describes the age-related changes in metabolite clearance well.

The weight related changes in clearance were described by a power-model with an estimated exponent of 1.44 in model 5 (Knibbe et al., 2009). This model was not able to predict the concentrations in the new dataset, suggesting that the single power function did not describe the maturational changes in morphine and metabolite PK in this population. The model accounted for high and low glucuronidation capacity, by dividing the infants into two age groups, suggesting a sudden increase in clearance approximately 10 days after birth. Including age-related changes as a continuous factor is however more physiologically plausible (Alcorn and McNamara, 2002).

The inability of the literature models to describe the information in the new dataset does not necessarily mean that these models are inadequate models. Part of the explanation can be due to differences in the population, since all these models were derived from index datasets that were slightly different from the current dataset. It was therefore found necessary to propose a new meta-model that combined the best structural and stochastic characteristics of the literature models, to give a better description of morphine and its metabolite pharmacokinetics of the population in the new dataset.

4.2 Pharmacokinetic meta model

The published models failed to sufficiently describe the new population’s time-courses of morphine and its metabolite concentrations, resulting in the need to estimate parameters using the observed dataset. The results of the estimation analysis incorporating prior information showed that having a prior of the morphine volume of distribution (V) gave the best fit. This work confirms that the prior functionality of NONMEM enabled a method for dealing with sparse and/or uninformative data.

The maturation model implemented in the meta-model allowed clearance to reach 78% of adult clearance at a PNA of 1 year for an infant born after 23 weeks of gestation, and 85% of adult clearance at a PNA of 1 year for an infant born after 32 weeks of gestation. The influence of the maturation model becomes negligible for PNAs of more than 1 year, and the increase in clearance is only dependent on bodyweight after this age.

The relative M3G formation accounted for 72.9% of morphine clearance, and relative formation of M6G accounted for 27.1% of morphine clearance. The ratio between the two metabolites is consistent with values previously reported (M3G: 45–85%, M6G: 10–15%; Barrett et al., 1996; Christrup, 1997). However, the values for the metabolite volume of distribution was set to a previously published value, so this formation ratio reflects apparent values, normalised for the true unknown volumes. Clearance through other pathways than glucuronidation, such as sulphation, was not examined since the current data only provided information about the glucuronide metabolites.

BSV on the structural parameters is relatively high for the meta-model. This was expected because the study was performed during routine clinical practice, which increased the variability in dosing and sampling schemes. Underlying diseases in these neonates makes the population very heterogeneous. Hydration state is not stable in this population, which affects renal function and hence clearance mechanisms (Roy and Sinclair, 1975). Use of concomitant medications, such as anaesthetics could have affected morphine pharmacokinetics (Allegaert et al., 2007). These factors are all expected to contribute to the variability in this population.

Nonetheless, the new meta-model was able to describe and predict morphine, M3G, and M6G concentrations for the population included in this study as determined by the evaluation criteria and predictive performance of the model. The model effectively described the effect of bodyweight and the influence of age-related changes on both morphine and metabolite clearance.

The meta-model potentially provides a method for supporting decisions about dosing regimens and dose adjustments. This meta-model takes into account features from previously published PK models. Hence a systematic covariate analysis in this particular situation is different to the norm, and a full covariate screening was not performed in this study. However, there is a risk of over-parameterising a model by introducing covariates based on assumptions, and not by formally testing them. To determine whether this model is representative of other populations of preterm neonates, external validation is needed. Nonetheless, given that it is a meta-model, this might be of less concern than if it were conventionally derived from modelling the current dataset alone, since the model encloses properties already validated in other populations.

4.3 Pharmacodynamic analysis

Despite peak morphine plasma concentrations that are high enough to be in the range expected to cause pain relief (Chay et al., 1992) it was not possible to establish a clear concentration-effect relationship in this study. This is consistent with previous studies that also failed to establish the concentration-effect relationship for morphine in neonates (Anand et al., 2008; Cignacco et al., 2008).

The absence of placebo data was expected to be a confounding issue obstructing the determination of the drug effect in this study. Previous studies with morphine and other opioids in experimental pain settings have demonstrated the importance of detecting placebo effect (Juul et al., 2014; Ravn et al., 2014; Sverrisdottir et al., 2014), which is generally only possible in an experimental pain setting. These studies showed that it is essential to understand the placebo response and the changes in baseline over time in order to identify the true drug effect. However, there are several ethical concerns associated with conducting a clinical trial and the use of placebo, such as the discomfort of not receiving pain relief, the possible long-term impact of untreated pain in neonates, and issues regarding consent for participation in placebo-controlled trials (Cignacco et al., 2008; Walker, 2013).

The pain in this study was very unstable, and the neonates had multiple different serious clinical conditions, which make it impossible to determine whether the changes in pain score were related to progression of their illness or to the effect of morphine. Without an underlying placebo model, describing what would happen in the absence of morphine, and thus accounting for the highly variable trajectory or the clinical condition including pain, it was found to be difficult to identify a plausible concentration-effect relationship for morphine in this clinical setting.

Measuring pain in neonates and the validity of tools used for pain assessment has been debated extensively (Mattsson et al., 2011; Ramelet et al., 2004; Stevens et al., 1995). Neonates are by nature unable to communicate their experience of pain, and pain assessment typically includes subjective observations. Difficulties in differentiating between pain and other constructs such as stress, agitation, and sedation, and insufficient knowledge about how severe or critical illness affects neonates’ signs of pain are issues that contribute to the challenges of managing pain in neonates. The interpretations of the PIPP score can be confounded by sedation, because this score includes behavioural indicators, such as sleep and body movements. The estimated baseline score is 5.94 on the PIPP scale, which means the neonates have minimal or no pain (Stevens et al., 1996). This may provide another explanation for the difficulties in detecting the concentration-effect relationship; the pain score was at a level where only minor decrease was possible. However, morphine’s analgesic properties generally occur at lower concentrations than sedation (Graves et al., 1985; Kart et al., 1997), and the neonates were potentially given a dose that was high enough to produce sedation and obliterate all pain (Chay et al., 1992). Therefore, the PIPP score might reflect sedation rather than pain, or fail to characterize pain in the presence of sedation.

Physiological parameters used as indicators for pain, such as blood pressure, respiratory rate, and heart rate, have previously been found unreliable in detecting pain in neonates (Büttner and Finke, 2000). Identification and use of objective biomarkers for the analgesic effect of morphine, e.g. prolactin concentration and pupil diameter (Brokjaer et al., 2015; Staahl et al., 2011), could possibly be a step towards describing the concentration-effect relationship of morphine in neonates to improve individual drug therapy (Anderson and van den Anker, 2014).

5 CONCLUSION

Five published models of morphine pharmacokinetics in neonates were compared using a new dataset as input. The published models were unable to predict morphine and metabolite PK profiles without bias. Therefore, structural and stochastic characteristics of the model that had the best predictive performance with respect to morphine were combined with structural and stochastic characteristics of the best metabolite model. Parameter estimates were updated to propose a new meta-model of the pharmacokinetics of morphine and its metabolites in preterm neonates. The model implements maturation models that describe the clearance of morphine, including the formation of M3G and M6G, and the elimination of the morphine metabolites. This study did not suggest a significant relationship between plasma concentrations and the analgesic effect of morphine, using PIPP score as the effect measure. This is most likely the result of the underlying highly variable trajectory of the clinical conditions and resultant pain experienced by the preterm neonate.

Supplementary Material

Acknowledgments

David J.R. Foster and Richard N. Upton acknowledge that the Australian Centre for Pharmacometrics is an initiative of the Australian Government as part of the National Collaborative Research Infrastructure Strategy.

John van den Anker was supported by a NIH grant (R01HD048689) to conduct this neonatal study; Patients were included at two institutions (Kosair Children’s Hospital in Louisville, Kentucky and Children’s National Health System in Washington, DC) and Drs. Jan Sullivan and Williams have been instrumental for the conduct of the study

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

  • Alcorn J, McNamara PJ. Ontogeny of hepatic and renal systemic clearance pathways in infants: part I. Clinical pharmacokinetics. 2002;41:959–998. [PubMed]
  • Allegaert K, Peeters MY, Verbesselt R, Tibboel D, Naulaers G, de Hoon JN, Knibbe CA. Inter-individual variability in propofol pharmacokinetics in preterm and term neonates. British journal of anaesthesia. 2007;99:864–870. [PubMed]
  • Anand KJ, Anderson BJ, Holford NH, Hall RW, Young T, Shephard B, Desai NS, Barton BA, Group NTI. Morphine pharmacokinetics and pharmacodynamics in preterm and term neonates: secondary results from the NEOPAIN trial. British journal of anaesthesia. 2008;101:680–689. [PMC free article] [PubMed]
  • Anderson BJ, van den Anker J. Why is there no morphine concentration-response curve for acute pain? Paediatric anaesthesia. 2014;24:233–238. [PubMed]
  • Barrett D, Barker D, Rutter N, Pawula M, Shaw P. Morphine, morphine-6-glucuronide and morphine-3-glucuronide pharmacokinetics in newborn infants receiving diamorphine infusions. British journal of clinical pharmacology. 1996;41:531–537. [PMC free article] [PubMed]
  • Beal S, Sheiner LB, Boeckmann A, Bauer RJ. NONMEM User’s Guides (1989–2009) ICON Development Solutions; Ellicott City, MD, USA: 2009.
  • Bhalla T, Shepherd E, Tobias JD. Neonatal pain management. Saudi journal of anaesthesia. 2014;8:S89–97. [PMC free article] [PubMed]
  • Bouw MR, Gårdmark M, Hammarlund-Udenaes M. Pharmacokinetic-Pharmacodynamic Modelling of Morphine Transport Across the Blood-Brain Barrier as a Cause of the Antinociceptive Effect Delay in Rats—A Microdialysis Study. Pharmaceutical Research. 2000;17:1220–1227. [PubMed]
  • Bouwmeester N, Van Den Anker J, Hop W, Anand K, Tibboel D. Age - and therapy - related effects on morphine requirements and plasma concentrations of morphine and its metabolites in postoperative infants. British journal of anaesthesia. 2003;90:642–652. [PubMed]
  • Bouwmeester NJ, Anderson BJ, Tibboel D, Holford NH. Developmental pharmacokinetics of morphine and its metabolites in neonates, infants and young children. British journal of anaesthesia. 2004;92:208–217. [PubMed]
  • Brendel K, Comets E, Laffont C, Laveille C, Mentré F. Metrics for external model evaluation with an application to the population pharmacokinetics of gliclazide. Pharmaceutical research. 2006;23:2036–2049. [PMC free article] [PubMed]
  • Brokjaer A, Olesen AE, Kreilgaard M, Graversen C, Gram M, Christrup LL, Dahan A, Drewes AM. Objective markers of the analgesic response to morphine in experimental pain research. Journal of pharmacological and toxicological methods. 2015;73C:7–14. [PubMed]
  • Bueva A, Guignard JP. Renal function in preterm neonates. Pediatric research. 1994;36:572–577. [PubMed]
  • Büttner W, Finke W. Analysis of behavioural and physiological parameters for the assessment of postoperative analgesic demand in newborns, infants and young children: a comprehensive report on seven consecutive studies. Pediatric Anesthesia. 2000;10:303–318. [PubMed]
  • Canadian Paediatric Society. Prevention and management of pain and stress in the neonate. In: American Academy of Pediatrics, Committee on Fetus and Newborn, Committee on Drugs, Section on Anesthesiolog, Section on Surgery, Fetus and Newborn Committee, editor. Pediatrics. 2000. pp. 454–461. 2000/02/02 ed. [PubMed]
  • Carbajal R, Lenclen R, Jugie M, Paupe A, Barton BA, Anand KJ. Morphine does not provide adequate analgesia for acute procedural pain among preterm neonates. Pediatrics. 2005;115:1494–1500. [PubMed]
  • Chay PC, Duffy BJ, Walker JS. Pharmacokinetic-pharmacodynamic relationships of morphine in neonates. Clinical pharmacology and therapeutics. 1992;51:334–342. [PubMed]
  • Choonara IA, McKay P, Hain R, Rane A. Morphine metabolism in children. British journal of clinical pharmacology. 1989;28:599–604. [PMC free article] [PubMed]
  • Christrup LL. Morphine metabolites. Acta anaesthesiologica Scandinavica. 1997;41:116–122. [PubMed]
  • Cignacco E, Hamers JP, van Lingen RA, Zimmermann LJ, Muller R, Gessler P, Nelle M. Pain relief in ventilated preterms during endotracheal suctioning: a randomized controlled trial. Swiss medical weekly. 2008;138:635–645. [PubMed]
  • Comets E, Brendel K, Mentré F. Computing normalised prediction distribution errors to evaluate nonlinear mixed-effect models: the npde add-on package for R. Computer methods and programs in biomedicine. 2008;90:154–166. [PubMed]
  • de Wildt SN, Kearns GL, Leeder JS, van den Anker JN. Glucuronidation in humans. Pharmacogenetic and developmental aspects. Clinical Pharmacokinetics. 1999;36:439–452. [PubMed]
  • Graves DA, Arrigo JM, Foster TS, Baumann TJ, Batenhorst RL. Relationship between plasma morphine concentrations and pharmacologic effects in postoperative patients using patient-controlled analgesia. Clinical pharmacy. 1985;4:41–47. [PubMed]
  • Holford N. The Visual Predictive Check Superiority to Standard Diagnostic (Rorschach) Plots. 14th Meeting of the Population Approach Group in Europe; Pamplona, Spain. 2005.
  • Inturrisi C, Colburn W. Application of pharmacokinetic-pharmacodynamic modelling to analgesia. In: Foley K, Inturrisi C, editors. Advances in pain research and therapy. Opioid analgesics in the management of clinical pain. Raven Press; New York: 1986. pp. 441–452.
  • Johnston CC, Stevens BJ. Experience in a neonatal intensive care unit affects pain response. Pediatrics. 1996;98:925–930. [PubMed]
  • Juul RV, Foster DJ, Upton RN, Andresen T, Graversen C, Drewes AM, Christrup LL, Kreilgaard M. Pharmacodynamic modelling of placebo and buprenorphine effects on event-related potentials in experimental pain. Basic & clinical pharmacology & toxicology. 2014;115:343–351. [PubMed]
  • Karlsson MO, Holford NH. A Tutorial on Visual Predictive Checks. 17th Annual Meeting of the Population Approach Group in Europe; Marseille, France. 2008.
  • Kart T, Christrup LL, Rasmussen M. Recommended use of morphine in neonates, infants and children based on a literature review: part 2–clinical use. Pediatric Anesthesia. 1997;7:93–101. [PubMed]
  • Kilpatrick GJ, Smith TW. Morphine-6-glucuronide: actions and mechanisms. Medicinal research reviews. 2005;25:521–544. [PubMed]
  • Knibbe CA, Krekels EH, van den Anker JN, DeJongh J, Santen GW, van Dijk M, Simons SH, van Lingen RA, Jacqz-Aigrain EM, Danhof M. Morphine glucuronidation in preterm neonates, infants and children younger than 3 years. Clinical pharmacokinetics. 2009;48:371–385. [PubMed]
  • Mannan MA, Shahidulla M, Salam F, Alam MS, Hossain MA, Hossain M. Postnatal development of renal function in preterm and term neonates. Mymensingh medical journal : MMJ. 2012;21:103–108. [PubMed]
  • Mattsson JY, Forsner M, Arman M. Uncovering pain in critically ill non-verbal children: nurses’ clinical experiences in the paediatric intensive care unit. Journal of child health care : for professionals working with children in the hospital and community. 2011;15:187–198. [PubMed]
  • Mazoit JX, Butscher K, Samii K. Morphine in postoperative patients: pharmacokinetics and pharmacodynamics of metabolites. Anesthesia & Analgesia. 2007;105:70–78. [PubMed]
  • Meng QC, Cepeda MS, Kramer T, Zou H, Matoka DJ, Farrar J. High-performance liquid chromatographic determination of morphine and its 3- and 6-glucuronide metabolites by two-step solid-phase extraction. Journal of chromatography. B, Biomedical sciences and applications. 2000;742:115–123. [PubMed]
  • Mould DR, Upton RN. Basic concepts in population modeling, simulation, and model-based drug development-part 2: introduction to pharmacokinetic modeling methods. CPT: pharmacometrics & systems pharmacology. 2013;2:e38. [PMC free article] [PubMed]
  • Owen JS, Fiedler-Kelly J. NONMEM Overview and Writing an NM-TRAN Control Stream, Introduction to Population Pharmacokinetic/Pharmacodynamic Analysis with Nonlinear Mixed Effects Models. John Wiley & Sons, Inc; 2014. pp. 28–65.
  • Pacifici GM, Sawe J, Kager L, Rane A. Morphine glucuronidation in human fetal and adult liver. European journal of clinical pharmacology. 1982;22:553–558. [PubMed]
  • Penson RT, Joel SP, Clark S, Gloyne A, Slevin ML. Limited phase I study of morphine-3-glucuronide. Journal of Pharmaceutical Sciences. 2001;90:1810–1816. [PubMed]
  • R Development Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing; Vienna, Austria: 2014.
  • Ramelet AS, Abu-Saad HH, Rees N, McDonald S. The challenges of pain measurement in critically ill young children: a comprehensive review. Australian critical care : official journal of the Confederation of Australian Critical Care Nurses. 2004;17:33–45. [PubMed]
  • Ravn P, Foster DJ, Kreilgaard M, Christrup L, Werner MU, Secher EL, Skram U, Upton R. Pharmacokinetic-pharmacodynamic modelling of the analgesic and antihyperalgesic effects of morphine after intravenous infusion in human volunteers. Basic & clinical pharmacology & toxicology. 2014;115:257–267. [PubMed]
  • Romberg R, Olofsen E, Sarton E, den Hartigh J, Taschner PE, Dahan A. Pharmacokinetic-pharmacodynamic modeling of morphine-6-glucuronide-induced analgesia in healthy volunteers: absence of sex differences. Anesthesiology. 2004;100:120–133. [PubMed]
  • Romberg R, Olofsen E, Sarton E, Teppema L, Dahan A. Pharmacodynamic effect of morphine-6-glucuronide versus morphine on hypoxic and hypercapnic breathing in healthy volunteers. Anesthesiology. 2003;99:788–798. [PubMed]
  • Roy RN, Sinclair JC. Hydration of the low birth-weight infant. Clinics in perinatology. 1975;2:393–417. [PubMed]
  • RStudio. RStudio: Integrated development environment for R. RStudio, Inc; Boston, MA: 2009–2013. 0.98.1028 ed.
  • RStudio Inc. shiny: Web Application Framework for R, R package version 0.11.1. 2014 CRAN.R-project.org.
  • Savic RM, Karlsson MO. Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions. The AAPS journal. 2009;11:558–569. [PMC free article] [PubMed]
  • Simons SH, van Dijk M, Anand KS, Roofthooft D, van Lingen RA, Tibboel D. Do we still hurt newborn babies? A prospective study of procedural pain and analgesia in neonates. Archives of pediatrics & adolescent medicine. 2003;157:1058–1064. [PubMed]
  • Soetaert K, Petzoldt T, Setzer RW. Solving differential equations in R: package deSolve. Journal of Statistical Software. 2010;33
  • Stevens B, Johnston C, Petryshen P, Taddio A. Premature Infant Pain Profile: development and initial validation. The Clinical journal of pain. 1996;12:13–22. [PubMed]
  • Stevens BJ, Johnston CC, Grunau RV. Issues of assessment of pain and discomfort in neonates. Journal of obstetric, gynecologic, and neonatal nursing : JOGNN/NAACOG. 1995;24:849–855. [PubMed]
  • Staahl C, Krarup AL, Olesen AE, Brock C, Graversen C, Drewes AM. Is electrical brain activity a reliable biomarker for opioid analgesia in the gut? Basic & clinical pharmacology & toxicology. 2011;109:321–327. [PubMed]
  • Sumpter AL, Holford NH. Predicting weight using postmenstrual age–neonates to adults. Paediatric anaesthesia. 2011;21:309–315. [PubMed]
  • Sverrisdottir E, Foster DJ, Upton RN, Olesen AE, Lund TM, Gabel-Jensen C, Drewes AM, Christrup LL, Kreilgaard M. Modelling concentration-analgesia relationships for morphine to evaluate experimental pain models. European journal of pharmaceutical sciences : official journal of the European Federation for Pharmaceutical Sciences. 2014;66c:50–58. [PubMed]
  • Højsgaard Søren, Halekoh Ulrich, Robison-Cox Jim, Wright Kevin, Leidi AA. doBy - Groupwise summary statistics, LSmeans, general linear contrasts, various utilities. 2014 CRAN.R-project.org.
  • Tibboel D, Anand KJ, van den Anker JN. The pharmacological treatment of neonatal pain. Seminars in fetal & neonatal medicine. 2005;10:195–205. [PubMed]
  • van den Anker JN. Pharmacokinetics and renal function in preterm infants. Acta paediatrica (Oslo, Norway : 1992) 1996;85:1393–1399. [PubMed]
  • Walker SM. Biological and neurodevelopmental implications of neonatal pain. Clinics in perinatology. 2013;40:471–491. [PubMed]
  • Wang C, Sadhavisvam S, Krekels EH, Dahan A, Tibboel D, Danhof M, Vinks AA, Knibbe CA. Developmental changes in morphine clearance across the entire paediatric age range are best described by a bodyweight-dependent exponent model. Clinical drug investigation. 2013;33:523–534. [PubMed]
  • Wickham H. Reshaping Data with the {reshape} Package. Journal of Statistical Software. 2007;21:1–20.
  • Wickham H. ggplot2: elegant graphics for data analysis. Springer; New York: 2009.
  • Wickham H. plyr - The Split-Apply-Combine Strategy for Data Analysis. Journal of Statistical Software. 2011;40:1–29.
  • Wojciechowski J, Hopkins AM, Upton RN. Interactive Pharmacometric Applications Using R and the Shiny Package. CPT: pharmacometrics & systems pharmacology. 2015;4 [PMC free article] [PubMed]
  • Zhang L, Beal SL, Sheiner LB. Simultaneous vs. sequential analysis for population PK/PD data I: best-case performance. Journal of pharmacokinetics and pharmacodynamics. 2003;30:387–404. [PubMed]