Results

Descriptive summary of the study cohort

From an initial cohort of 300 million patients, we developed a study cohort of 295,728 patients from January 1, 2015, to December 31, 2023 (Table 1), including patients from the South (33.3%), Northeast (28.6%), Midwest (21.4%), and West (16.6%). The cohort selection pipeline is summarized in Fig. 1b. Of these, 85,535 (28.9%) patients were the medical-source fentanyl exposure cohort, while 210,193 (71.1%) patients were included in the illicit-source cohort.

Variables | Total | Medical-source | Illicit-source | P valuea
Demographics |  |  |  | 
Patients, No. (%) | 295,728 | 85,535 (28.9%) | 210,193 (71.1%) | NA
Male, No. (%) | 176,107 (59.6%) | 51,163 (59.8%) | 124,944 (59.4%) | 0.061
Female, No. (%) | 119,621 (40.4%) | 34,372 (40.2%) | 85,249 (40.6%) | 0.061
Married, No. (%) | 52,419 (17.7%) | 22,458 (26.3%) | 29,961 (14.3%) | <0.001
Age, mean (SD) | 43.229 (16.32) | 46.079 (18.39) | 41.99 (15.168) | <0.001
Age, median (IQR) | 41.0 (31.0, 55.0) | 46.0 (31.0, 60.0) | 39.0 (31.0, 53.0) | NA
Race, No. (%) |  |  |  | 
White | 181,228 (61.3%) | 51,283 (60.0%) | 129,945 (61.8%) | NA
Black or African American | 65,776 (22.2%) | 20,393 (23.8%) | 45,383 (21.6%) | NA
American Indian or Alaska native | 8609 (2.9%) | 2026 (2.4%) | 6583 (3.1%) | NA
Asian | 3358 (1.1%) | 1377 (1.6%) | 1981 (0.9%) | NA
Native Hawaiian or other Pacific Islander | 974 (0.3%) | 268 (0.3%) | 706 (0.3%) | NA
Other race | 14,518 (4.9%) | 4073 (4.8%) | 10,445 (5.0%) | NA
SVI, mean (SD) |  |  |  | 
Overall | 0.681 (0.271) | 0.677 (0.272) | 0.683 (0.27) | <0.001
Socioeconomic | 0.635 (0.288) | 0.635 (0.288) | 0.635 (0.288) | 0.523
Household characteristics | 0.612 (0.288) | 0.622 (0.282) | 0.608 (0.291) | <0.001
Housing type transportation | 0.69 (0.251) | 0.672 (0.256) | 0.697 (0.248) | <0.001
Racial ethnic minority | 0.692 (0.241) | 0.685 (0.246) | 0.695 (0.239) | <0.001
Fentanyl |  |  |  | 
Medication record count, mean (SD) | 14.85 (43.757) | 32.723 (63.352) | 7.577 (29.625) | <0.001
Urine drug testing (UDT) count, mean (SD) | 2.465 (6.045) | 1.405 (1.858) | 2.897 (7.026) | <0.001
Avg UDT value, mean (SD) | 77.661 (101.559) | 47.542 (61.187) | 87.428 (112.875) | <0.001
Clinical conditionsb |  |  |  | 
Cancer, No. (%) | 6059 (2.0%) | 2732 (3.2%) | 3327 (1.6%) | <0.001
Cardiovascular disease | 44,709 (15.1%) | 14,660 (17.1%) | 30,049 (14.3%) | <0.001
Hypertension | 48,309 (16.3%) | 15,775 (18.4%) | 32,534 (15.5%) | <0.001
Diabetes | 20,726 (7.0%) | 7381 (8.6%) | 13,345 (6.3%) | <0.001
Depression | 32,883 (11.1%) | 6260 (7.3%) | 26,623 (12.7%) | <0.001
Anxiety | 43,634 (14.8%) | 8287 (9.7%) | 35,347 (16.8%) | <0.001
Bipolar | 17,145 (5.8%) | 2402 (2.8%) | 14,743 (7.0%) | <0.001
Stimulant use disorder | 23,957 (8.1%) | 2573 (3.0%) | 21,384 (10.2%) | <0.001
Other psychoactive | 27,012 (9.1%) | 2431 (2.8%) | 24,581 (11.7%) | <0.001
Sedative use | 3408 (1.2%) | 383 (0.4%) | 3025 (1.4%) | <0.001
Other substance use disorder | 48,853 (16.5%) | 4887 (5.7%) | 43,966 (20.9%) | <0.001
Opioid overdose | 9176 (3.1%) | 525 (0.6%) | 8651 (4.1%) | <0.001
Opioid abuse | 18,574 (6.3%) | 1232 (1.4%) | 17,342 (8.3%) | <0.001
Opioid dependence | 34,605 (11.7%) | 1999 (2.3%) | 32,606 (15.5%) | <0.001

