Pharmacokinetic modeling of R and S-Methadone and their metabolites to study the effects of various covariates in post-operative children.

Methadone is a synthetic opioid used as an analgesic and for the treatment of opioid abuse disorder. The analgesic dose in the pediatric population is not well-defined. The pharmacokinetics (PKs) of methadone is highly variable due to the variability in alpha-1 acid glycoprotein (AAG) and genotypic differences in drug-metabolizing enzymes. Additionally, the R and S enantiomers of methadone have unique PK and pharmacodynamic properties. This study aims to describe the PKs of R and S methadone and its metabolite 2-ethylidene-1,5-dimethyl-3,3-diphenylpyrrolidine (EDDP) in pediatric surgical patients and to identify sources of inter- and intra-individual variability. Children aged 8-17.9 years undergoing orthopedic surgeries received intravenous methadone 0.1 mg/kg intra-operatively followed by oral methadone 0.1 mg/kg postoperatively every 12 h. Pharmacokinetics of R and S methadone and EDDP were determined using liquid chromatography tandem mass spectrometry assays and the data were modeled using nonlinear mixed-effects modeling in NONMEM. R and S methadone PKs were well-described by two-compartment disposition models with first-order absorption and elimination. EDDP metabolites were described by one compartment disposition models with first order elimination. Clearance of both R and S methadone were allometrically scaled by bodyweight. CYP2B6 phenotype was a determinant of the clearance of both the enantiomers in an additive gene model. The intronic CYP3A4 single-nucleotide polymorphism (SNP) rs2246709 was associated with decreased clearance of R and S methadone. Concentrations of AAG and the SNP of AAG rs17650 independently increased the volume of distribution of both the enantiomers. The knowledge of these important covariates will aid in the optimal dosing of methadone in children.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

R and S methadone have varied pharmacodynamic effects and undergo differential CYP metabolism to 2‐ethylidene‐1,5‐dimethyl‐3,3‐diphenylpyrrolidine (EDDP). Pharmacogenomic variation in drug metabolism and concentrations of the binding protein AAG contributes to the interindividual variability in methadone drug.

WHAT QUESTION DID THIS STUDY ADDRESS?

This study aimed to identify significant sources of variability in the pharmacokinetics (PKs) of R and S methadone, including physiological parameters and the relevant genotypes in children undergoing major surgeries.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

The concentration of AAG and its ORM1 genotype independently affect the volume of distribution of both methadone enantiomers. Novel genotypes in CYP2B6, in addition to the haplotype‐based activity score, were found to contribute to variability in the clearance of R and S methadone.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

Knowledge of the effects of these genotypes on methadone’s PKs could help determine the precise dose of methadone for surgical pain relief. In addition, this will promote research into these intronic single‐nucleotide polymorphisms on the metabolism of various drugs.

INTRODUCTION

Postoperative care demands potent analgesia with opioids for the adequate control of acute surgical pain. Despite the use of around‐the‐clock analgesics along with rescue medication, about 40% of children still experience moderate to severe pain. Methadone is a synthetic opioid with a long half‐life, which is indicated for analgesia in addition to its use in opioid use disorder. Methadone‐based multimodal analgesia has been shown to provide best‐in‐class analgesia in acute surgical pain. Unfortunately, methadone possesses wide variability in its exposure in children. Although this variability has largely been attributed to the binding protein AAG (alpha 1 acid glycoprotein) and the metabolizing CYP enzymes, the knowledge of the effect of these and other covariates have not been well‐studied. ,

Administered as a racemic mixture, the R and the S enantiomers of methadone differ in their pharmacodynamic and pharmacokinetic (PK) properties. The R enantiomer is mainly responsible for analgesia, sedation, and respiratory depression, whereas the S enantiomer is responsible for QT prolongation. Methadone is stereoselectively oxidized by various cytochrome P450 (CYP) enzymes to the primary metabolite 2‐ethylidene‐1,5‐dimethyl‐3,3‐diphenylpyrrolidine (EDDP), which is largely inactive. R and S methadone are metabolized primarily by CYP3A (67.9% and 70.2%, respectively), CYP2B6 (14.3% and 23.0%, respectively), and CYP2C19 (13.0% and 3.9%, respectively), with CYP2C9 and CYP2D6 contributing to a lesser extent. , ABCB1 (P‐glycoprotein transporter) and POR (cytochrome P450 reductase) pharmacogenetics may also contribute to variability in methadone PKs. ,

