Physiologically based pharmacokinetic modelling of atomoxetine with regard to CYP2D6 genotypes

Atomoxetine is a norepinephrine reuptake inhibitor indicated in the treatment of attention-deficit/hyperactivity disorder. It is primarily metabolized by CYP2D6 to its equipotent metabolite, 4-hydroxyatomoxetine, which promptly undergoes further glucuronidation to an inactive 4-HAT-O-glucuronide. Clinical trials have shown that decreased CYP2D6 activity leads to substantially elevated atomoxetine exposure and increase in adverse reactions. The aim of this study was to to develop a pharmacologically based pharmacokinetic (PBPK) model of atomoxetine in different CYP2D6 genotypes. A single 20 mg dose of atomoxetine was given to 19 healthy Korean individuals with CYP2D6*wt/*wt (*wt = *1 or *2) or CYP2D6*10/*10 genotype. Based on the results of this pharmacokinetic study, a PBPK model for CYP2D6*wt/*wt individuals was developed. This model was scaled to those with CYP2D6*10/*10 genotype, as well as CYP2D6 poor metabolisers. We validated this model by comparing the predicted pharmacokinetic parameters with diverse results from the literature. The presented PBPK model describes the pharmacokinetics after single and repeated oral atomoxetine doses with regard to CYP2D6 genotype and phenotype. This model could be utilized for identification of appropriate dosages of atomoxetine in patients with reduced CYP2D6 activity to minimize the adverse events, and to enable personalised medicine.

Atomoxetine is the first non-stimulant medication approved for the treatment of attention-deficit/hyperactivity disorder (ADHD) in children, adolescents, and adults. It has a potent and selective inhibitory effect on norepinephrine reuptake 1,2 . Guidelines from several organizations, including the American Academy of Pediatrics and the National Institute for Health and Care Excellence in the United Kingdom, recommend atomoxetine for the treatment of ADHD in patients of all ages.
One in vitro study found that the main metabolic pathway of atomoxetine is its oxidation to 4-hydroxyatomoxetine (4-HAT) by CYP2D6. To a minor degree, CYP2C19 is responsible for the biotransformation of atomoxetine to N-desmethylatomoxetine (N-DAT) 3 . 4-HAT is equipotent to atomoxetine, but promptly undergoes further glucuronidation to 4-HAT-O-glucuronide and therefore circulates in the plasma at very low concentrations 4 . The same study also indicated the effect of the CYP2D6 phenotype on the pharmacokinetics of atomoxetine. CYP2D6 is a highly polymorphic enzyme with more than 113 known variants (https://www. pharmvar.org/gene/CYP2D6). Among others, CYP2D6*3, *4, *5, and *6 are defective alleles that lack CYP2D6 enzyme activity. These alleles occur in ~20% of Caucasians. Approximately 7% of Caucasians are known to be poor metabolizers (PMs) due to presence of two non-functional alleles 5,6 . Oral administration of atomoxetine resulted in prolonged half-life (t 1/2 ) and 10-fold lower clearance (CL/F) in PMs compared to extensive metabolizers (EMs) with normal enzyme activity. These pharmacologic changes led to a 5.7-fold higher maximum plasma concentration (C max ) of atomoxetine and 7.8-fold greater area under the curve (AUC) in PMs compared to that in EMs 4 . Interestingly, the C max of the equipotent 4-HAT was only ~1.3% of that of atomoxetine 4 . Therefore, the pharmacodynamic effects of atomoxetine are primarily due to the parent drug, and the contribution of the metabolite is exiguous. Higher exposure to atomoxetine is not only correlated with its efficacy, but also with higher risk of adverse events, including tachycardia and elevated diastolic blood pressure 7 . The most common nonspecific adverse reactions of atomoxetine in CYP2D6 PMs include dry mouth, depression, tremor and insomnia 8 . Therefore, these patients are recommended to take lower doses of atomoxetine than EMs. In East Asia, the frequency of nonfunctional alleles is low. The most frequent variant in East Asians including Korean, Chinese, and Japanese populations is CYP2D6*10 (45%, 55%, and 38%, respectively) [9][10][11] . This allele is known to have decreased enzyme activity both in vitro and in vivo. Previous studies have demonstrated decreased biotransformation of atomoxetine in Chinese and Japanese patients with the CYP2D6*10/*10 genotype 12,13 . We also verified a significant influence of the CYP2D6*10 allele on the pharmacokinetic parameters of atomoxetine in a study with a larger sample size 14 . In that study, the C max and AUC from time 0 to infinity (AUC inf ) of the CYP2D6*10/*10 genotype group were 1.7-fold and 3.4-fold higher, respectively, than those of EMs with the CYP2D6*wt/*wt (*wt = *1 or *2) genotype (P < 0.001 for both). The t 1/2 and CL/F were also significantly different between the two groups (P < 0.001 for both). However, the increase in atomoxetine exposure in the CYP2D6*10/*10 group was not to the extent observed in the CYP2D6 PMs. Nevertheless, the risk of a higher incidence of concentration-related adverse events cannot be excluded in the CYP2D6*10/*10 group. Physiologically based pharmacokinetic (PBPK) modelling is a mechanistic approach to simulate the pharmacokinetics of xenobiotics based on their specific properties and mammalian physiology. A PBPK model is advantageous compared to conventional models because it considers individual anatomical and physiological parameters, including enzyme expression, and physicochemical compound data 15 . Therefore, PBPK modelling can be used to individualise therapies for special patient groups (such as children and the elderly) and to predict drug-drug interactions [16][17][18] . A PBPK model can also be useful to predict the effects of genetic polymorphism on the pharmacokinetics of xenobiotics 19 . In this study, we aimed to develop a PBPK model to predict the pharmacokinetics of atomoxetine in the different CYP2D6 genotype groups. This model was intended to demonstrate the possibility of simulating the pharmacokinetics of atomoxetine, taking into account the factors that affect it, including demographic characteristics, and CYP2D6 genotype 20 .

Results
Pharmacokinetic study. The major alleles found in the 399 Korean subjects were CYP2D6*1 and *2 with normal enzyme activity and *10 with decreased enzyme activity (34.3%, 13.4%, and 47.3%, respectively). The frequency of the nonfunctional allele CYP2D6*5 was 3.5%. Fourteen different genotypes were detected. The frequency of CYP2D6*wt/*wt, heterozygous CYP2D6*10, and homozygous CYP2D6*10 carriers were 23.3%, 45.1%, and 22.1%, respectively. There was only one PM with two defective alleles, while 15 subjects had the CYP2D6*5*/10 genotype. Among 19 healthy subjects enrolled in the pharmacokinetic study, 11 subjects had the CYP2D6*wt/*wt genotype (age range 19-25 years, BMI 18-26 kg/m 2 , weight 49-73 kg) and 8 subjects had the CYP2D6*10/*10 genotype (age range 19-25 years, BMI 18-23 kg/m 2 , weight 52-72 kg). The demographic characteristics were not significantly different between the groups, and all participants completed the study. There were no reported adverse events during and up to 10 days after the study completion. The plasma concentration-time curves and pharmacokinetic parameters of atomoxetine are summarized in Fig. 1 and Table 1. The C max , AUC 0-24 , AUC inf , t 1/2 , and CL/F were all significantly different between those with CYP2D6*wt/*wt and CYP2D6*10/*10 genotypes (P value < 0.0001 for all). Compared to individuals with CYP2D6*wt/*wt, those with CYP2D6*10/*10 had 1.5-fold higher C max , 3.1-fold higher AUC 0-24 and AUC inf , and 2.0-fold higher half-life. Oral clearance was 3.0-fold lower in the group with homozygous CYP2D6*10 than in those with the wildtype alleles. The t max was also significantly different between these groups (P value = 0.02).
Development of a PBPK model. Raw data from the pharmacokinetic study were used to develop the PBPK model. Values for physicochemical and pharmacokinetic (absorption, distribution, metabolism, excretion) parameters, both of which were used in the PBPK model, are listed in Table 2. Optimal value for absorption and distribution permeability was calculated by PK-Sim ® based on physicochemical data and plasma concentration profile of atomoxetine. The procedures on which the calculation is based on, are described by Thelen et al. and Kawai et al. [21][22][23] . Initially, all of the data were entered as provided in the literature. Specific renal clearance was calculated by the software from plasma clearance as input value. For the paramerisation of clearance, the Levenberg-Marquardt algorithms in the Parameter Identification tool of PK-Sim ® were used. The plama clearance values were first inputted individually in order to obtain similar total plasma clearance results to those observed in our pharmacokinetic study, and for the final model, the mean value was used. According to the Pharmaceuticals and Medical Devices Agency in Japan, the dissolution rate of atomoxetine capsules over 30 minutes, which was determined using the paddle method (described in the Japan Pharmacopoeia), was 85% or more. Therefore, 30 minutes was set as the dissolution time 24 . The first simulation showed a slightly lower plasma concentration-time profile than was previously observed. We calculated the organ-plasma partition coefficient based on the method developed by Berezhkovskiy, which is also available in PK-Sim ® 25 . This method calculates the organ-plasma coefficients using the volume fractions of the aqueous and organic subcompartments of the respective organ and plasma and peripheral drug eliminations are also considered. A sensitivity analysis for the parameters logP, fraction unbound, pKa, plasma clearance, amongst others, was also performed with the Sensitivity Analysis tool provided by the PK-Sim ® software. The sensitivity was determined as the mean of several sensitivities based on different relative variations, which were defined by multiplication of the value used for the simulation with variation factors. Both C max and AUC of atomoxetine were identified to be significantly sensitive to logP. After performing the parameter identification as described above, it was adjusted appropriately (logP = 3.8 instead of 3.9). The received curves reflected the observed data better after the adjustment, but they were shifted to the right. This was corrected by reducing the dissolution time to 15 minutes. Finally, the PBPK simulation data adequately fit the experimental data (Fig. 2). The visual comparison of the simulated and observed plasma concentration profiles demonstrated that the shape of the curves were similar. The simulated pharmacokinetic parameters AUC, C max and CL were also in accordance with the experimental results in Table 1. The absorption and distribution were in agreement with the literature. Atomoxetine is known to be well absorbed and the volume of distribution following oral administration is 2.33 L/kg 4 . In our simulation, the oral bioavailability was 96% and the volume of distribution at steady state (V ss ) was 2.18 L/kg. Apparent volume of distribution V d , which is calculated from the plasma curve according to V d = CL/λ, where λ is the terminal elimination rate, was 3.21 L/kg. In the next step, the achieved model was scaled to CYP2D6*10/*10 individuals. In addition to the individual demographic data and plasma clearance value, we modified the in vitro metabolic rate of CYP2D6 in the presence of liver microsomes. The value was reduced by around 90% of CYP2D6*1 activity as reported  previously 26 . The modified parameters used in the simulation of this genotype group are specified in Table 3. The plasma concentration-time profiles for CYP2D6*10/*10 are shown in Fig. 3. The observed and simulated pharmacokinetic parameters of CYP2D6*wt/*wt and *10/*10 are compared in Fig. 4.
Validation of the PBPK model. From previous published reports, human clinical PK data for atomoxetine in genotype-or phenotype-charaterised subjects were collected and used to validate the developed PBPK model of atomoxetine. The pharmacokinetics were simulated with the created populations (CYP2D6*10/*10, CYP2D6 EMs, PMs) and the adapted dosage regimens from each study. For the simulations of each genotype and phenotype group, only the in vitro metabolic rate was adjusted according to CYP2D6 genotype (Table 3). For the prediction of CYP2D6 EM, plasma clearance of 0.25 L/h/kg was used, which is slightly lower than that in our simulation     Table 4. All the prediction errors of the simulations were within two-fold error range.

Discussion
Atomoxetine is mainly metabolised by CYP2D6 to 4-HAT, and its pharmacokinetics are strongly affected by CYP2D6 genotype. The C max and AUC were substantially increased in CYP2D6 PMs compared to those in EMs.
The half-life was also prolonged by a factor of 3.7, and clearance was decreased by over 90% in PMs compared to EMs 4,20 . The efficacy and adverse events of a drug are both strongly dependent on drug exposure. In clinical trials, including open-label and long-term studies, PMs experienced more adverse events than EMs and subsequently had to discontinue the medication more frequently than those with adequate metabolism (11.2% of PMs vs. 6.3% of EMs) 8 . In contrast, more EMs (than PMs) stopped taking the drug because of lack of efficacy 7 . Therefore, in the United States and European countries, drug labels of atomoxetine recommend different dosage regimens based on CYP2D6 phenotype. In contrast, there is a low rate of CYP2D6 PMs in East Asian countries, and dosing adjustment recommendations are not present on atomoxetine drug labels in Korea and Japan, for instance.
Nevertheless, several studies have demonstrated the effect of the CYP2D6*10 allele, which is found mostly in East Asians, on the pharmacokinetics of atomoxetine [12][13][14] . The parameters of the pharmacokinetic study conducted in this study were considerably different in the two genotype groups, CYP2D6*wt/*wt and CYP2D6*10/*10.
The C max and AUC of individuals homozygous for CYP2D6*10 were significantly higher (1.53-fold in C max and 3.14-fold in AUC) than those of the CYP2D6*wt/*wt genotype group (P < 0.001 for both). The oral clearance was lower and the half-life was longer for those homozygous for CYP2D6*10 compared to those with the wildtype allele (P < 0.001 for both). While the C max and AUC increases in CYP2D6*10/*10 individuals were not to the extent of 5.7-fold and 7.8-fold respective increases observed in CYP2D6 PMs 4 , an increased risk of concentration-related adverse reactions in these individuals cannot be excluded. There were two other reports 12,13 comparing the pharmacokinetic parameters of atomoxetine in CYP2D6 EM and CYP2D6*10/*10, and their results also revelaled the significant increase in atomoxetine exposure in CYP2D6*10/*10 than in CYP2D6 EM although the differences in C max and AUC between two genotypes their studies were slightly less than those in this study 12,13 . PBPK models can be used to predict adequate drug dosing in order to avoid excessive difference in drug exposure: for instance, in individuals of different genotypes. Some previous studies have been conducted regarding the PBPK modelling of atomoxetine. Ball et al. 27 reported a PBPK modelling approach to predict the unbound CNS drug concentration in rats during preclinical drug development. As this PBPK modelling was conducted in rats, it has a limited application in humans. Dinh et al. 28 investigated the biotransformation of atomoxetine by CYP2D6 and other CYP isoforms as a preparation for developing a paediatric PBPK model. However, a PBPK model of of atomoxetine in humans was not reported in their study. Recently, after the first submission of our manuscript, Huang et al. reported a human PBPK model of atomoxetine 29 . Their PBPK model was developed using datasets obtained from literature, including human clinical pharmacokinetic data. However, their PBPK model failed to predict atomoxetine disposition in 100% of East Asian populations with CYP2D6 EM or CYP2D6*10/*10 genotypes or phenotypes, drug interaction studies, and specific population with renal disease or hepatic impairment. They validated the PBPK model using drug-specific acceptance criteria but did not succeed in predicting atomoxetine disposition in 100% of East Asian populations with CYP2D6 EM or CYP2D6*10/*10 genotypes or phenotypes, even if they validated by applying two-fold error criteria. Since we carried out our own pharmacokinetic studies, and thus all individual raw data was available, it was possible, in contrast to Huang et al., to perform individual simulations for model development. Our simulation results are not only consistent with the data from our pharmacokinetic study, but also with the previously published pharmacokinetic studies (Table 4). Among the nearly 400 genotyped individuals recruited for this study, there was only one CYP2D6 PM with the presence of two nonfunctional alleles (i.e. CYP2D6*5/*5), and this individual gave up participation in the study. Because of the low frequency, it was not possible to recruit PM subjects for the pharmacokinetic study. Therefore, we were unable to perform simulations on PM individuals. Instead, population simulations were performed, and the data were compared to that from the literature 4,30-32 . The developed PBPK model of atomoxetine is not coupled with models of its metabolites. Since the equipotent 4-HAT is immediately metabolised to a glucuronide form, its exposure is expected to only exert a minor effect on ADHD therapy. We were able to develop a suitable model in relation to CYP2D6 genotype after several adjustments. Regardless, there are a few limitations that must be considered. Only few papers have genotyped the subjects on CYP2D6, in which would describe the enzyme activity more accurately, whereas the other papers have phenotyped them. Thus, a group of CYP2D6 EMs might include lots of different combinations of CYP2D6 alleles, including those with decreased enzyme activity. Therefore, the potential variability of a phenotyped group is larger than a genotype group. Nevertheless, simulations with a plasma clearance of 0.25 L/h/kg, which is smaller than the plasma clearance used for CYP2D6*wt/*wt, described the pharmacokinetics of EM patients accurately. Inter-individual variability in pharmacokinetics is not only dependent on pharmacogenetics. As shown in our study, plasma clearance and/or protein binding vary greatly between individuals. Both of these factors can affect the pharmacokinetic characteristics of atomoxetine. Individual clearance was considered in this study. For protein binding, we used the mean value from the literature. Therefore, our model will not necessarily be able to predict the pharmacokinetics of each individual exactly, but by using the mean values for all parameters, the results of all individuals were within the two-fold error range. Thus the model might be helpful to identify the appropriate doses at the start of therapy by simply inputting the demographics of the patient, which is a factor that influences the pharmacokinetics of atomoxetine, with the presumption that the genotype of the patient is known. The presented model can also be further adjusted to predict drug-drug interactions. Belle et al. found that the C max and AUC of atomoxetine in EM subjects increased by 3.5-fold and 6.5-fold, respectively, after co-administration with paroxetine, a strong CYP2D6 inhibitor 30 . These increases are similar to those observed in PMs. Therefore, the dose adjustments in PM patients may also be similarly applied to those with concomitant use of potent CYP2D6 inhibitors. Atomoxetine is a drug that is primarily administered to children and adolescents. Therefore, it is important that data are generalisable to this population. The prediction error of the simulation with paediatric patients was larger than for the other simulations, but still within the two-fold error range. But there was only one dataset in the literature which provided pharmacokinetics of paediatric patients in relation to CYP2D6 33 . Further simulations are needed for validation. In conclusion, we have developed a PBPK model of atomoxetine, which predicts the plasma concentration-time profiles of CYP2D6 EMs with CYP2D6*wt/*wt genotype and also of those with reduced enzyme activity due to CYP2D6*10 alleles and those with defective alleles. This model appropriately represents the pharmacokinetics of atomoxetine, considering the factors that affect the pharmacokinetics of atomoxetine, including demographic characteristics, and CYP2D6 genotype 20 and therefore may serve as a basis for further model development, i.e. for the simulation of the drug's pharmacokinetics in special patient groups, such as those with hepatic impairment, or for the prediction of the effects of co-administered drugs.

Subjects.
A total of 399 healthy Koreans subjects were recruited, and their genomic DNA was extracted from 10 mL whole blood samples using the Wizard ® Genomic DNA Purification Kit (Promega, Madison, WI, USA).
Isolated DNA was used to determine the CYP2D6 and CYP2C19 genotypes. Detection of CYP2D6*2, *5, *10, *XxN and CYP2C19*2, *3, *17 alleles were performed via polymerase chain reaction restriction fragment length polymorphism (PCR-RFLP) or long range PCR using previously described methods 34 . Nineteen healthy Korean subjects with CYP2D6*wt/*wt or CYP2D6*10/*10 genotype were randomly selected for the pharmacokinetic study of atomoxetine. None of them had a CYP2C19 PM (CYP2C19*2/*2, *2/*3, *3/*3) genotype. Subjects were healthy according to medical history, physical examination, and routine laboratory tests (blood chemistry, hematology, and urine analysis). They were not permitted to use any medications, alcohol, or caffeine for 10 days prior to or during the pharmacokinetic study. All subjects provided verbal and written informed consent.
Study protocol. The  Determination of plasma concentration. Plasma concentrations of atomoxetine and its metabolites (4-HAT and N-DAT) were determined and validated using liquid chromatography/tandem mass spectrometry following a previously described protocol 35,36 . Briefly, after extraction with methyl t-butyl ether, analytes in human plasma were chromatographically separated using a Luna C 18 column ( . The C max values and time to reach C max (t max ) were derived from the experimentally measured values. Area under the curve was calculated using the log-linear trapezoidal approximation. AUC inf was extrapolated by adding C t /k e to AUC 0-t , where C t is the last measured plasma concentration, and k e is the elimination rate constant. The elimination rate constant was estimated from the terminal phase using log-linear regression. The half-life was assessed using the formula t 1/2 = ln 2 /k e . The apparent oral clearance was calculated with the following formula: CL/F = dose/AUC inf . All data are presented as mean ± standard deviation (SD), except t max , which is expressed as median and range. The differences in pharmacokinetic parameters between the two genotype groups were assessed using the Mann-Whitney rank sum test. SigmaPlot ® version 12 (Systat Software, Chicago, IL, USA) was used for the statistical analysis. P values < 0.05 were considered statistically significant (α = 0.05). The power of the study was at least 80%. The sample size was estimated to be sufficient to detect 100% greater AUC of atomoxetine among the two genotype groups. The sample sizes were calculated with the Power and Sample Size Program, PS (version 3.0.17) 37 .
PBPK model development workflow. The modelling software PK-Sim ® version 7.1.0 (Open Systems Pharmacology Suite) was used to develop this PBPK model. Demographic information on age, weight, and height were collected during the recruitment for the pharmacokinetic study. Compound specific physicochemical data were obtained from the literature and integrated into the simulation. Values describing the absorption, distribution, metabolism and excretion (ADME) were obtained from the literature or calculated in silico with PK-Sim ® .
All other parameters and settings were maintained at default values. If necessary, the input parameters were modified in the appropriate dimensions after performing a sensitivity analysis and parameter identification to better fit the simulated and experimented data. Once the PBPK model for the CYP2D6*wt/*wt genotype group was generated, the CYP2D6*10/*10 genotype model was scaled by adjusting the individual biometric data, the metabolic enzyme activity, and the clearance value. Both simulations were validated with the experimental results from diverse studies. The populations (n = 10) were created using PK-Sim ® . Different ethnic groups were addressed by selecting the ethnicity provided by the software. The demographic characteristics (age, BMI, weight) provided in the literature were set as upper and lower limits of the population. The proportion of female volunteers was defined equally to each study. The software then randomly selected parameters to be included. We included all observed data from the available literature, which included the pharmacokinetics of adult and paediatric CYP2D6 EM, CYP2D6*10/*10 and CYP2D6 PM. The administration protocols were equally adapted. We created populations of poor metabolisers with two defective CYP2D6 alleles as described above and similarly compared these data. Paediatric CYP2D6 EM and PM were simulated individually, as all demographic characteristics, their CYP2D6 genotype and administered dosage were known 33 . This model was evaluated by visual comparison and in terms of the prediction error (PE), calculated as PE = (simulated-observed)/observed × 100. A relative error smaller than a two-fold error was accepted, which is widely applied by the pharmaceutical industry in drug discovery and development 38 . If necessary, model parameter was refined during model verification, as recommended by the United States Food and Drug Administration (FDA) guidance.
Data availability. The datasets generated and analysed in the current study are available from the corresponding author on reasonable request.