The illicit-source cohort was younger on average than the medical-source cohort (mean age, 42.0 versus 46.1 years; P < 0.001). A greater proportion of patients in the medical-source cohort were aged 50 years or older (43.9% versus 30.2%; Supplementary Fig. S1d). The illicit-source cohort was less commonly married (14.3% versus 26.3%; P < 0.001) and had a slightly higher SVI overall score (0.683 versus 0.677; P < 0.001). For state-level distributions (Fig. 1d), medical-source exposure was most common in North Carolina (10.6%), Ohio (10.5%), and California (10.3%), whereas illicit-source exposure was highly concentrated in Ohio (12.2%), followed by Massachusetts (10.5%) and California (8.8%).

Notably, patients in the medical-source cohort had a substantially higher mean number of fentanyl medication records (32.7 versus 7.58; P < 0.001) and lower UDT values (47.5 ng/mL versus 87.4 ng/mL; P < 0.001). Chronic diseases, including cancer, cardiovascular disease, hypertension, and diabetes, were more prevalent in the medical-source cohort, with cancer rates two-fold higher. In contrast, mental health disorders, substance use disorders, and opioid-related harmful outcomes, such as opioid overdose (nonfatal), abuse, and dependence, were more common in the illicit-source cohort. Rates were more than three-fold higher for stimulant use disorder and six-, five-, and six-fold higher for overdose, abuse, and dependence, respectively, compared with the medical-source cohort. All clinical conditions are presented in Supplementary Table S4, and region-stratified statistics for the medical-source cohort are provided in Supplementary Tables S6 and S7.

Region-stratified subgroups in the illicit-source fentanyl exposure cohort

Among the 210,193 patients with inferred illicit-source fentanyl exposure (Table 2), the largest proportion resided in the Northeast (32.2%), followed by the South (30.2%), the Midwest (21.3%), and the West (16.1%). The Midwest had the highest mean number of fentanyl medication records (9.59) and the highest Avg UDT (136 ng/mL). The West had the lowest number of medication records (5.74) but the second highest Avg UDT (104 ng/mL). The Northeast showed relatively higher UDT testing frequency (mean count 4.87).

