Data from a prospective clinical study with 1000 patients undergoing ambulatory surgery were analysed. The study was registered to EudraCT (no. 2015-005561-23) and approved by the Coordinating Ethics Committee of the Helsinki and Uusimaa Hospital District (3/2016) and the Finnish National Agency for Medicines (KL no. 4/2016). The data were collected at the day-surgery unit of Jorvi Hospital, University of Helsinki, and Helsinki University Hospital, Espoo, Finland from August 2016 to March 2018. When the patients were recruited to the study, medication used at home before surgery was recorded. Use of strong inhibitors of CYP2D6 was an exclusion criterion. All study patients received oral premedication with paracetamol (1000 mg) and codeine (60 mg) for preoperative prophylactic analgesia. After pre-medication, two blood samples were obtained through a cannula placed in the forearm vein. The first blood sample was drawn after 20–60 min and the second sample was drawn after 180–360 min after receiving the paracetamol-codeine combination tablets. The time of blood sampling was adapted to the progress of the treatment process, and the blood samples were commonly taken during the induction of anaesthesia and at the time of discharge from hospital.

The adverse effects of pain medication, such as nausea and vomiting (yes or no), dizziness, sleepiness and constipation (score from 0 to 3) were asked repeatedly in the recovery room after surgery and after discharge from hospital twice daily until 2 days after the surgery. Use of antiemetics in the recovery room was also recorded. Detailed information about the study protocol and conduct has been provided in Online Resource 1.

Plasma concentrations of codeine, morphine and the glucuronide conjugates of codeine and morphine were determined with high-performance liquid chromatography-tandem mass spectrometry. Morphine, codeine, morphine 6-glucuronide, codeine-6-glucuronide and stable labelled morphine-D6 and codeine-D6 were purchased from Toronto Research Chemicals (Toronto, Canada). Prior to analysis, the plasma samples were pre-treated using a Strata X-C solid phase extraction in a 96-well format (Phenomenex, Torrance, CA, USA). Briefly, the samples (0.1 mL) were diluted with 0.15 mL of 4% phosphoric acid containing the internal standards and loaded into the pre-conditioned extraction wells. The wells were then washed with 0.15 mL of 2% formic acid, followed by 0.15 mL of methanol, and the analytes were eluted using 0.15 mL of 5% ammonium hydroxide in methanol. Finally, the sample extracts were dried using a centrifugal evaporator (GeneVac, Thermo Fisher Scientific, USA) and reconstituted in 0.1 mL of 5% methanol.

All measurements were carried out using a Sciex 6500 Qtrap LC-MS system interfaced with an electrospray ion source (Sciex, Toronto, Ontario, Canada). The analytes were separated on a Kinetex biphenyl column (2.6 µm particle size, 2.1 × 100 mm internal diameter; Phenomenex, Torrance, CA, USA) under gradient conditions and a 0.28 mL/min flow rate. The mobile phase consisted of a mixture of 0.1% acetic acid (solvent A) and acetonitrile-methanol (20:80 v/v) (solvent B), and the solvent gradient was set as follows: a linear ramp from 5 to 30% B over 3.5 min, a second linear ramp from 70 to 90% B over 1.5 min and 3 min equilibration back to the initial eluent composition. The mass spectrometer was operated in positive polarity mode, and the mass-to-charge (m/z) transitions used for quantification were 286–152 for morphine, 462–286 for morphine 6-glucuronide, 300–152 for codeine and 476–300 for codeine-6-glucuronide. The limit of quantification was 0.05 ng/mL, except for morphine-6-glucuronide (0.1 ng/mL). The day-to-day coefficients of variation (CV) values were below 11% at relevant concentrations for all analytes.

Genotyping was carried out using TaqManR genotyping and copy number assays on the QuantStudio™ 12K Flex Real-Time PCR System or targeted next-generation sequencing using the Ion GeneStudio™ S5 Prime system (Thermo Fisher Scientific, Waltham, MA). All samples were genotyped for the clinically relevant CYP2D6 sequence variations, defining the *1-, *2-, *3-, *4-, *5-, *6-, *9-, *10-, *15-, *17-, *35-, *39- and *41-alleles, as well as for CYP2D6 copy number variation (CNV). The CYP2D6-genotype distribution in the study population is shown in additional data provided in Online Resource 2. CYP2D6 phenotypes were inferred from the genotypes using the activity scores on the basis of recent consensus [5, 6, 12]. A detailed description of the genotyping is found in the Online Resource 1.

