Methods

Study design and ethics

We performed a retrospective cohort study using large-scale, longitudinal, individual-level healthcare data from the Epic Cosmos,20 a U.S. national clinical dataset created and managed by the Epic Systems Corporation that serves as the EHR vendor. The Cosmos integrates both inpatient and outpatient charts into a single, comprehensive record, including as patients move between health systems. It includes over 300 million patients from 1762 hospitals and 40,700 clinics from all 50 U.S. states, and is broadly representative of the U.S. Census for age, race, ethnicity, coverage, and the Social Vulnerability Index (SVI).20

We identified an initial cohort of 2,350,441 patients with at least one fentanyl UDT record from Jan 1, 2015 to Dec 31, 2023. Of these, about 295,728 patients met the study inclusion criteria (Fig. 1b). The data utilized for this study included comprehensive patient-level information, including demographics, medication records, diagnosis records, and laboratory results from inpatient, outpatient, and emergency department encounters.

Fig. 1: Study design and demographic characteristics of medical- and illicit-source fentanyl exposure cohorts. A. identification strategy of the source of fentanyl exposure; B. a flowchart of patient inclusion and exclusion; C. age distributions at the index date; and D. geographic distributions for the medical- and illicit-source fentanyl exposure cohorts. The index date was defined as the initial fentanyl exposure, either the date of the first fentanyl medication record or the first positive UDT result, whichever occurred first. UDT = Urine drug testing, Medical = medical-source fentanyl exposure cohort, Illicit = illicit-source fentanyl exposure cohort.

This study used secondary, de-identified data and was therefore exempt from ethics committee approval. The requirement for written informed consent was waived owing to the retrospective study design and the use of de-identified secondary data.

Key variables

Urine drug testing (UDT) data, a key diagnostic tool in addiction medicine,15 was used to identify and track fentanyl exposure. All fentanyl UDT records were extracted, including fentanyl or norfentanyl test results by screen or confirmatory methods. Records were removed if the specimen collection date (or the best available proxy date) was missing or if both the quantitative results with units and qualitative results (i.e., positive or negative, or abnormal) were unavailable. Only records that met all completeness requirements were retained for analysis. All quantitative measurements were standardized into a consistent unit (ng/mL) and classified using a predefined threshold of 3 ng/mL,21 with values at or above this threshold classified as positive for fentanyl exposure and lower values classified as negative.

To infer medical-source and potential illicit-source (non-medical) fentanyl exposure, UDT results were linked to documented fentanyl medication records using a sequential time-window screening method (Fig. 1a). We identified fentanyl medication records whose service date and days’ supply (when available) overlapped a 72-h window preceding the positive UDT date, consistent with the typical 2–4 day urine detection period for opioids.22 Positive UDTs with a linked medication record were attributed to medical-source exposure, whereas positive UDTs without a linked record were categorized as potential illicit-source exposure, indicating fentanyl exposure outside healthcare settings. The index date was defined as the initial fentanyl exposure, either the date of the first fentanyl medication record or the first positive UDT result, whichever occurred first. Fentanyl medication records included medication orders and medication administration records (Supplementary Table S3).

Participants

The study cohort included all individuals with at least one positive fentanyl UDT result between 2015 and 2023. They were classified into two subcohorts based on the inferred source of fentanyl exposure; individuals for whom all positive UDT results were temporally linked to fentanyl medication records were categorized as the medical-source fentanyl exposure cohort, whereas individuals for whom any positive UDT result lacked a temporally linked record were classified as the illicit-source fentanyl exposure cohort. The illicit-source cohort may therefore include individuals with exclusively illicit fentanyl exposure as well as those with both medical- and illicit-source exposures over time, including transitions between exposure sources.

Cohort characteristics

Baseline cohort characteristics were examined to assess their association with the source of fentanyl exposure. Demographic variables included age on the index date, sex, marital status, race/ethnicity, and state of residence. States were grouped into four major U.S. Census regions—Northeast, Midwest, South, and West—for all stratified analyses (Supplementary Table S1). Social context was also evaluated using the Social Vulnerability Index (SVI), which is based on percentile rankings by the ZIP code and ranges from 0 to 1, with higher values indicating greater vulnerability, and serves as a standardized, multidimensional measure of area-level social context. Fentanyl-related characteristics included the number of fentanyl medication records, the number of positive UDTs, and the mean of quantitative UDT values (Avg UDT). To ensure data reliability, the top 1% of UDT values were excluded from summary statistics to minimize the influence of potential outliers.