Variables | Total | West | Midwest | South | Northeast | P valuea
Demographics |  |  |  |  |  | 
Patients, No. (%) | 210,193 | 33,760 (16.1%) | 44,857 (21.3%) | 63,475 (30.2%) | 67,775 (32.2%) | NA
Male, No. (%) | 124,944 (59.4%) | 20,699 (61.3%) | 25,228 (56.2%) | 36,794 (58.0%) | 42,004 (62.0%) | <0.001
Female, No. (%) | 85,249 (40.6%) | 13,061 (38.7%) | 19,629 (43.8%) | 26,681 (42.0%) | 25,771 (38.0%) | <0.001
Married, No. (%) | 29,961 (14.3%) | 4946 (14.7%) | 6405 (14.3%) | 10,420 (16.4%) | 8182 (12.1%) | <0.001
Age, mean (SD) | 41.99 (15.168) | 41.688 (16.115) | 40.47 (14.597) | 42.37 (15.698) | 42.888 (14.353) | <0.001
Age, median (IQR) | 39.0 (31.0, 53.0) | 38.0 (30.0, 53.0) | 38.0 (30.0, 50.0) | 40.0 (31.0, 53.0) | 41.0 (32.0, 53.0) | NA
Race, No. (%) |  |  |  |  |  | 
White | 129,945 (61.8%) | 20,721 (61.4%) | 28,751 (64.1%) | 36,461 (57.4%) | 43,918 (64.8%) | NA
Black or African American | 45,383 (21.6%) | 3832 (11.4%) | 10,614 (23.7%) | 18,959 (29.9%) | 11,944 (17.6%) | NA
AIAN | 6583 (3.1%) | 2744 (8.1%) | 1721 (3.8%) | 683 (1.1%) | 1430 (2.1%) | NA
Asian | 1981 (0.9%) | 830 (2.5%) | 306 (0.7%) | 325 (0.5%) | 520 (0.8%) | NA
NHPI | 706 (0.3%) | 190 (0.6%) | 115 (0.3%) | 111 (0.2%) | 289 (0.4%) | NA
Other race | 10,445 (5.0%) | 3031 (9.0%) | 526 (1.2%) | 2329 (3.7%) | 4547 (6.7%) | NA
SVI, mean (SD) |  |  |  |  |  | 
Overall | 0.683 (0.27) | 0.722 (0.242) | 0.662 (0.279) | 0.726 (0.233) | 0.641 (0.297) | <0.001
Socioeconomic | 0.635 (0.288) | 0.65 (0.272) | 0.624 (0.301) | 0.698 (0.246) | 0.58 (0.308) | <0.001
Household characteristics | 0.608 (0.291) | 0.572 (0.292) | 0.621 (0.295) | 0.648 (0.268) | 0.58 (0.301) | <0.001
Housing type transportation | 0.697 (0.248) | 0.761 (0.235) | 0.677 (0.232) | 0.685 (0.248) | 0.693 (0.256) | <0.001
Racial ethnic minority | 0.695 (0.239) | 0.804 (0.162) | 0.618 (0.254) | 0.72 (0.219) | 0.672 (0.25) | <0.001
Fentanyl |  |  |  |  |  | 
Medication record count, mean (SD) | 7.577 (29.625) | 5.739 (24.085) | 9.585 (33.697) | 7.129 (23.789) | 7.615 (33.858) | <0.001
UDT count, mean (SD) | 2.897 (7.026) | 1.986 (2.632) | 2.186 (3.359) | 1.786 (2.921) | 4.868 (11.331) | <0.001
Avg UDT value, mean (SD) | 87.428 (112.875) | 104.517 (126.956) | 136.315 (163.416) | 63.979 (76.402) | 76.386 (103.092) | <0.001
Clinical conditionsb |  |  |  |  |  | 
Cancer, No. (%) | 3327 (1.6%) | 536 (1.6%) | 727 (1.6%) | 1032 (1.6%) | 1031 (1.5%) | 0.42
Cardiovascular disease | 30,049 (14.3%) | 4550 (13.5%) | 6578 (14.7%) | 9781 (15.4%) | 9124 (13.5%) | <0.001
Hypertension | 32,534 (15.5%) | 4305 (12.8%) | 7066 (15.8%) | 11,060 (17.4%) | 10,095 (14.9%) | <0.001
Diabetes | 13,345 (6.3%) | 1966 (5.8%) | 2736 (6.1%) | 4309 (6.8%) | 4333 (6.4%) | <0.001
Depression | 26,623 (12.7%) | 2942 (8.7%) | 6011 (13.4%) | 6688 (10.5%) | 10,976 (16.2%) | <0.001
Anxiety | 35,347 (16.8%) | 4144 (12.3%) | 7797 (17.4%) | 8713 (13.7%) | 14,683 (21.7%) | <0.001
Bipolar | 14,743 (7.0%) | 1541 (4.6%) | 3478 (7.8%) | 3891 (6.1%) | 5828 (8.6%) | <0.001
Stimulant use disorder | 21,384 (10.2%) | 3703 (11.0%) | 4369 (9.7%) | 5662 (8.9%) | 7629 (11.3%) | <0.001
Other psychoactive | 24,581 (11.7%) | 2395 (7.1%) | 5359 (11.9%) | 6369 (10.0%) | 10,444 (15.4%) | <0.001
Sedative use | 3025 (1.4%) | 318 (0.9%) | 564 (1.3%) | 826 (1.3%) | 1312 (1.9%) | <0.001
Other substance use disorder | 43,966 (20.9%) | 6110 (18.1%) | 9113 (20.3%) | 10,980 (17.3%) | 17,728 (26.2%) | <0.001
Opioid overdose | 8651 (4.1%) | 1096 (3.2%) | 2016 (4.5%) | 2556 (4.0%) | 2975 (4.4%) | <0.001
Opioid abuse | 17,342 (8.3%) | 2340 (6.9%) | 3714 (8.3%) | 3797 (6.0%) | 7475 (11.0%) | <0.001
Opioid dependence | 32,606 (15.5%) | 4262 (12.6%) | 6171 (13.8%) | 5456 (8.6%) | 16,699 (24.6%) | <0.001

The Northeast had the highest prevalence of mental health and substance use disorders and most opioid-related harmful outcomes, except for opioid overdose. In contrast, the South showed the lowest prevalence of opioid abuse (11.0% versus 6.0%; P < 0.001) and dependence (24.6% versus 8.6%; P < 0.001). All conditions are presented in Supplementary Table S5.

Impact of illicit-source fentanyl exposure

Illicit-source fentanyl initiation was associated with increased risk of 30-day opioid-related harmful outcomes, including nonfatal opioid overdose, opioid abuse, and opioid dependence. In the entire cohort (Table 3), adjusted HRs for illicit-source initiation were 2.99 (95% confidence interval [CI], 2.71–3.29; P < 0.001) for opioid overdose, 1.96 (95% CI, 1.84–2.10; P < 0.001) for abuse, and 3.05 (95% CI, 2.89–3.22; P < 0.001) for dependence, compared with medical-source initiation, after confounding adjustment.