Non-linear mixed effects modelling was conducted with NONMEM® (version 7.4.3 or above) (ICON Development Solutions, Ellicott city, MD, USA) [17], assisted by the Perl-Speaks-NONMEM (PsN)-toolkit [18].

The pharmacokinetic modelling of concentration–time data was conducted in a stepwise manner. First, a pharmacokinetic model was developed and validated for codeine. Next, we included morphine data in the dataset and tested parent-metabolite disposition models. The C6G and M3G data were then added in the analysis data in a sequential manner to improve model performance in accounting for the complex enzymatic interplay defining codeine and morphine elimination (Fig. 1). Finally, the influence of CYP2D6 activity scores (AS) was added to the model, and several parameterisations were tested. The influence of the remaining observable patient characteristics (covariates) on codeine/morphine metabolism and disposition profiles was tested using the stepwise-covariate-modelling tool in the PsN toolkit [18]. Detailed description of the modelling process and the differential equations derived from the final model are given in the Online Resource 1.

Schematic of the final structural model for codeine and its metabolite morphine. The pharmacokinetics of codeine and morphine are described with one-compartmental mammillary models with the first-order absorption of codeine. The total clearance of codeine (CLcod) is divided into CYP2D6- (ke × fmor) and non-CYP2D6-mediated [ke × (1 – fmor)] clearance. CYP2D6-mediated clearance is used as an input in the C6G model, of which 60% is subsequently metabolised to M3G (ke,mor × fM3G) and 40% is metabolised via other pathways [ke,mor × (1 − fM3G )]. I(t), oral dose; F, oral bioavailability; fmor, ratio metabolised to morphine; ke, elimination rate constant for codeine; MPR, metabolite-to-parent molecular weight ratio; ke,mor, elimination rate constant for morphine; fM3G, fraction of morphine metabolised to M3G; ke,C6G, elimination rate constant for C6G; ke,M3G, elimination rate constant for M3G

Models were coded using differential equations and solved by the ADVAN13 subroutine of the NONMEM® software. The model parameters were modelled in the log-domain using MU-referencing and computed using the first-order conditional estimation with interaction estimation method.

Exponential model was used to characterise between-subject variability in model parameters, and the covariance between the random effect parameters was also tested. The residual variability in the model predictions was specified using an additive error model. Both were assumed to be normally distributed random variables with a mean of zero and a SD of one.

Improvement in model fit was measured on the basis of significant differences in objective function value (OFV) and Akaike information criteria (AIC), the numerical stability of the model parameters, standard goodness-of-fit plots and visual predictive checks [19].

The pharmacokinetic parameters defining clearance and the volume of the central compartment for both codeine and morphine were allometrically scaled to the population’s median weight. The covariate effects of CYP2D6 activity on the ratio of codeine metabolised to morphine were tested sequentially using CYP2D6 phenotypes predicted from genotypes or activity scores to explain between-subject variability in codeine metabolism. Both categorical and continuous modelling approaches were tested, as described in Online Resource 1.

The stepwise covariate modelling protocol included in the PsN toolkit was used to test the influence of the additional patient covariates (age, cigarette smoking, and American Society of Anesthesiologists classification) on the fixed effects. For continuous covariates, power and exponential models were tested, as was a linear model specification for categorical covariates. Both forward and backward searches were performed, with search probabilities of 0.05 and 0.01, respectively.

Model uncertainty was evaluated with the sampling importance resampling protocol [20] to determine the robustness of the final pharmacokinetic models. The protocol was run with 1000 samples and 500 resamples (M/m = 2). Secondly, the final pharmacokinetic model was used to simulate data for 1000 virtual study subjects who had designated activity scores and were receiving a clinically plausible range of codeine doses (from 30 mg o.d. for 1 day up to 60 mg q.i.d. for 4 days), using a typical population median body weight of 80 kg. Additionally, a new patient collective of 10,000 new patients was simulated using the covariate set of study subjects at codeine dose levels of 30 mg or 60 mg (t.i.d. or q.i.d.) to evaluate the effect of CYP2D6 phenotype on codeine/morphine disposition kinetics. Each phenotype class was simulated against varying codeine dosages. The non-compartmental analysis protocol in the PsN toolkit and ncappc R package was used to compute non-compartmental pharmacokinetic analysis.