Baseline clinical conditions were assessed during the three months before the index date (the initial fentanyl exposure). A total of 40 clinical conditions were examined, including 24 chronic diseases (e.g., cancer, diabetes, and hypertension), five mental health disorders (e.g., anxiety, bipolar, and depression), three opioid-related harmful outcomes (overdose, abuse, and dependence), and eight substance use disorders (e.g., stimulant, alcohol, cannabis, and tobacco), all defined using International Classification of Diseases (ICD) codes (Supplementary Table S2). Opioid overdose was defined as nonfatal events using all ICD-10 subcodes under T40.0–T40.4 and T40.6, excluding adverse effect and underdosing codes. Fatal overdose events were not captured, as cause-of-death information (ICD: X40–44, X60–64, X85, and Y10–14) is not available in the dataset. Clinical conditions were identified using Epic Cosmos terminology-based mappings to ensure consistent condition definitions throughout the study period.

Statistical analysis

We compared baseline characteristics between the medical- and illicit-source cohorts and across region-stratified subgroups within each cohort. Variables included demographics, SVI measures, fentanyl-related factors, and the prevalence of the baseline clinical conditions. Most continuous and count variables were non-normally distributed; therefore, non-parametric tests were prespecified. We used two-sided Mann–Whitney U tests for pairwise comparisons and Kruskal–Wallis tests for regional comparisons. All tests were considered statistically significant at P < 0.05.

To evaluate the causal association between illicit-source fentanyl initiation and 30-day opioid-related harmful outcomes, including nonfatal opioid overdose, opioid abuse, and opioid dependence, we estimated hazard ratios (HRs) using Cox proportional hazards models.23 The exposure was illicit-source fentanyl initiation at the index date (the initial fentanyl exposure date). Individuals with a positive UDT at the index date were classified as having illicit-source fentanyl initiation (exposure), whereas those with a fentanyl medication record at the index date were classified as having medical-source fentanyl initiation (no exposure). The outcome was the occurrence of an opioid-related harmful outcome (e.g., opioid dependence) within 30 days after the index date. Individuals who experienced the outcome were coded as events at the observed time-to-event, whereas those without the outcome were coded as censored at day 30. Individuals with any prior opioid-related harmful outcomes were excluded for the association analysis to ensure assessment of incident outcomes.

Confounding was adjusted using stabilized inverse probability weighting (IPW) based on demographics, SVIs, and 37 baseline conditions (excluding opioid-related outcomes) measured before the index date. Other opioid prescriptions were not included as covariates, as the primary exposure of interest was the source of fentanyl exposure. Propensity scores were estimated by logistic regression of exposure on these covariates, with 3-fold cross-validation for training and evaluation. Post-weighting covariate balance was assessed using absolute standardized mean differences (SMD), a commonly used metric for evaluating between-group differences, with |SMD| ≤0.10 indicating adequate balance. Analyses were considered balanced if the proportion of unbalanced covariates was ≤ 2% of all covariates.24

An IPW-weighted Cox model was fitted in the entire study cohort to estimate adjusted HRs for illicit-source fentanyl initiation, after confirming covariate balance. We additionally conducted analyses restricted to the illicit-source cohort alone to further evaluate associations within a higher-risk subpopulation; all patients in the medical-source cohort were classified as having non-illicit-source fentanyl initiation.

Post-hoc subgroup analyses were conducted to examine whether associations varied by underlying clinical conditions, including substance use disorders, mental health conditions, and chronic conditions (ICD-based definitions in Supplementary Table S2). For each clinical category, patients were grouped by the presence or absence of any prior condition, and analyses were performed within each subgroup using the same analytic framework as in the primary analysis. Post-hoc sensitivity analyses evaluated alternative UDT screening windows (48, 72, and 96 h) and follow-up periods (30, 60, 90, and 180 days) to assess sensitivity to exposure and outcome definitions (Supplementary Tables S9 and S10). All statistical analyses were conducted using Python (version 3.12), SciPy, and Causallib25 libraries.

Temporal and regional patterns

We examined longitudinal patterns in disease prevalence within the medical- and illicit-source fentanyl cohorts from 2015 to 2023 to enable year-over-year comparisons of disease burden and its evolution over time. Region-stratified analyses were additionally performed to characterize geographic variation in these patterns. For each year, we calculated the percentage of patients diagnosed with each outcome of interest. To ensure temporal validity, we included only individuals whose index date occurred within the corresponding year, thereby assessing disease prevalence among individuals with confirmed fentanyl initiation during that year.

Role of the funding source

The funder had no role in the study design, data collection, data analysis, data interpretation, or writing of the manuscript.