Subgroup | # Patients | Overdose | Abuse | Dependence
All | 194,348 | 2.99 (2.71–3.29); <0.001 | 1.96 (1.84–2.10); <0.001 | 3.05 (2.89–3.22); <0.001
Substance use disorder |  |  |  | 
Positive | 77,499 | 3.00 (2.62–3.44); <0.001 | 2.07 (1.88–2.27); <0.001 | 3.32 (3.07–3.59); <0.001
Negative | 116,849 | 3.04 (2.66–3.48); <0.001 | 2.08 (1.89–2.29); <0.001 | 3.17 (2.94–3.41); <0.001
Mental health conditions |  |  |  | 
Positive | 64,029 | 2.80 (2.38–3.29); <0.001 | 2.49 (2.21–2.80); <0.001 | 3.26 (2.96–3.58); <0.001
Negative | 130,319 | 3.09 (2.75–3.48); <0.001 | 1.78 (1.64–1.92); <0.001 | 2.97 (2.78–3.17); <0.001
Chronic conditions |  |  |  | 
Positive | 117,891 | 3.54 (3.11–4.02); <0.001 | 2.47 (2.26–2.71); <0.001 | 3.50 (3.25–3.76); <0.001
Negative | 76,457 | 2.54 (2.18–2.96); <0.001 | 1.58 (1.43–1.75); <0.001 | 2.79 (2.57–3.03); <0.001

In analyses restricted to the illicit-source cohort alone, associations remained significant, with adjusted HRs of 2.84 (95% CI, 2.46–3.28; P < 0.001) for opioid overdose, 1.53 (95% CI, 1.40–1.68; P < 0.001) for abuse, and 2.25 (95% CI, 2.09–2.42; P < 0.001) for dependence.

Across all clinical subgroups in the entire cohort (Table 3), illicit-source initiation was consistently associated with increased risk of 30-day opioid-related harmful outcomes. For nonfatal opioid overdose, adjusted HRs were 3.00 versus 3.04 for patient subgroups with and without substance use disorders, respectively, and 2.80 versus 3.09 (mental health disorders), and 3.54 versus 2.54 (chronic conditions). The largest differences in HRs were observed for chronic conditions, with consistently higher HRs among patients with chronic conditions than among those without. Similar findings were observed in subgroup analyses restricted to the illicit-source cohort alone (Supplementary Table S8), in which illicit-source initiation remained consistently associated with increased risk of these outcomes across clinical subgroups.

Sensitivity analyses using alternative UDT screening windows demonstrated consistent findings (Supplementary Table S9). Increased risks were also observed across alternative follow-up outcome windows, with decreasing HRs as the prediction window lengthened (Supplementary Table S10). Covariate balance (SMDs) is presented in Supplementary Tables S11 and S12 and Supplementary Fig. S2, and details of the propensity score model are provided in Supplementary Tables S13 and S14 and Supplementary Fig. S3.

Temporal patterns in opioid-related harmful outcomes

From 2015 to 2023, nonfatal opioid overdose prevalence was consistently higher in the illicit-source fentanyl exposure cohort (Fig. 2a). Overdose prevalence in the illicit-source cohort increased markedly over time, reaching 18.6%, whereas only a modest increase was observed in the medical-source cohort, reaching 4.3%. For opioid dependence and abuse, the illicit-source cohort exhibited substantially higher prevalence and steeper year-over-year increases than the medical-source cohort until 2020. After 2020, prevalence in the illicit-source cohort declined, particularly for dependence, but remained consistently higher than in the medical-source cohort throughout the study period.

Fig. 2: Temporal patterns in opioid overdose (nonfatal), dependence, and abuse from 2015 to 2023. A. the medical- and illicit-source fentanyl exposure cohorts; and B. the four regions within medical-source (left) and illicit-source fentanyl exposure cohorts (right). Medical = medical-source fentanyl exposure cohort, Illicit = illicit-source fentanyl exposure cohort.

Region-stratified temporal patterns followed the overall cohort-level patterns (Fig. 2b). Within the medical-source cohort, regional trajectories were relatively stable and comparable over time. In contrast, the illicit-source cohort showed regional variation. The Northeast generally exhibited higher prevalence levels of opioid abuse and dependence, with sharp increases beginning in 2016 and peaking around 2019. Across all regions, prevalence remained consistently higher in the illicit-source cohort than in the medical-source cohort throughout the study period.