- 1School of Pharmacy, China Pharmacokinetic University, Nanjing, China
- 2Office of China National Narcotics Control Commission, China Pharmacokinetic University Joint Laboratory on Key Technologies of Narcotics Control, Beijing, China
- 3Beijing Narcotics Control Technology Center, Beijing, China
- 4Key Laboratory of Drug Monitoring and Control, Drug Intelligence and Forensic Center, Ministry of Public Security, Beijing, China
Introduction: Fentanyl analogs, as emerging new psychoactive substances (NPS), pose a global public health threat due to widespread abuse, high toxicity, and frequent overdose fatalities. However, their structural diversity and scarce experimental pharmacokinetic (PK) data hinder hazard and abuse risk assessment. Conventional physiologically based pharmacokinetic (PBPK) models for these analogs are limited by reliance on time-consuming in vitro experiments or error-prone interspecies extrapolation for key parameters (e.g., tissue/blood partition coefficient, Kp).
Methods: To address this, we developed and validated a QSAR-integrated PBPK framework (QSAR: Quantitative Structure-Activity Relationship) for predicting human PK of fentanyl analogs. The workflow included: (1) Validating the framework via intravenous β-hydroxythiofentanyl in Sprague-Dawley rats (QSAR-predicted Kp via Lukacova method, GastroPlus® modeling); (2) Comparing Kp accuracy (literature in vitro data, QSAR, interspecies extrapolation) in rat/human fentanyl PBPK models; (3) Predicting PK and tissue distribution (plasma +10 organs including brain/heart) of 34 human fentanyl analogs.
Results: Key results: (1) For β-hydroxythiofentanyl, all predicted rat PK parameters (area under the plasma concentration-time curve from time zero to the last measurable time point [AUC0-t], teady-state volume of distribution [Vss], and elimination half-life [T1/2]) of rats fell within a 2-fold range of the experimental values; (2) In human fentanyl models, QSAR-predicted Kp improved accuracy (Vss error: >3-fold [extrapolation] vs. <1.5-fold [QSAR]) (3) Among 34 analogs, eight (e.g., p-fluorofentanyl); had brain/plasma ratio >1.2 (vs. fentanyl’s 1.0), indicating higher CNS penetration and abuse risk.
Discussion: This study demonstrates that the QSAR-PBPK framework enables rapid prediction of human pharmacokinetics (PK) for understudied fentanyl analogs without relying on scarce experimental data. For structurally similar, clinically characterized analogs (e.g., sufentanil, alfentanil), predictions of key PK parameters (e.g., T1/2, Vss) fall within 1.3–1.7-fold of clinical data, supporting the framework’s utility for generating testable hypotheses about the PK of understudied analogs. It not only fills the data gap for fentanyl analog hazard assessment but also provides a scalable modeling strategy for PK evaluation of other NPS or illicit drugs.
1 Introduction
Fentanyl, a potent synthetic opioid receptor agonist, was initially developed to exert robust analgesic effects. Clinically, pharmaceutical fentanyl is widely employed in the management of various pain conditions, post-operative analgesia, and as an adjunct to anesthetics during surgical procedures (Patocka et al., 2024). Alarmingly, the 2024 World Drug Report highlights that fatalities attributed to fentanyl overdoses in North America surged to unprecedented levels during the COVID-19 pandemic (UNODC, 2024). Compounding this challenge, the relative ease of synthesis and structural modification of fentanyl has led to the proliferation of numerous analogs with diverse biological activities. Unfortunately, the pace of research has failed to keep abreast of the emergence and rapid spread of these new fentanyl analogs. For example, after Sweden classified acetylfentanyl and butyrfentanyl as controlled substances and implemented bans, these derivatives disappeared from online marketplaces, only to be replaced by variants such as furanylfentanyl and 4-methoxybutyrfentanyl (Helander et al., 2017). Consequently, a rapid understanding of the pharmacokinetic profiles of emerging analogs is of paramount importance for formulating relevant therapeutic strategies and establishing effective regulatory policies.
Physiologically based pharmacokinetic (PBPK) models represent advanced computational tools that predict the concentration-time profiles of compounds in plasma across different species. This is accomplished by integrating the physicochemical properties of the compound with the physiological characteristics of the target species (Sager et al., 2015). Key parameters such as logD, pKa, and the unbound fraction to plasma proteins (Fup) can be obtained either through in vitro measurements or predicted using quantitative structure-activity relationship (QSAR) models (Miller et al., 2019). Moreover, by leveraging these input parameters, PBPK models can predict chemical concentration changes in various tissues, enabling the rapid screening of substances with a quick onset of action and facilitating the elucidation of the pharmacological properties of specific compounds.
However, traditional methods for evaluating model performance rely on pharmacokinetic (PK) data specific to the target chemical, rendering them inadequate for fentanyl derivatives that have not been sufficiently studied (EPA, 2006). Previous research has explored the hypothesis that PBPK models fully developed for a target chemical (with unavailable PK data) can be evaluated using PK data from its structural or functional analogs (with available PK data) (Ellison, 2018). Fentanyl, being extensively documented in clinical literature, provides a rich dataset that can be used to assess the accuracy of models for fentanyl analogs.
Notably, our approach constitutes a notable advancement in the field of PBPK modeling for fentanyl analogs compared with prior investigations. Earlier studies were often constrained by two key limitations: they either relied exclusively on in vitro assays for parameterization—a process inherently characterized by high time and resource consumption—or restricted the application of QSAR predictions to the estimation of opioid-opioid receptor affinity (Floresta et al., 2019), rather than leveraging QSAR for the prediction of PBPK-essential critical parameters (e.g., logD, pKa, plasma protein unbound fraction). In distinct contrast, our work integrates the predictive capacity of QSAR with PBPK modeling frameworks. Through a systematic comparative analysis of parameters derived from in vitro measurements versus QSAR predictions, we successfully identified the optimal parameter source, which directly contributes to enhanced model accuracy. Furthermore, the QSAR integration in our study represents a significant improvement over existing models. By harnessing QSAR to predict PBPK-essential parameters, we reduce reliance on scarce in vitro data and accelerate the modeling process. This not only addresses the data gap for understudied fentanyl analogs but also enhances the speed of pharmacokinetic prediction without compromising accuracy—a claim validated by comparisons with measured experimental data.
The objectives of this study are threefold: 1) to evaluate the accuracy of the PBPK model based on QSAR-predicted data using measured PK data of beta-hydroxythiofentanyl in rats; 2) to develop PBPK models for fentanyl in both rats and humans, with modeling parameters derived from in vitro measurements and QSAR predictions respectively, and to compare the accuracy of these parameters during the modeling process; 3) to predict the pharmacokinetics of 34 fentanyl analogs in human plasma and various tissues, thereby bridging the existing data gap in this field. The study design and workflow are illustrated in Figure 1.