Methadone is 80–90% bound to the acute‐phase protein AAG in plasma. Increased concentrations of AAG during stressful conditions like inflammation, infections, burns, and surgeries can affect the effective exposure of methadone by reducing its unbound concentrations. The orosomucoid genes ORM1 and ORM2 encode for AAG. Although there are some data available on the binding of methadone to AAG, the effect of its genetic polymorphism on the binding of the drug is unknown. The baseline concentrations of AAG in a pediatric population and the changes in a peri‐operative setting are largely unknown. This problem is compounded in the pediatric population by the paucity of available PK data for the drug.

The objective of this study is to evaluate the variability in the PKs of methadone in a postoperative pediatric population and study the effect of various covariates on the PKs.

MATERIALS AND METHODS

Study population

This study was approved by the Indiana University institutional review board (IRB #1707525204) and registered at clinicaltrials.gov (NCT03495388). Children and adolescents undergoing pectus excavatum or posterior spinal fusion surgeries aged 8–17.9 years were included in the study after informed consent from a parent or a legal guardian and written assent from the child when applicable. Children were excluded from the study if allergic to methadone, diagnosed with developmental delay, neurologic disorder, liver or renal disease, or had preoperative pain requiring analgesics.

Study design

The first dose of racemic methadone was administered intra‐operatively by intravenous infusion and every 12 h postoperatively as an oral dose for a total of four to six doses at 0.1 mg/kg. Oral doses were administered either as a tablet or a suspension depending on the patient’s convenience. Blood samples (5–8) were collected at specific time windows for each of the individuals across the first three interdose intervals.

Analytic methods

R and S methadone and R‐ and S‐EDDP was quantified from plasma by the IU Clinical Pharmacology Analytical Core Laboratory using liquid chromatography tandem mass spectrometry. The standard curve was linear from 0.015–150 ng/ml. Coefficient of variance were less than 15% for all the quality controls for all these analytes. AAG was quantified from plasma by high‐performance liquid chromatography‐UV as previously described. The standard curve was linear from 20–1500 µg/ml and inter‐ and intra‐day coefficients of variability were less than 10%. Additional details are provided in Supplementary Methods.

Genotyping

Whole‐genome sequencing libraries were prepared with the Illumina DNA PCR‐Free Library Prep Kit according to the manufacturer’s instructions. The prepared libraries were sequenced and the high‐quality reads were aligned to human reference genome hg19 using BWA‐MEM (version 0.7.15). , Sentieon version 201911.01 (Sentieon, Inc., https://www.sentieon.com/) was applied for variants detection. Subjects’ whole‐genome sequencing BAM and corresponding VCF files were uploaded into LifeOmic’s (Indianapolis, IN, USA) Precision Health Cloud (PHC) for PGx analysis. Aldy version 3.0 was operationalized in the PHC to genotype pharmacogenomic star alleles and diplotypes for CYP2B6, CYP2C8, CYP2C19, CYP2D6, CYP3A4, and CYP3A5. In the PHC, the JupyterLab notebook was used to extract genotypes of interest not covered by Aldy from the VCF files. Additional details are provided in Supplementary Methods.

Pharmacokinetic modeling

R and S methadone plasma concentration‐time profiles were modeled by nonlinear mixed‐effects modeling using NONMEM 7.4 using first‐order conditional estimation method with η–ε interaction. The minimum objective function value (OFV; −2 log‐likelihood) calculated by NONMEM was used to discriminate between nested models, with a change in OFV (dOFV) of 3.84 considered significant (χ2, one degree of freedom, p < 0.05). Model assessment was based on successful minimization, OFV, successful run of the covariance step, relative standard error of the parameter estimates, goodness of fit plots, and eta‐shrinkage. The predictive performance of the models was evaluated using prediction corrected visual predictive checks (n = 1000). The robustness of the final parameter estimates and the stability of the final model were validated using bootstrapping by generating 1000 resampled datasets with replacement.

Models were developed using a step‐wise process. The PK model for the drug was first built and parent drug parameter estimates fixed prior to fitting the metabolite model. The final model structure is described in Figure 1.

FIGURE 1

Structural pharmacokinetic model for both R and S methadone. Oral drug was dosed into the gastrointestinal dosing compartment. F1 is bioavailability of oral dose. V2 is the central compartment and V4 the peripheral compartment of methadone. The metabolite EDDP was fit to a single compartment, V3, which was parameterized as V2*Vf, where Vf is a scaling factor between the central volume of distribution of the drug and the metabolite. Vf was fixed to 1 in the model due to the unidentifiability of the metabolite parameters. Q represents the intercompartmental clearance. CL is the total clearance of methadone and CLF is the fraction of CL to EDDP (CL2). CL1 represents the clearance of the drug by all other routes and CL3 was the clearance of the metabolite

Covariate models

The covariates analyzed included biological sex, body weight, fat‐free mass, body mass index, race, AAG concentrations, and genotypic variants in CYP enzymes (CYP3A4, CYP3A5, CYP2B6, CYP2C19, CYP2C8, CYP2D6, and CYP2C18), P‐glycoprotein, AAG, pregnane X receptor (PXR), and constitutive androstane receptor (CAR). Covariate selection was performed using a step‐wise covariate approach with a forward selection (p = 0.05) and a more stringent backward elimination step (p = 0.01). AAG was measured at every time point the methadone concentrations were measured and were included as a time‐varying covariate in the model. Genotypic variants in CYP enzymes where allelic phenotype information was available were classified as their corresponding phenotypes – ultra‐rapid metabolizer, rapid metabolizer, normal metabolizer, intermediate metabolizer, and poor metabolizer based on available literature and were given a corresponding activity score (2, 1.5, 1, 0.5, or 0). The relationships between the PK parameters and the activity score of the CYP enzymes were checked using a linear, power, and exponential relationship. Single‐nucleotide polymorphisms (SNPs) without phenotype information were evaluated using dominant, additive, and recessive models. Additional details are provided in the Supplementary Methods.

RESULTS

Sixty‐one children and adolescents aged 11–17 years participated in the study. Demographic details are given in Table 1. Some methadone (0.46%) and metabolite concentrations (9.4%) were below the lower limit of quantification, leaving 430 drug and 381 metabolite concentrations available for analysis. Nine children had missing genotype information.

TABLE 1

Demographic information in the study population

Demographic	Median (IQR) N = 61
Dose of methadone (mg/kg body weight)	0.087 mg/kg (IQR 0.069–0.094)
Age (years)	14.74 (13.62–15.66)
Weight (kg)	53.60 (47.90–60.10)
Height (cm)	164.50 (158.00–171.50)
BMI	19.40 (17.61–22.50)
Baseline AAG (ng/ml)	84 (62–109)
Biological sex	Females 31, Males 30
Race	White 49, African American 7, Hispanics 3, Native American 1, unknown 1
Type of surgery	Pectus excavatum repair = 25 Idiopathic scoliosis spinal fusion surgery = 36

Abbreviations: BMI, body mass index; IQR, interquartile range.

Structural model of R methadone

The R methadone concentration‐time data were best described by a two‐compartment disposition model (Figure 1) with interindividual variability on the total clearance of the drug (CL), the volume of distribution of the central compartment (V2), and the peripheral compartment volume (V4). The absorption phase was well‐described by a first‐order process estimated with unique absorption rate constants (Katablet and Kasusp) for the tablet and the suspension formulations. This resulted in a better fit than using a single absorption parameter (dOFV = −8.3). Oral bioavailability was estimated without the estimation of interindividual variability. Individual bioavailability of the tablet and the suspension formulations could not be estimated separately and therefore were assumed to be equivalent.

The metabolite model was subsequently added to the drug model after fixing the drug parameters to their model estimates. The volume of distribution of the metabolite was nonidentifiable and therefore was set equal to the central volume of distribution of the drug. The metabolite concentrations were well‐described by a one‐compartment model with interindividual variability on clearance of the metabolite (CL3). The η and the ε shrinkage were too low to graphically explore the covariate relationships.

Covariate modeling on R methadone

The covariate relationships evaluated for R methadone are shown in Table 2. Bodyweight was added allometrically to the clearance of the drug with a fixed exponent of 0.75. Allometric scaling of the volume parameters with bodyweight worsened the model fit as evidenced by an increase in the OFV (dOFV = 1.2) and a lack of correlation in the diagnostic plots. Therefore, only clearance was allometrically scaled. Estimation of scaling factor of clearance improved the fit with an estimated exponent of 1.2, but the factor was fixed to a standard value of 0.75 to prevent false estimation or bias in the estimation of scaling factor.

TABLE 2

Final population pharmacokinetic and pharmacodynamic parameters for R and S methadone with the parameter estimates obtained from bootstrapping (N = 1000)

Parameter	Population estimates (%RSE)	95% CI	%CV for BSV (%RSE)	95% CI	Shrinkage (%)
R methadone
Kasusp (h−1)	0.318 (28.3)	0.151–0.499
Katablet (h−1)	0.123 (47.2)	0.0706–0.351
CL (L*h−1)	15.7 (27.1)	7.58–24.3	72.1 (43.9)	17.2–129	22.5
V2 (L)	176 (16.6)	113–225	79.4 (19.9)	52.6–117	20.3
Q (L*h−1)	69.2 (39.3)	49.7–162
V3 (L)	335 (23.6)	212–511	62.23 (36.4)	13.8–116	41.6
F	0.718 (13.2)	0.537–0.910
V2AAG	−0.00291 (31.3)	−0.00323 to −0.000428
V2rs17650	−0.443 (35.4)	−0.700 to −0.119
CLF	0.217 (31.9)	0.133–0.429	65.0 (24.6)	26.7–96.8	37.9
VF	1 FIX
CL3 (L*h−1)	25.7 (37.3)	14.5–56.3	49.8 (27.2)	18.3–79.8	37
CLFCYP2B6	0.745 (17.4)	0.370 −0.863
CLFrs2246709	0.450 (33.9)	0.157–0.708
RUVdrug	0.165 (5.17)	0.132–0.201			8.4
RUVmetabolite	0.207 (5.38)	0.167–0.252			9.1
S methadone
Kasusp (h−1)	0.432 (22.8)	0.201–00601
Katablet (h−1)	0.257 (42.4)	0.122–0.579
CL (L*h−1)	13.0 (16.7)	9.35–17.9	40.9 (23.2)		31.7
V2 (L)	98.3 (12.9)	75.4–126	63.6 (25.3)		26.6
Q (L*h−1)	105 (18.9)	45.2–139
V3 (L)	139 (19.3)	90.4–196	116 (23.6)		15.4
F	0.606 (14.8)	0.468–0.808
V2AAG	−0.00192 (51.8)	−0.00317 to −0.000116
V2rs17650	0.526 (28.4)	−0.709 to −0.173
CLF	0.135 (23.7)	0.0924−0.202	47.9 (24.6)	15.5–67.8	26.9
VF	1 FIX
CL3 (L*h−1)	7.97 (23.7)	5.37–13.0	33.7 (43.3)	8.01–70.1	57.3
CLFCYP2B6	0.636 (25.2)	0.198–0.771
CLFrs2246709	1.68 (57.8)	0.198–4.27
RUVdrug	0.165 (6.15)	0.130–0.212			9.2
RUVmetabolite	0.194 (5.47)	0.157–0.239			8.6

Kasusp and Katablet were the first‐order absorption rate constant of the oral formulations suspension and tablet, CL was the total clearance of the drug, V2 was the central volume of distribution of the drug, Q was the intercompartmental clearance of the drug, V3 was the central volume of distribution of the metabolite, F was the common bioavailability term for both the oral formulations, V2AAG and V2rs17650 denote the covariate relationship of V2 with AAG concentrations and genotype of ORM1 rs17650, respectively, as described in the Results section. CLF was the fraction of CL that contributes to CL2, the clearance toward the formation of the metabolite. VF was the fixed factor by which V2 and V3 were related. CL3 was the clearance of the metabolite. CLFCYP2B6 and CLFrs11884424 denote the covariate relationship of CLF with the activity score of CYP2B6 and intronic SNP on CYP3A4 rs2246709, respectively. RUVdrug and RUVmetabolite were the residual unexplained variability of the drug and the metabolite respectively.

Abbreviations: BSV, between subject variability; CI, confidence interval; CL, clearance; CV, coefficient of variation; RSE, relative standard error.

AAG concentration was a significant covariate on the volume of distribution of the central compartment V2 using method one as described in methods (dOFV = −7.1; Figure S1). AAG was also analyzed using method two. However, this was not maintained as the RSE of the estimated ΘDAAG and improvement compared to method one was marginal (dOFV = −3.3). After addition of the effect of AAG concentration, an SNP of ORM1 (gene encoding for AAG) rs17650 was also found to affect the V2 in an additive gene model (dOFV = −6.5; Figure S2). This effect is independent of the concentrations of AAG (Figure S1).

The CYP2B6 activity score exhibited a linear relationship with the fractional clearance (CLF) of R methadone to R‐EDDP (dOFV = −8.8), as shown in Figure 2. The addition of the CYP3A4 SNP rs2246709 as a covariate on CLF further improved the model (dOFV = −7.1; Figure S3). None of the other relationships which were found to be significant by univariate analysis (Table S2) significantly reduced the OFV when incorporated into the PK model (Table 3).

FIGURE 2

Relationship between CYP2B6 phenotype and CLF of R and S methadone. The CYP2B6 activity score was associated with the clearance of R and S methadone to their primary metabolites R and S‐EDDP (p = 0.01, R 2 = 0.11 for R methadone, and p = 0.01, R 2 = 0.10 for S methadone)

TABLE 3

List of covariates which were positive in base model (p < 0.05) but were insignificant (p > 0.01) with addition of other covariates

R methadone
SNP	Gene	PK parameter	Gene model	dOFV from final model
rs7250601	CYP2B6	CLF	Recessive	−1.375
rs4646437	CYP3A4	CL	Recessive	−3.635
rs1042194	CYP2C18	CL3	Additive	−4.285
rs2229109	ABCB1	CL3	Recessive	−10.2
Activity score of CYP2C19	CYP2C19	CL3	Additive	−6.279

Each of the single nucleotide polymorphism (SNP) belong to the corresponding gene in the adjacent column. Each CYP enzyme has its corresponding gene named according to the enzyme name. Other SNPs, ABCB1‐p glycoprotein; CAR‐constitutive androstane receptor (CAR) also known as nuclear receptor subfamily 1, group I, member 3 is responsible for the regulation of multiple enzymes including members of the CYP2B and CYP3A subfamilies; rs17180299 is an intergenic SNP earlier shown to be associated with the plasma concentration of R methadone.

Abbreviations: dOFV, objective function value difference; PK, pharmacokinetic.

The final model for R methadone has the following covariate relationships. where CYP indicates the CYP2B6 activity score and rs2246709 indicates the number of active alleles rs2246709. V2 = 176 L * (1–0.443 * (number of active alleles rs17650–1)) * (1–0.00291 * (AAG–94.76)).

Where rs17650 is the number of active alleles of rs17650 and AAG is the plasma concentration of AAG.

The CLF for a 70 kg individual with CYP2B6 activity score of one and a heterozygous rs2246709 is 0.217 L/h. The V2 in an individual with a heterozygous rs17650 and an average AAG concentration of 94.76 µg/ml is 176 L.

A total of nine (14.8%) out of the 61 individuals did not have any genotype information. An estimation method using the mixture modeling for the identification of the missing genotype information was undertaken. However, this method was not successful in identifying the mixture population (probability of mixture estimate remained above 0.99 for any one of the populations) for most of the covariates. Therefore, a simple imputation method with mode values was done for all of the missing genotypes. The goodness of fit plots are shown in Figure S4.

Structural model of S methadone

The S methadone concentration‐time data were best described by a two‐compartment disposition model (Figure 1) with interindividual variability on CL, volume of the central compartment V2, and volume of the peripheral compartment V4. The absorption phase was well‐described by a first‐order process estimated with unique absorption rate constants (Katablet and Kasusp) for the tablet and the suspension formulations (dOFV = −3.9, compared to a single ka). As with the model for R methadone absolute oral bioavailability was estimated without the estimation of its interindividual variability and separate bioavailability of the tablet and the suspension could not be estimated.

The S‐EDDP metabolite model was developed after establishing the model for S methadone and fixing parameters to the final estimates. The volume of distribution of the metabolite was not identifiable and therefore was set equal to the central volume of distribution of the drug. The metabolite concentrations were well‐described by a one‐compartment model. Interindividual variability was added for the clearance of the metabolite (CL3). The η and the ε shrinkage were too low to graphically explore the covariate relationships.

Covariate modeling on S methadone

The covariates relationships evaluated for S methadone are shown in Table 2. Body weight was added allometrically to the clearance of the drug with a fixed exponent of 0.75 despite the improvement with the estimated scaling factor of 0.67. AAG was a significant covariate on V2 (dOFV = −6.3). Additionally, the ORM1 SNP rs17650 was significantly associated with V2 and incorporation into the PK model led to improved model fit (dOFV = −7.2).

As with R methadone, the CYP2B6 activity score was found to have a linear relationship with S methadone CLF (dOFV = −9.1). The intronic SNP rs1188424 further improved the model fit when used as a covariate on CLF (dOFV = −13.3). The CYP3A4 SNP rs2246709 was found to be significantly (dOFV = −9.0) related to CLF in the absence of rs11882424 in the model. The genotypic covariates whose relationships were tested in the model but were not included to the final model are given in Table 3.

The final model for S methadone has the following covariate relationships. where CYP indicates CYP2B6 activity score and rs11882424 is 1 for homozygous recessive rs11882424. With this model, the CLF in a 70 kg individual with CYP2B6 activity score of one and dominant rs11882424 genotype is 0.135 L/h. Another model was also considered in parallel because of the questionable plausibility of rs11882424 effects on clearance: where CYP indicates the CYP2B6 activity score and rs2246709 is the number of active alleles rs2246709. In this model, the CLF in a 70 kg individual with a CYP2B6 activity score of one and a heterozygous rs2246709 is 0.128 L/h. where rs17650 indicates the number of active alleles rs17650 and AAG is the plasma concentration of AAG. The V2 in an individual with a heterozygous rs17650 and an average AAG concentration of 94.76 µg/ml is 98.3 L.

This second model was used as the final model for the subsequent analyses. An estimation method using the mixture modeling for the identification of the missing genotype information was performed but was not successful in estimating the subpopulation. Therefore, a simple imputation method with mode values was done for all of the missing genotypes. The goodness of fit plots are shown in Figure S4.

Model qualification

The visual predictive checks indicate that the simulated concentrations of R and S methadone and their respective metabolites agree with the observed data. The 5th, 50th, and 95th percentiles of the observed data fall within the respective 95% confidence intervals of the simulated data (Figure 3). A bootstrapping procedure (n = 1000) found the parameter estimates to be reliable with acceptable relative standard errors (Table 2).

FIGURE 3

Prediction‐corrected visual predictive checks of R and S methadone and their corresponding metabolites. Grey circles represent the observed concentrations, and the solid red line represents 50th percentile, whereas the dotted lines represent the 5th and 95th percentiles of the observed data. The shaded areas represent the 95% confidence intervals around the simulated 5th, 50th, and 95th percentiles

DISCUSSION

In this work, the plasma concentration of the enantiomers of methadone and its metabolite EDDP were modeled and several potential covariates on their PK parameters were explored in a pediatric population undergoing major surgeries. We found that the CYP2B6 activity score had a significant effect on CL and AAG concentrations had a significant effect on the volume of distribution of the drug. Interestingly, the volume of distribution of both R and S methadone was not only related to AAG concentrations but also to ORM1 rs17650 genotype. This may be due to differences in binding affinity associated with the rs17650 variant. This is a novel finding concerning the PKs of methadone. We also report for the first time that rs2246709, an intronic SNP of CYP3A4, was associated with reduced clearance of both R and S methadone. The effects of these covariates on the PKs of R and S methadone may in part explain variability in analgesic effect observed in pediatric surgery patients.

The study data consisted of both intravenous and oral doses, either tablet or suspension, and therefore enabled the estimation of absolute bioavailability and the accurate estimation of the PK parameters. Although unique absorption rate constants could be estimated for the two oral formulations, individual bioavailability terms could not be determined. In addition, the introduction of absorption lag time or a transit compartment absorption model did not improve the fit. This was probably due to the sparsity of concentrations in the absorption phase and the use of two different oral formulations in our population. Other studies which have estimated clearance in a two‐compartment model have found it to range from 7 to 9 L/h, , , which was much lower than the estimates in this study (13–15 L/h). This difference could be attributed to the pediatric population in this study and the fluctuating hemodynamics and metabolic changes seen in the immediate postoperative period.

Allometric weight scaling was done on the clearance of both R and S methadone. Estimation of the allometric scaling factor improved the fit, but the factor was fixed to a standard value of 0.75 to prevent false estimation or bias in the estimation of scaling factor with a smaller sample size as shown in a recent publication. Inclusion of a maturation model led to increased model instability with both the enantiomers. As we recruited primarily adolescents (age range 10.8 to 17.9), there is likely little effect of maturation in our population.

The below the limit of quantification (BLOQ) concentrations were predominantly samples collected toward the end of the dosing period in the case of the parent drug, and at early timepoints for the metabolites. The BLOQ concentrations were omitted according to the M1 method due to the relatively small numbers of BLOQ data.

The concentrations of the binding protein AAG have been reported to increase up to two to seven‐fold with physiological stress. AAG concentrations may also be impacted by diurnal variations and sex. In a peri‐operative study in infants, AAG concentrations increased postoperatively, peaking at 3–4 days and declining to the baseline in 2–4 weeks. In this study, AAG increased in most patients in the postsurgical period with a median baseline concentration of 84 µg/ml and range of 24.48 to 205.38 µg/ml across the study. This variability increased through the period of in‐hospital stay in most patients, consistent with previous studies (Figure S5). , The median concentration of AAG at the end of the sample collection was 115.75 µg/ml (range 30.2–354.9 over 14.5–53.4 h postsurgery). The effect of AAG concentration was captured in the time‐varying covariate relationship with the volume of distribution of the drug (V2). AAG is encoded by two loci ORM1 and ORM2. ORM2 is widely preserved and is not polymorphic, whereas ORM1 exhibits polymorphisms. ORM1 genotype has previously been reported to be an important predictor of dose required for drugs like warfarin and telmisartan. , Rs17650 determines the phenotype of AAG as F and S and rs1126801 further classifies F as F1 and F2. The variant F phenotypes have lesser binding for methadone than the S phenotype. An interesting finding in this work was the relationship between rs17650 and V2 in the study, even after accounting for AAG concentration. The addition of rs17650 genotype to V2 of methadone resulted in a significant improvement in model fit. Rs1126801 was poorly represented in our study population (n = 2) and therefore was not evaluated in the model.

CYP enzymes, including CYP3A4, CYP2B6, CYP2C19, CYP2C8, CYP2D6, and CYP2C18 are involved in the metabolism of methadone in vitro. , CYP2B6 has been consistently shown to contribute to the variability in the metabolism of methadone. , , , Due to the vast number of CYP enzymes involved in the metabolism of methadone, it is difficult to isolate the effect of an individual enzyme on clearance in vivo. We used whole genome sequencing to determine CYP variants. Therefore, in addition to common * haplotypes, we identified SNPs that have not been integrated into haplotype nomenclature. , A two‐staged process was used in the selection of the covariates, a univariate regression analysis followed by step‐wise covariate model building. CYP2B6, CYP2D6, CYP2C19, CYP3A4, and CYP3A5 were tested using their activity scores as these alleles are well studied and the phenotypic classifications or activity scores are reasonably known for multiple drugs from published literature. , , , Most of these were not found to be related to the clearances of methadone, contrary prior in vitro and in vivo studies. The CYP2B6 activity score was found to be strongly related to the metabolism of the parent drug to the primary metabolite EDDP for both the enantiomers. The rs11882424, an intronic SNP of CYP2B6, was previously found to affect the clearance of methadone, although the direction of influence is unclear. We found expression of the variant allele at rs1182424 to significantly increase (dOFV = −13.2) the CLF of S methadone. The variant allele of the same was found to increase the CLF of S methadone, which has not been reported in literature elsewhere in the past. As rs1182424 is not commonly genotyped, we included the more commonly assessed rs2246709 (dOFV = −9.0) in our final model.

CYP2C19 activity scores were associated with the clearance of R and S‐EDDP but were not retained in the final model because of the use of a stricter backward elimination step and a marginal increase in RSE of the clearance of the metabolite CL3. Wang et al. had demonstrated the effect of the CYP2C19 polymorphism on the clearance of R methadone. This relationship was not seen in our study. CYP3A4 and CYP3A5 activity scores were tested individually and as a cumulative score and were not found to be significant. Although CYP3A4 has the maximum role in the metabolism of the drug, it is generally less polymorphic and had a low frequency of heterozygous *22 (n = 6) and no homozygous *22 in our population. Interestingly, one of the intronic SNPs of CYP3A4 rs2246709 was found to be a better predictor of R methadone CLF and was therefore retained in the model. This relationship has been earlier shown in a study with a univariate analysis comparing concentrations at autopsy following a fatal overdose and the victim’s genotype. CYP2C18 has been shown earlier in in vitro studies to metabolize R and S methadone. Its SNP rs1042194 may have an effect on CL3, the clearance of both R and S‐EDDP.

The small sample size included in this study and strict inclusion/exclusion criteria may limit our ability to detect pharmacogenetic associations observed in other studies. For instance, rs2229109 (1199G>A) a possibly relevant variant of ABCB1 , when added as a covariate on clearance (CL3) of R and S‐EDDP, showed an improvement in OFV but had very few individuals with the variant that it was dropped from the final model (n = 2). A similar decision was taken with rs17180299, an important variant found to influence the clearance in a genomewide association study study. The NR1I subfamily receptors, the CAR and the PXR regulates numerous metabolic enzymes, including the CYP enzymes studied. The rs3003596, an SNP on CAR, shows a trend in relation with CLF of S methadone. Many of these relationships could indeed turn significant with a larger population, where the effects of multiple CYPs could be meaningfully separated. All the relationships that were significant from the univariate analysis and the covariate modeling are therefore given in Table 2 and Table S2.

In conclusion, this study reports an extensive analysis of R and S methadone and EDDP PK and relevant pharmacogenetic data in children undergoing major surgeries. CYP2B6 activity score and an intronic SNP on CYP3A4 were found to influence clearance. Plasma AAG concentrations and ORM1 genotype independently affected the volume of distribution of methadone. This research lays the foundation to further study these novel genotypes with a larger sample size and the addition of pharmacodynamic targets to best optimize the dose of methadone in children.

CONFLICT OF INTEREST

The authors declared no competing interests for this work.

AUTHOR CONTRIBUTIONS

B.W.A., S.K.Q., and B.R.O wrote the manuscript. S.K.Q., B.R.O., B.W.A., S.P., and S.S. designed the research. B.W.A., S.K.Q., B.R.O., M.A.H., A.R.M., H.G., and R.C.L performed the research. B.W.A., S.K.Q., B.R.O., and M.A.H. analyzed the data.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Fig S4

Click here for additional data file.

Code S1

Click here for additional data file.

Supplementary Material

Click here for additional data file.