Figure 1. Experimental flowchart for predicting human ADME/PK of fentanyl analogs using PBPK modeling.
2 Materials and methods
2.1 Materials
2.1.1 In silico programs
Fentanyl and its derivatives were identified through a comprehensive search of peer-reviewed journal publications (PubMed, National Library of Medicine). The structural formulas of fentanyl analogs were obtained from the PubChem Website based on the compound name (https://pubchem.ncbi.nlm.nih.gov). QSAR models predicted the physiochemical and pharmacokinetic properties of fentanyl analogs in ADMET Predictor v. 10.4.0.0 (AP, Simulations Plus, Inc.). PBPK modeling and simulations were established using GastroPlus v. 9.8.3 (GP, Simulations Plus, Inc. Lancaster, CA). The PK parameters were estimated using Phoenix WinNonlin software (version 8.3).
2.1.2 Animals
Male Sprague-Dawley (SD) rats, aged 6–8°weeks old, were procured from SPF Biotechnology Co.,Ltd (Beijing, China). The animals were housed in a controlled environment with humidity maintained at 50% ± 10% and temperature at 25 C ± 2 C, under a 12 h light/dark cycle They had ad libitum access to food and water. All experiments were conducted according to protocols approved by the Institutional Welfare and Ethics Committee for Laboratory Animals, Key Laboratory of Drug Monitoring and Control.
2.1.3 Chemicals and reagents
Beta-hydroxythiofentanyl (content ≥98%) was provided by the Drug Intelligence and Forensic Center of the Ministry of Public Security, China. 0.9% saline solution was obtained from Shandong Qidu Pharmaceutical Co., Ltd. Acetonitrile and formic acid were sourced from Sigma-Aldrich (Saint Louis, MO). All other reagents utilized in the experiment were commercially available and of analytical or HPLC grade. The triple quadrupole mass spectrometric detector (6,500+; AB SCIEX Technologies) was used for LC-MS/MS analysis.
2.2 Methods
2.2.1 Development and validation of the β-hydroxythiofentanyl PBPK model in rats
To assess the modeling approach based on QSAR predictions, our research developed the PBPK model for intravenous β-hydroxythiofentanyl in rats and compared the experimentally obtained concentration and time profiles and pharmacokinetic parameters with the predicted data.
Time points for the study were 0 min, 15 min, 45 min, 60 min, 90 min, 120 min, 180 min, and 240 min post-dosing with 7 μg/kg iv. dose of β-hydroxythiofentanyl. For the pharmacokinetic assessment, 400 μL of blood was collected and subsequently centrifuged at 4,000 rpm for 5 min. The plasma was then carefully transferred into microcentrifuge tubes and stored in a freezer at −20 °C until further analysis. On the day of testing, the plasma samples were analyzed using LC-MS. Non-compartmental analyses were conducted utilizing Phoenix WinNonlin software to estimate pharmacokinetic parameters.
For the PBPK model of β-hydroxythiofentanyl constructed using QSAR predictions, the molecular structure of β-hydroxythiofentanyl was input as the core descriptor. The tissue/blood partition coefficient (Kp) shown in Table 1, a critical parameter governing tissue distribution in PBPK modeling, was predicted using the Lukacova method. It was a structure-driven QSAR approach widely applied for estimating tissue partition coefficients of small-molecule compounds based on their structural features. The QSAR-predicted Kp values were incorporated into GastroPlus® software to construct the PBPK model for β-hydroxythiofentanyl. For β-hydroxythiofentanyl, the systemic clearance (CLsys), a critical parameter for PBPK model parameterization, was derived directly from in vivo experimental data, rather than QSAR prediction. The CLsys value obtained in the β-hydroxythiofentanyl model prediction section is derived from the CLsys value predicted by GastroPlus® software after inputting the experimentally measured CLsys value and the QSAR-predicted Kp value.

Table 1. The tissue/blood partition coefficient (Kp) for β-hydroxythiofentanyl and fentanyl in rats and humans.
Finally, to validate the accuracy of this QSAR-based PBPK model, the model-predicted pharmacokinetic (PK) profiles of β-hydroxythiofentanyl were systematically compared with experimentally measured PK data of the compound. This comparative analysis allowed quantification of the agreement between predicted and observed PK outcomes, thereby verifying whether the QSAR-based parameterization (Lukacova method-derived Kp) could support reliable PBPK modeling of β-hydroxythiofentanyl.
2.2.2 Development and validation of the fentanyl PBPK model in rats and humans
PBPK models for fentanyl were established to characterize the pharmacokinetic (PK) profiles of fentanyl following intravenous (IV) administration in rats and humans, respectively. In PBPK modeling, the Kp is recognized as a critical parameter for accurate prediction of tissue concentrations, as it directly governs the distribution of fentanyl between systemic circulation and peripheral tissues. Conventionally, Kp values are determined via multiple approaches, including in vivo animal PK studies, in vitro equilibrium dialysis (for measuring plasma protein binding and tissue affinity), and QSAR modeling (Yau et al., 2020). For the rat PBPK model, Kp values were initially derived from previously published literature data to ensure consistency with established experimental findings (Björkman et al., 1990). To further evaluate the reliability of QSAR-predicted Kp for model parameterization, an additional rat PBPK model was constructed using Kp values predicted via the Lukacova method—a structure-based QSAR approach that estimates tissue partition coefficients based on the molecular structure of the compound. This dual-parameter-source design allowed direct comparison of model performance between literature-derived and QSAR-predicted Kp.
A major challenge in studying fentanyl and its analogs is the scarcity of clinical PK data, as direct human studies are often constrained by ethical and safety considerations. To address this gap, PBPK models enable interspecies extrapolation of PK data from animal models (e.g., rats) to predict human absorption, distribution, metabolism, and excretion (ADME) profiles. In the present study, two distinct human fentanyl PBPK models were developed, differing only in their Kp parameter sources: QSAR-predicted Kp: Human Kp values were generated via the same Lukacova QSAR method applied to the rat model; Human Kp values were predicted from rat Kp values using the following formula:
Physicochemical properties of fentanyl (logP, pKa, and solubility profiles) were predicted by GastroPlus software [fup (rat) = 8.3, fup (human) = 25.5, BP(rat) = 1.01, BP(human) = 1.01]. The human PBPK model was administered at a dose of 0.1 mg at 70 kg.
For the CLsys value of fentanyl, the observed value was obtained from experimentally measured data reported in the literature. Subsequently, the experimentally measured CLsys values corresponding to rats and humans, along with the Kp values derived from QSAR (Quantitative Structure-Activity Relationship) predictions and interspecies extrapolation, were separately input into the model for prediction. The corresponding predicted CLsys values were then generated as outputs. The performance of the models was evaluated by comparing the predicted plasma concentration-time (C-t) profiles with experimental clinical data (for humans) and literature-derived PK data (for rats). The primary objective of this validation was to verify whether QSAR-based Kp parameterization yields a human PBPK model with higher predictive accuracy than the model relying on interspecies extrapolation of Kp (Lim et al., 2012, YK. et al., 2024).
2.2.3 Development and application of QSAR-Based PBPK models for 34 fentanyl analogs in humans
Leveraging the established QSAR-based PBPK modeling framework, PBPK models were constructed for 34 fentanyl analogs to characterize their disposition following intravenous administration in humans.
The model development workflow was standardized across all analogs: First, QSAR modeling was employed to generate two core sets of input parameters: (1) tissue/blood partition coefficients (Kp) for each analog in human, predicted via the Lukacova method (consistent with the fentanyl model); and (2) key physicochemical properties of the analogs (e.g., logP, aqueous solubility), predicted using ADMET Predictor™. These QSAR-derived parameters were then integrated into the human PBPK modeling platform (GastroPlus® software) to construct species-specific PBPK models for each fentanyl analog. It is noteworthy that the predictions for 34 fentanyl analogs in this study all utilized the experimentally measured CLsys value of fentanyl in humans, as reported in the literature. This approach was adopted to achieve the effect of rapidly predicting the relevant data of this class of substances based on one analog. After inputting fentanyl’s CLsys value and other parameters into the model, the predicted CLsys values adjusted by the model can be generated as outputs.
This modeling subsequently enabled the prediction of their concentration and time profiles in plasma and ten tissues and organs, including the brain, heart, and adipose. Finally, the predictive capability of the model was evaluated using clinical data pertaining to fentanyl.
3 Results
3.1 QSAR-based PBPK results show good agreement with experimental measurements by β-hydroxythiofentanyl
According to the findings presented in 3.1, a PBPK model was developed for rats that were administered β-hydroxythiofentanyl via intravenous injection, using Kp predicted through QSAR. The data obtained from the experiment versus predicted plasma concentration and time profiles after intravenous administration are plotted in Figure 2. Based on visual inspection, the intravenous administration of β-hydroxythiofentanyl demonstrated great agreement between the data from experiment and the predicted. The model parameters and pharmacokinetic parameters from experimented and predicted are listed Table 2. The experimented and predicted values for Vss and AUC0-t are compared and illustrated in Figure 2. The predicted values for all pharmacokinetic parameters fell within a 2-fold error margin of the experimented values.

Figure 2. (A) Predicted and experimented concentration-time profiles of β-hydroxythiofentanyl following intravenous administration. (B,C) The experimented versus predicted graphs for pharmacokinetic parameters of β-hydroxythiofentanyl, including Vss and AUC0-t. The orange dot represents data of Kp from QSAR.

Table 2. The model parameters and pharmacokinetic parameters from experimented and predicted of β-hydroxythiofentanyl.
3.2 QSAR-based PBPK results in human exhibit higher accuracy than those from interspecies extrapolation in human PBPK model for fentanyl
The PBPK model of fentanyl was respectively developed in rats and humans, with the injectable dose derived from the commonly utilized analgesic dosage of fentanyl in humans, which is 0.1 mg. This human dosage was converted to a corresponding amount based on body surface area, resulting in an administration of 1.75 μg to the rat subjects (Nair and Jacob, 2016). Fentanyl possesses a low molecular weight and high lipid solubility, which facilitate its diffusion across cellular membranes. Consequently, the model was developed using the perfusion rate-limited tissue model approach.
The observed versus predicted plasma concentration and time profiles after intravenous administration are plotted in Figure 3. Based on visual inspection, improved concordance between observed data from the report and predicted data constructed using Kp reported in the research in rat model, as well as between the observed data from the report and predicted data developed utilizing QSAR predicted Kp in human model. The model parameters, clinical study information and pharmacokinetic parameters are listed Table 3. The observed and predicted values for Vss and AUC0-t are compared and illustrated in Figure 4. The predicted values for most pharmacokinetic parameters fell within a 2-fold error margin of the observed values. However, the Vss derived from human models through interspecies extrapolation exceeded this range, although it remained within a 3-fold error margin.

Figure 3. Predicted and observed concentration-time profiles of fentanyl following intravenous administration (A) rat model of Kp from report; (B) rat model of Kp from QSAR prediction; (C) human model of Kp from interspecies extrapolation; (D) human model of Kp from QSAR prediction.

Table 3. The model parameters, clinical study information and pharmacokinetic parameters of fentanyl in rats and humans.

Figure 4. The observed versus predicted graphs for pharmacokinetic parameters of fentanyl (A) Vss of rat model; (B) area under the plasma concentration-time curve from time zero to time t, AUC0-t of rat mode; (C) Vss of human model; (D) AUC0-t of human model. The solid line denotes the unity line, where the ration of predicted to observed values equals 1. The dotted line indicates a two-fold error margin. The green triangle represents data of Kp from report; The red square represents of data Kp from QSAR prediction; The orange dot represents data of Kp from interspecies extrapolation.
3.3 Human PBPK modeling of 34 fentanyl analogs derived from QSAR-predicted results
The PBPK model for 34 fentanyl analogs following intravenous administration in the human was developed using the Kp predicted by QSAR. The predicted model parameters and pharmacokinetic parameters predicted are listed in Table 4. When the AUC of the 34 fentanyl analogs were compared to plasma data for fentanyl in humans, all AUCs fell within a 1.3-fold range of one another. Figure 5 illustrates the Cmax of fentanyl analogs in the brain, heart, adipose and liver, and brain/plasma ratio. The results showed that drug concentrations in the brain and heart were significantly higher than in plasma.

Figure 5. Predicted distribution in tissues of fentanyl analogs following intravenous administration (A) Cmax of adipose, brain, heart and liver; (B) the ratio of brain and plasma, green bars indicate ratios lower than fentanyl, blue bars indicate ratios higher than fentanyl.
4 Discussion
4.1 Assessment of the viability of the modeling approach
An important factor that hinders the application of PBPK models in drug research is the substantial amount of data required for model construction. This encompasses the experimental determination of tissue affinities, particularly Kp, which can be both costly and time-consuming. The most recent research on the modeling of PBPK for fentanyl and its analogs employs Kp derived from mathematical formulas or assumes that the same data applies to both humans and rats (Han et al., 2025; Rodgers et al., 2005). In addition, the extrapolation of animal data to humans is a widely accepted method in PBPK modeling. One study has successfully developed a PBPK model for cocaine and its metabolite benzoylcholine in humans based on extrapolated data from rats (Bravo-Gómez et al., 2019). Therefore, our study employed experimental measurements, QSAR predictions, and interspecies extrapolation to derive Kp for modeling purposes. To verify the applicability of the QSAR prediction modeling method for Kp to fentanyl analogs, rather than being restricted solely to fentanyl, our study employed this approach to develop a rat PBPK model of β-hydroxythiofentanyl, then evaluated the model against experimentally determined drug-time profiles and pharmacokinetic parameters. By comparing several fentanyl models that we have developed, it is observed that the experimentally determined Kp is preferred in the rat model, whereas the QSAR-predicted Kp is favored in the human model. The results demonstrated that the model exhibited high predictive accuracy and confirmed the feasibility of this approach. Eventually, our research constructed 34 fentanyl human PBPK models for fentanyl analogs following this method.
4.2 Analysis of model results
4.2.1 Comparison of data on fentanyl analogs with existing studies
Intravenous fentanyl and Sufentanil are primarily utilized for general anesthesia, while alfentanil and remifentanil serve the purposes of analgesia and sedation. Notably, remifentanil is particularly well-suited for short-term or outpatient procedures (Ziesenitz et al., 2018). Given the specific medical applications of these three fentanyls, their pharmacokinetics have been extensively investigated in existing studies. In contrast, there is a significant lack of experimental data regarding other illicitly manufactured fentanyl analogs in both animal models and human volunteers. Consequently, we selected alfentanil, Sufentanil, and remifentanil as representative drugs to evaluate the predictive accuracy of our developed model using PK parameters documented in the current research (details shown in Table 5).

Table 5. The reported observed and predicted values of selected PK (logP and Fup) and QSAR (CL and T1/2) parameters for Sufentanil, Alfentanil, and Remifentanil.
In a study on sufentanil, the pharmacokinetic characteristics of 10 surgical patients who received an intravenous dose of 5 μg/kg sufentanil were reported. The results showed a mean half-life (T1/2) of 164 ± 22 min (converted to hours: 2.73 ± 0.37 h) and a mean steady-state volume of distribution (Vdss) of 1.7 ± 0.2 L/kg (Bovill et al., 1984). Based on a standard body weight of 70 kg, the absolute Vdss value in this study was calculated as 119 ± 14 L. The T1/2 of sufentanil predicted in this study was 3.652 h, with an error ratio of 1.34 (3.652/2.73) relative to the literature value, within the acceptable error range of 1.3–1.7 fold. The predicted steady-state volume of distribution (Vss) was 308.601 L, resulting in an error ratio of approximately 2.59 when compared to the literature-derived absolute value (119 ± 14 L). This discrepancy may be attributed to two main factors: first, the standardized setting of human physiological parameters in the model (e.g., uniform 70 kg body weight, standard tissue blood perfusion rates) versus inherent inter-individual variability in patient physiology (e.g., age, body weight distribution, organ function status); second, differences between the experimental measurement method for Vdss in the literature and the model’s prediction algorithm (derived from QSAR-predicted tissue partition coefficients, Kp). Despite the numerical difference, both values reflect the large volume of distribution of sufentanil, indicating its extensive tissue distribution in vivo, which confirms consistency in qualitative trends.
For alfentanil, a study on healthy volunteers reported that following intravenous administration of 170 μg alfentanil (equivalent to approximately 2.4 μg/kg for a 70 kg individual), the mean T1/2 was 1.21 h (Bower and Hull, 1982). The T1/2 of alfentanil predicted in this study was 0.93 h, with an error ratio of 1.75 (1.63/0.93) relative to the literature value—approaching the upper limit of the acceptable 1.7-fold error range. The predicted systemic clearance (CLsys) was 53.931 L/h. Although no direct measurement of alfentanil clearance was reported in the literature, clearance was estimated using the relationship between AUC0-t (area under the plasma concentration-time curve) and dose (CL = dose/AUC), yielding an estimated clearance of approximately 31.44 L/h (Bower and Hull, 1982). The error ratio between the predicted and estimated clearance values was 1.71 (53.931/31.44), which also fell within the acceptable range. Both values consistently indicate the relatively rapid clearance of alfentanil.
In a study on remifentanil, surgical patients received intravenous doses of 2, 5, 15, and 30 μg/kg remifentanil, and blood samples were analyzed. The results showed that T1/2 increased with dose, measuring 10.19, 14.35, 15.67, and 20.47 min (converted to hours: 0.17, 0.24, 0.26, and 0.34 h), respectively (Westmoreland et al., 1993). The T1/2 of remifentanil predicted in this study was 0.571 h, with error ratios ranging from 1.68 to 3.36 (0.571/0.34 to 0.571/0.17) relative to the literature values for each dose group. It is important to note that the subjects in this literature study were surgical patients, and remifentanil—an ultra-short-acting opioid—undergoes metabolism that is significantly influenced by esterase activity, hepatic/renal function, and surgical stress. In contrast, the model in this study was constructed based on physiological parameters of healthy individuals and did not account for the effects of pathological conditions or stress responses on metabolic enzyme activity, which likely contributed to the larger discrepancy between predicted and literature values. Additionally, regarding the AUC and Cmax (peak plasma concentration) of remifentanil, only one 2004 study reporting data from ICU patients with renal dysfunction was identified (Pitsiu et al., 2004). These patients had impaired renal function; although remifentanil is primarily metabolized by non-specific esterases (with minimal impact from renal function), the overall physiological state of the patients (e.g., circulatory stability, metabolic enzyme expression levels) differed significantly from that of healthy individuals. This made the data poorly comparable to the predictions of the healthy human-based model in this study, so these parameters were not included in quantitative comparisons.
Beyond the three aforementioned drugs, recent literature has reported pharmacokinetic data for other fentanyl analogs, but most focus on non-intravenous administration routes (e.g., sublingual, subcutaneous implantation). For example, one study (van de Donk et al., 2018) reported the bioavailability and AUC of sublingual fentanyl wafers; however, this route involves drug absorption through the oral mucosa, which differs significantly from the intravenous route (no first-pass effect, simple absorption process) used in this study in terms of pharmacokinetic characteristics. Another study measured the PK parameters of analogs such as acetylfentanyl and butyrylfentanyl via subcutaneous injection (Canfield and Sprague, 2024), which involves slow drug diffusion and absorption in subcutaneous tissue—leading to differences in the timing of peak plasma concentration and AUC calculation methods compared to intravenous injection. These differences prevent direct quantitative comparison with the predictions of this study. Therefore, such literature data from non-intravenous routes were not included in the comparative analysis.
To further verify the reliability of the model’s input parameters, this study supplemented a comparison between QSAR-predicted key physicochemical properties (e.g., logP, Fup) and experimentally measured values. As shown in Table 5: 1) For sufentanil, the experimental logP value was 3.95, and the QSAR-predicted value was 3.85, with an error of only 2.5%. 2) For alfentanil, the experimental logP value was 2.16, and the QSAR-predicted value was 2.26, with an error of 4.6%. 3) For remifentanil, the experimental logP value was 1.4, and the QSAR-predicted value was 1.95, with an error of 39.3% (primarily due to the presence of ester and morpholine groups in the remifentanil molecule, which cause variability in experimentally measured logP values across different solvent systems).
In terms of Fup: 1) The experimental Fup value for sufentanil was 7.5%, and the QSAR-predicted value was 8.91%, with an error of 18.8%. 2) The experimental Fup value for alfentanil was 12.75%, and the QSAR-predicted value was 15.79%, with an error of 23.8%.
Overall, the QSAR-predicted values for logP (except for remifentanil) and Fup showed minimal deviation from experimental values, particularly for logP, which exhibited high predictive accuracy. As a key parameter influencing drug tissue partition coefficients (Kp) and transmembrane transport capacity, the accurate prediction of logP directly ensures the reliability of the PBPK model’s simulations of drug distribution and clearance processes. Even for remifentanil—where logP prediction showed a larger deviation—the QSAR-predicted Fup value (26.57%) was consistent with the qualitative characteristic of “low plasma protein binding of remifentanil” reported in the literature, still providing a reasonable basis for input parameters in the model.
In summary, the errors between the predicted PK parameters (T1/2, CLsys) of sufentanil, alfentanil, and remifentanil in this study and the literature experimental values mostly fell within the 1.3–1.7 fold range. Additionally, the QSAR-predicted values of key physicochemical properties showed minimal deviation from experimental values, indicating high predictive accuracy of the model. For parameters with moderate deviations (e.g., Vss of sufentanil), the causes can be reasonably explained by factors such as differences in administration routes, standardized physiological parameter settings, and experimental conditions—with consistent qualitative trends with the literature. This result significantly enhances the credibility of the PK predictions for other fentanyl analogs lacking experimental data in this study, providing a reliable model basis for subsequent assessments of the abuse potential of these analogs (e.g., evaluating central nervous system penetration based on brain/plasma concentration ratios).
4.2.2 Impact of the structural characteristics of fentanyl analogs on pharmacokinetics
New psychoactive substances (NPS) are frequently synthesized in clandestine laboratories with the intention of chemically modifying controlled drugs to circumvent legal regulations. Illegally manufactured substances resembling fentanyl are among the contributing factors to the global rise in fatalities associated with the NPS. Fentanyl analogs are frequently marketed as fentanyl itself, offered as substitutes for other substances, or incorporated into counterfeit prescription medications (Armenian et al., 2018). Due to the insufficient PK and PD evaluation of these analogs, potential users and others who may be exposed to such drugs remain unaware of their effects, hazards, and potency. Our study aimed to predict the PK of 34 fentanyl analogs and sought to establish a relationship between their structural characteristics and pharmacokinetic profiles. A hypothesis presented in the research suggests that pharmacokinetic parameters increase as fentanyl analogs become more lipophilic. The pharmacokinetics of acetylfentanyl, butyrylfentanyl, cyclopropylfentanyl, and valerylfentanyl were evaluated in rats following s. c injection at a dosage of 300 μg/kg. The results indicated that acetylfentanyl exhibited the shortest carbon side chain along with the lowest T1/2 and Cmax (Canfield and Sprague, 2024). It has been proposed that the incorporation of the functional group 3-carbomethoxy into fentanyl analogs may have the potential to reduce the duration of action by modifying their pharmacokinetic properties. This is attributed to the fact that more hydrophilic groups tend to accumulate minimally, if at all, in adipose tissue and are rapidly excreted. Additionally, it may be due to the susceptibility of the carbomethoxy group to rapid hydrolysis by non-specific esterases (Vucković et al., 2009). Therefore, it can be concluded that the physicochemical property exerting a significant influence on the pharmacokinetics of fentanyl analogs is lipophilicity. In conjunction with the predictive results, we summarize the relationship between the structural characteristics of fentanyl analogs and their pharmacokinetic profiles.
• Influence of carbon side chain on the in vivo pharmacokinetic characteristics of fentanyls
1. Pharmacokinetic parameters, including T1/2 and Vss, exhibit an increase with longer carbon side chain lengths, while AUC demonstrates a decrease as carbon side chain length increases.
2. Functional groups present on the carbon side chain, such as methoxyacetyl, are associated with a reduction in both T1/2 and Vss. In contrast, the presence of furanyl, tetrahydrofuranyl, cyclopropyl, and phenyl groups tends to increase both T1/2 and Vss to varying extents. The order is tetrahydrofuranyl < cyclopropyl < furanyl < phenyl.
• Influence of the position 3, or four of the piperidine ring on the in vivo pharmacokinetic characteristics of fentanyls
1. Both 3-carbomethoxy and methoxymethyl substituents at the four-position of the piperidine ring are associated with a reduction in T1/2, whereas a methyl group at the three-position is linked to an increase in T1/2.
• Influence of the directly connected to the nitrogen on the vivo pharmacokinetic characteristics of fentanyls.
1. The introduction of hydroxyl group to the carbon chain that is directly bonded to the nitrogen results in a reduction of both the T1/2 and Vss. Conversely, the incorporation of methyl group leads to an enhancement in these two parameters.
2. The incorporation of functional groups such as thienyl, methyl and phenyl can decrease the T1/2 and Vss to varying extents. The order is phenyl < thienyl < methyl.
4.2.3 Distribution of fentanyl analogs in tissues and organs
The brain-blood ratio serves as a crucial metric for estimating the pharmacokinetics of the central nervous system (CNS). Among various methodologies, the brain-blood ratio is favored over more complex techniques such as in situ brain perfusion and microdialysis due to its simplicity and practicality. Compounds exhibiting a higher brain-blood ratio demonstrate an enhanced capacity to traverse the blood-brain barrier and exert effects on the CNS. For any given compound, a ratio that approaches or slightly exceeds one suggests that it can readily cross the blood-brain barrier, whereas a significantly elevated ratio indicates a strong likelihood of accumulation within brain tissue (Kulkarni et al., 2016). A study has demonstrated that the brain-blood ratio for p-fluorofentanyl were significantly higher than that for fentanyl in several critical regions of the brain. This result indicates that the heightened toxicity associated with p-fluorofentanyl may be attributed to its enhanced permeability across the blood-brain barrier and increased exposure within brain tissue (Canfield and Sprague, 2025). This observation aligns with our prediction that the brain-blood ratio for p-fluorofentanyl is greater than that of fentanyl. The blue bars in Figure 5B highlight compounds exhibiting higher brain-blood ratios compared to fentanyl, which, as analogs of fentanyl, may possess a greater potential for abuse and pose increased risks relative to fentanyl. In addition, the use of opioids may be associated with several cardiovascular disorders, including cardiac arrest, tachycardia, bradycardia, and palpitations (Dai et al., 2024). The elevated Cmax of fentanyl and its derivates observed in the heart, as depicted in Figure 5A, indicates a potential for cardiotoxicity.
4.3 Comparison with the existing PBPK reports on fentanyl analogs
High-throughput (HT)-PBPK modeling for fentanyl analogs is grounded in a growing body of domain-specific research, whose findings both support and constrain the present QSAR-PBPK framework.
Key insights from fentanyl PBPK literature align with our study’s design and results. Björkman et al. (Björkman, 2003) demonstrated that simplifying fentanyl’s PBPK model—while prioritizing core parameters (tissue/blood partition coefficient, Kp; plasma protein unbound fraction, Fup)—retains predictive accuracy, validating our focus on QSAR-derived Kp (a critical driver of tissue distribution). Notably, Björkman also highlighted that cross-species Kp extrapolation (rat-to-human) introduces ≥30% error for lipophilic opioids, which echoes our observation that QSAR-predicted human Kp reduced steady-state volume of distribution (Vss) error to <1.5-fold (vs. >3-fold for extrapolation, Section 3.2). Population-specific PBPK studies further contextualize our model’s scope. Alsmadi (Alsmadi, 2023) and Kovar et al. (Kovar et al., 2020) showed that age-dependent physiology (e.g., pediatric CYP3A4 activity, neonatal tissue perfusion) alters fentanyl PK predictions by 25%–40%. While our use of standardized adult parameters (70 kg, healthy physiology) enables reliable analog-to-analog comparisons, these studies confirm that population-specific adjustments would be required for clinical translation—a key limitation.
Literature on fentanyl PBPK also supports QSAR’s utility for parameterization. Ni et al. (Ni et al., 2024) used QSAR-predicted enzyme inhibition constants (Ki) to model fentanyl-ritonavir interactions, yielding AUC fold-changes within 1.2-fold of clinical data—reinforcing our use of QSAR for Kp and physicochemical properties (logP, solubility) when experimental data are scarce. In contrast, Shankaran et al. (Shankaran et al., 2013) showed non-intravenous fentanyl models require route-specific absorption parameters (e.g., nasal permeability) to avoid ≥2-fold error, justifying our focus on intravenous administration (minimizing absorption uncertainty).
For novel analogs, Canfield and Sprague (Canfield and Sprague, 2024) reported a strong correlation (R2 = 0.89) between carbon side-chain length and T1/2 for illicit fentanyls—consistent with our predictions (e.g., valerylfentanyl, C5: T1/2 = 7.587 h; acetylfentanyl, C2: T1/2 = 2.766 h, Table 4). However, their observation that analogs with novel functional groups (e.g., furanyl) exhibit unexpected tissue binding underscores QSAR limitations for structurally divergent compounds.
4.4 Limitations of research
This study does have certain limitations. Firstly, given that fentanyl analogs represent an emerging class of controlled substances, clinical data are exceedingly scarce and cannot be utilized to evaluate the accuracy of our developed model. While our study proposes novel approaches for model assessment, the incorporation of additional clinical data would significantly enhance the validity of our model results. In addition, this study has established a PBPK model exclusively for fentanyl and its analogs following intravenous administration. However, other routes of administration, such as oral and nasal inhalation, have yet to be explored. To further investigate the pharmacokinetics of these substances, it is essential to develop additional routes of administration and conduct models after multiple administrations.
5 Conclusion
After a thorough evaluation of the modeling approach, QSAR-based PBPK models for 34 fentanyl analogs were developed. These models provide quantitative estimates of the pharmacokinetics of 34 fentanyl analogs in humans, with predictions for clinically validated analogs (e.g., sufentanil, alfentanil) aligning with available clinical data. These estimates help address critical data gaps in fentanyl analog hazard assessment and support preliminary prioritization of analogs for further experimental validation. It provided valuable insights that can not only guide further in vivo and in vitro experiments but also facilitate a preliminary assessment of their potential for abuse. This significantly addresses the existing data gap in this area.
Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: https://data.mendeley.com/datasets/92jxfc48gg/1.
Ethics statement
The animal study was approved by the Institutional Welfare and Ethics Committee for Laboratory Animals, Key Laboratory of Drug Monitoring and Control. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
SZ: Writing – original draft, Software, Visualization, Validation, Formal Analysis. YX: Writing – original draft, Data curation. XZ: Methodology, Writing – review and editing, Software. JR: Writing – review and editing, Project administration, Visualization. YC: Visualization, Writing – review and editing, Project administration. LK: Writing – review and editing, Formal Analysis, Data curation. KL: Formal Analysis, Data curation, Writing – review and editing. PX: Methodology, Writing – review and editing, Conceptualization, Supervision. FY: Writing – review and editing, Conceptualization, Supervision. DW: Project administration, Supervision, Writing – review and editing, Funding acquisition.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This project was supported by National Natural Science Foundation of China, No. 82130056.
Conflict of interest
The authors declare that they have no known economic interests or personal relationships that could potentially influence the work reported in this article.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
Alsmadi, M. M. (2023). Evaluating the pharmacokinetics of fentanyl in the brain extracellular Fluid, Saliva, Urine, and plasma of Newborns from Transplacental exposure from Parturient Mothers dosed with Epidural fentanyl utilizing PBPK modeling. Eur. J. Drug Metab. Pharmacokinet. 48, 567–586. doi:10.1007/s13318-023-00842-8
Armenian, P., Vo, K. T., Barr-Walker, J., and Lynch, K. L. (2018). Fentanyl, fentanyl analogs and novel synthetic opioids: a comprehensive review. Neuropharmacology 134, 121–132. doi:10.1016/j.neuropharm.2017.10.016
Björkman, S. (2003). Reduction and lumping of physiologically based pharmacokinetic models: prediction of the disposition of fentanyl and pethidine in humans by successively simplified models. J. Pharmacokinet. Pharmacodyn. 30, 285–307. doi:10.1023/a:1026194618660
Björkman, S., Stanski, D. R., Verotta, D., and Harashima, H. (1990). Comparative tissue concentration profiles of fentanyl and alfentanil in humans predicted from tissue/blood partition data obtained in rats. Anesthesiology 72, 865–873. doi:10.1097/00000542-199005000-00017
Bovill, J. G., Sebel, P. S., Blackburn, C. L., Oei-Lim, V., and Heykants, J. J. (1984). The pharmacokinetics of sufentanil in surgical patients. Anesthesiology 61, 502–506. doi:10.1097/00000542-198411000-00004
Bower, S., and Hull, C. J. (1982). Comparative pharmacokinetics of fentanyl and alfentanil. Br. J. Anaesth. 54, 871–877. doi:10.1093/bja/54.8.871
Bravo-Gómez, M. E., Camacho-GarcíA, L. N., Castillo-AlaníS, L. A., Mendoza-MeléNDEZ, M., and Quijano-Mateos, A. (2019). Revisiting a physiologically based pharmacokinetic model for cocaine with a forensic scope. Toxicol. Res. (Camb) 8, 432–446. doi:10.1039/c8tx00309b
Canfield, J. R., and Sprague, J. E. (2024). Influence of carbon side chain length on the in vivo pharmacokinetic and pharmacodynamic characteristics of illicitly manufactured fentanyls. Drug Test. Anal. 16, 1113–1121. doi:10.1002/dta.3636
Canfield, J. R., and Sprague, J. E. (2025). In vivo pharmacokinetic, pharmacodynamic and brain concentration comparison of fentanyl and para-fluorofentanyl in rats. Arch. Toxicol. 99, 287–297. doi:10.1007/s00204-024-03887-z
Dai, M., Dou, X., Chen, M., Yang, J., Long, J., and Lin, Y. (2024). Strong opioids-induced cardiac, neurologic, and respiratory disorders: a real-world study from 2004 to 2023 based on FAERS. Naunyn Schmiedeb. Arch. Pharmacol. 397, 4105–4121. doi:10.1007/s00210-023-02844-4
Ellison, C. A. (2018). Structural and functional pharmacokinetic analogs for physiologically based pharmacokinetic (PBPK) model evaluation. Regul. Toxicol. Pharmacol. 99, 61–77. doi:10.1016/j.yrtph.2018.09.008
EPA (2006). Approaches for the application of physiologically-based pharmacokinetic (PBPK) models and supporting data in risk assessment. Available online at: https://www.epa.ie.
Floresta, G., Rescifina, A., and Abbate, V. (2019). Structure-based approach for the prediction of Mu-opioid binding affinity of Unclassified designer fentanyl-Like molecules. Int. J. Mol. Sci. 20, 2311. doi:10.3390/ijms20092311
Han, J., Zhang, Z., Liu, X., Yang, H., and Liu, L. (2025). Prediction of pharmacokinetics for CYP3A4-metabolized drugs in pediatrics and Geriatrics using Dynamic age-dependent physiologically based pharmacokinetic models. Pharmaceutics 17, 214. doi:10.3390/pharmaceutics17020214
Helander, A., BäCKBERG, M., Signell, P., and Beck, O. (2017). Intoxications involving acrylfentanyl and other novel designer fentanyls - results from the Swedish STRIDA project. Clin. Toxicol. (Phila) 55, 589–599. doi:10.1080/15563650.2017.1303141
Kovar, L., Weber, A., Zemlin, M., Kohl, Y., Bals, R., Meibohm, B., et al. (2020). Physiologically-based pharmacokinetic (PBPK) modeling providing insights into fentanyl pharmacokinetics in adults and pediatric patients. Pharmaceutics 12, 908. doi:10.3390/pharmaceutics12100908
Kulkarni, A. D., Patel, H. M., Surana, S. J., Belgamwar, V. S., and Pardeshi, C. V. (2016). Brain-blood ratio: implications in brain drug delivery. Expert Opin. Drug Deliv. 13, 85–92. doi:10.1517/17425247.2016.1092519
Lim, C. B., Schug, S. A., Sunderland, V. B., Paech, M. J., and Liu, Y. (2012). A phase I pharmacokinetic and bioavailability study of a sublingual fentanyl wafer in healthy volunteers. Anesth. Analg. 115, 554–559. doi:10.1213/ANE.0b013e3182575cbf
Miller, N. A., Reddy, M. B., Heikkinen, A. T., Lukacova, V., and Parrott, N. (2019). Physiologically based pharmacokinetic Modelling for first-in-human predictions: an Updated model Building strategy illustrated with challenging Industry Case studies. Clin. Pharmacokinet. 58, 727–746. doi:10.1007/s40262-019-00741-9
Nair, A. B., and Jacob, S. (2016). A simple practice guide for dose conversion between animals and human. J. Basic Clin. Pharm. 7, 27–31. doi:10.4103/0976-0105.177703
Ni, L., Cao, Z., Jiang, J., Zhang, W., Hu, W., Zhang, Q., et al. (2024). Evaluating drug interactions between ritonavir and opioid analgesics: implications from physiologically based pharmacokinetic simulation. Pharm. (Basel) 17, 640. doi:10.3390/ph17050640
Patocka, J., Wu, W., Oleksak, P., Jelinkova, R., Nepovimova, E., Spicanova, L., et al. (2024). Fentanyl and its derivatives: pain-killers or man-killers? Heliyon 10, e28795. doi:10.1016/j.heliyon.2024.e28795
Pitsiu, M., Wilmer, A., Bodenham, A., Breen, D., Bach, V., Bonde, J., et al. (2004). Pharmacokinetics of remifentanil and its major metabolite, remifentanil acid, in ICU patients with renal impairment. Br. J. Anaesth. 92, 493–503. doi:10.1093/bja/aeh086
Rodgers, T., Leahy, D., and Rowland, M. (2005). Physiologically based pharmacokinetic modeling 1: predicting the tissue distribution of moderate-to-strong bases. J. Pharm. Sci. 94, 1259–1276. doi:10.1002/jps.20322
Sager, J. E., Yu, J., Ragueneau-Majlessi, I., and Isoherranen, N. (2015). Physiologically based pharmacokinetic (PBPK) modeling and simulation approaches: a systematic review of published models, applications, and model Verification. Drug Metab. Dispos. 43, 1823–1837. doi:10.1124/dmd.115.065920
Shankaran, H., Adeshina, F., and Teeguarden, J. G. (2013). Physiologically-based pharmacokinetic model for fentanyl in support of the development of Provisional Advisory levels. Toxicol. Appl. Pharmacol. 273, 464–476. doi:10.1016/j.taap.2013.05.024
UNODC (2024). The 2024 world drug report. Available online at: https://www.unodc.org/unodc/.
VAN DE Donk, T., Ward, S., Langford, R., and Dahan, A. (2018). Pharmacokinetics and pharmacodynamics of sublingual sufentanil for postoperative pain management. Anaesthesia 73, 231–237. doi:10.1111/anae.14132
Vucković, S., Prostran, M., Ivanović, M., Dosen-Mićović, L., Todorović, Z., Nesić, Z., et al. (2009). Fentanyl analogs: structure-activity-relationship study. Curr. Med. Chem. 16, 2468–2474. doi:10.2174/092986709788682074
Westmoreland, C. L., Hoke, J. F., Sebel, P. S., Hug, C. C., and Muir, K. T. (1993). Pharmacokinetics of remifentanil (GI87084B) and its major metabolite (GI90291) in patients undergoing elective inpatient surgery. Anesthesiology 79, 893–903. doi:10.1097/00000542-199311000-00005
Yau, E., Olivares-Morales, A., Gertz, M., Parrott, N., Darwich, A. S., Aarons, L., et al. (2020). Global Sensitivity analysis of the Rodgers and Rowland model for prediction of tissue: plasma partitioning coefficients: assessment of the key physiological and physicochemical factors that determine small-molecule tissue distribution. Aaps J. 22, 41. doi:10.1208/s12248-020-0418-7
Yk, R. Z. W., and Xd, X. (2024). Changes in pharmacokinetics of single dose of fentanyl in simulated high altitude in rats. J. Army Med. Univ. 46 (07), 732–737. doi:10.16016/j.2097-0927.202310018
Keywords: fentanyl analogs, PBPK model, human pharmacokinetics, QSAR model, abuse risk
Citation: Zhang S, Xu Y, Zeng X, Ran J, Chen Y, Kuai L, Li K, Xu P, Yan F and Wang D (2025) QSAR-based physiologically based pharmacokinetic (PBPK) modeling for 34 fentanyl analogs: model validation, human pharmacokinetic prediction and abuse risk insights. Front. Pharmacol. 16:1692293. doi: 10.3389/fphar.2025.1692293
Received: 25 August 2025; Accepted: 22 September 2025;
Published: 17 October 2025.
Edited by:
André Dallmann, Bayer, FranceReviewed by:
Damaris Albores-Garcia, Department of Pharmacobiology - CINVESTAV, MexicoRené Geci, University Hospital RWTH Aachen, Germany
Copyright © 2025 Zhang, Xu, Zeng, Ran, Chen, Kuai, Li, Xu, Yan and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fang Yan, MTAyMDA4MTg3OUBjcHUuZWR1LmNu; Dan Wang, ZG9uZ3d1ZmFuZzQ1NkAxNjMuY29t
†These authors have contributed equally to this work and share first authorship