Urban scaling of opioid overdose deaths in the USA: a cross-sectional study in three periods between 2005 and 2017.

OBJECTIVES: To describe the association between population size, population growth and opioid overdose deaths-overall and by type of opioid-in US commuting zones (CZs) in three periods between 2005 and 2017. SETTINGS: 741 CZs covering the entirety of the US CZs are aggregations of counties based on commuting patterns that reflect local economies. PARTICIPANTS: We used mortality data at the county level from 2005 to 2017 from the National Center for Health Statistics. OUTCOME: Opioid overdose deaths were defined using underlying and contributory causes of death codes from the International Classification of Diseases, 10th revision (ICD-10). We used the underlying cause of death to identify all drug poisoning deaths. Contributory cause of death was used to classify opioid overdose deaths according to the three major types of opioid, that is, prescription opioids, heroin and synthetic opioids other than methadone.
RESULTS: Opioid overdose deaths were disproportionally higher in largely populated CZs. A CZ with 1.0% larger population had 1.10%, 1.10%, and 1.16% higher opioid death count in 2005-2009, 2010-2014, and 2015-2017, respectively. This pattern was largely driven by a high number of deaths involving heroin and synthetic opioids, particularly in 2015-2017. Population growth over time was associated with lower age-adjusted opioid overdose mortality rate: a 1.0% increase in population over time was associated with 1.4% (95% CI: -2.8% to 0.1%), 4.5% (95% CI: -5.8% to -3.2%), and 1.2% (95% CI: -4.2% to 1.8%) lower opioid overdose mortality in 2005-2009, 2010-2014, and 2015-2017, respectively. The association between positive population growth and lower opioid mortality rates was stronger in larger CZs.
CONCLUSIONS: Opioid overdose mortality in the USA was disproportionately higher in mid-sized and large CZs, particularly those affected by declines in population over time, regardless of the region where they are located. © Author(s) (or their employer(s)) 2022. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

In this study of patterns of opioid overdose mortality in US commuting zones (CZs), counts and rates of opioid overdose deaths were disproportionately higher in mid-sized and large CZs, compared with small CZs.

The phenomenon described here is known as superlinear scaling, that is, disproportionally more deaths in large cities. However, our analysis showed a more complex pattern of mortality distribution, as CZs that experienced population growth over time had smaller rates of opioid overdose mortality, compared with CZs that experienced stagnant growth.

The use of CZs as the geographical unit of analysis reflects commuting patterns and social interactions that may be at the core of the mechanisms underlying the associations between population metrics and opioid overdose outcomes.

This study does not aim to examine potential causal mechanisms between population metrics and opioid overdose deaths. We recognise the complexity of the determinants of substance use disorders and the exploratory nature of our work.

Nonetheless, we believe that our approach can be a valuable way to understand broad patterns of opioid mortality. In particular, using population growth as a potential proxy for other macro-level determinants can be helpful to predict which areas may be at a greater risk of increases in opioid overdose mortality.

Introduction

Drug overdose deaths increased exponentially in the USA over the past 40 years.1 Currently, opioids account for almost 75% of all overdose deaths in 2020.2 Opioid overdose mortality has grown steadily despite increased attention from policymakers, public health professionals and the general public.3 4 National trends in opioid overdose mortality show that the number of opioid analgesic overdose deaths declined slightly between 2011 and 2013, coinciding with interventions to limit prescribing and dispensing of opioid analgesics.5–8 After a short period, mortality began to rise again, a trend that was in part driven by heroin overdose deaths,4 followed by sharp increases in overdoses from illicitly manufactured synthetic opioids.9 10

Substantial evidence points to geographical patterning in the distribution of opioid overdose deaths. Studies initially identified hot-spots for overdose deaths clustered in areas of the country with high rates of long-term unemployment, such as Appalachia and the Rust Belt.1 11 But the demographic profile and geographical distribution of opioid overdose deaths have changed over time.10–13 Initial reports showed that mortality from prescription opioids was higher among older (vs younger) age groups, white (vs non-white) and rural (vs urban) communities.12 However, further studies found an increase in opioid overdose deaths in urban centres, particularly deaths due to heroin overdose in younger age groups.1 Data from 2017 show that opioid mortality in urban counties has surpassed that of rural counties.14 However, little is known about how opioid overdose deaths are patterned across areas of different population sizes. Additional characterisation of the opioid epidemic may shed light on underlying dynamics of opioid overdose mortality across the urban–rural spectrum.

We adopted the framework of urban scaling to study the association between population size, population growth and opioid overdose deaths. According to this framework, cities share a set of universal characteristics that lead to non-linear patterning of outcomes among city residents.15–17 Specifically, superlinear and sublinear scaling are two potential responses to changes in city size. Superlinear scaling is observed when large cities have a disproportionately large number of outcomes compared with small cities. For example, large cities have higher economic outputs compared with smaller cities, even after accounting for population size.15 In other words, the economic output per capita is higher in large (vs small) cities. Sublinear scaling is observed when cities have a disproportionately small number of outcomes compared with small cities. Examples of sublinear scaling are common when dealing with features of the infrastructure of cities, such as the length of road networks; large cities need a relatively short road network compared with small cities due to economies of scale.15

Health outcomes also show scaling behaviours.18–22 Sexually transmitted infections (ie, chlamydia, gonorrhoea and syphilis) in US cities scale superlinearly—disproportionally more cases of sexually transmitted diseases occur in larger cities even after adjusting for population,19 potentially due to a relatively larger number of contacts among residents of larger cities.20 Conversely, motor vehicle crashes scale sublinearly, potentially due to shorter road network length.21 Additionally, while most studies on scaling have used population size as the main exposure, recent reports have also examined population growth. For example, UK districts that have shrunk over time have higher age-standardised mortality,23 while Latin American cities with higher population growth have higher life expectancy.24 Population growth can be thought of as a cause of improved living conditions, such as when the share of the working-age population is larger than the non-working-age share of the population,25 or as a consequence of it, as increases in population can result from changes in economic opportunities, with economically strong areas attracting people from other regions.22 Growing areas may receive more migrants, who tend to have improved health statuses,26 causing in part these growth–mortality associations.27

This study examines the urban scaling of opioid overdose deaths in the USA. The aims of this study are: (1) to estimate the scaling parameter of opioid overdose deaths in US commuting zones (CZs) between 2005 and 2017; and (2) to estimate the independent associations of population size, growth and their interaction, with opioid overdose deaths. We hypothesised that a large number of contacts resulting from larger population counts will be associated with a disproportionally higher count of opioid deaths (superlinear scaling) due to an increase in successful matches between susceptible people and sources of prescription and illicit opioids.28 The relationship between population growth and opioid mortality may be inverted, as population growth can be an indicator of communities with a thriving economy and potentially lower rates of mortality.29–31

Methods

Study setting

We used data on all CZs of the USA from 2005 to 2017. CZs are aggregations of counties based on commuting patterns that were developed by the US Department of Agriculture32 (section 1 in the online supplemental material provides further details on CZ’s definition and relationship with county boundaries). Unlike counties, which reflect political boundaries, CZs reflect local economies, which may be important to the understanding of the macro-determinants of opioid outcomes. There are 741 CZs that cover the entirety of the USA and provide a complete picture of the country, from rural to highly urbanised areas. These CZs have a perfect overlap with county boundaries allowing for computation of CZ-level measures by aggregating data from the counties that are part of the CZs.32 We analysed data in three periods roughly corresponding to different waves of the opioid epidemic: 2005–2009 characterised by an increase in deaths due to prescription opioids, 2010–2014 marked by an increase in the number of heroin deaths and 2015–2017 marked by an increase in deaths associated with synthetic opioids, in particular illicitly manufactured fentanyl.33

Data sources

We used mortality data at the county level from 2005 to 2017 from the National Center for Health Statistics34 and population estimates from the US Census Bureau for the same period.35

Patient and public involvement

No patients or members of the public were involved in this study.

Outcome variable

The main outcomes of this study are opioid overdose deaths, overall and by type of opioid involved. Opioid overdose deaths were defined using underlying and contributory causes of death codes from the International Classification of Diseases, 10th revision (ICD-10). Consistent with the methodology used by the Centers for Disease Control and Prevention (CDC), we used the underlying cause of death to identify all drug poisoning deaths coded as accidental poisoning (X40–X44), self-poisoning (X60–X64), homicidal poisoning (X85) and those of unknown intent (Y10–Y14). Then, we used the contributory cause of death to determine if the overdose was associated with opioids (ICD-10 codes: T40.0, T40.1, T40.2, T40.3, T40.4 and T40.6).36 37 In addition to the overall number of opioid overdose deaths, we also examined overdose deaths according to the three major types of opioid involved in the overdose, that is, prescription opioids, heroin and synthetic opioids other than methadone, also following the CDC classification.37 Each death certificate may include more than one contributory code. Therefore, for the analysis by opioid type, the same death could have been included in more than one group. Approximately 25% of all drug poisoning deaths in the data analysed were coded as ‘unspecified’ (T50.9) with regard to the type of drug involved. To address this issue, we followed the approach used by Ruhm38 to assign unspecified drug codes to either opioid or non-opioid overdose group (section 2 in the online supplemental material provides further details on this procedure).

For the first aim, we operationalised the main outcome as the sum of all opioid overdose deaths in each CZ in three periods: 2005–2009, 2010–2014 and 2015–2017. We also calculated counts associated with each opioid type. For the second aim, we calculated age-standardised opioid overdose mortality rates, and by type, in the three periods. We used the age distribution of the US standard population in 2000 and the direct method of standardisation to calculate age-standardised rates per CZ per period.

Exposure variables

For the first aim, we used population size, operationalised as the natural logarithm (ln) of the average population in each of the three periods, as the main exposure. For the second aim, the exposures were population size (operationalised as in the first aim) and population growth, measured as the relative change in population between the first and last years of each period (eg, for the first period (2005–2009), the difference between the population in 2009 and 2005, divided by the population in 2005). We adjusted the models for region where the CZ is located (Northeast, Midwest, South and West).

Analytical approach

For descriptive purposes, we first examined mortality rates by groups of CZs classified into deciles of the population distribution in 2005.

For the first aim, we estimated the scaling parameter for opioid overdose deaths by regressing the log of death counts due to opioid overdose on the log of the average population in each period. We used the following equation:

Where Yi is the death count for the i-th CZ and Sizei is the population of the CZ. β is the scaling coefficient: β<1 corresponds to sublinear scaling and β>1 corresponds to superlinear scaling. A total of 36, 32, and 42 CZs did not have any opioid overdose death in 2005–2009, 2010–2014 and 2015–2017, respectively. We also calculated the scaling parameter deaths associated with three opioid types (prescription opioid, heroin and synthetic opioid). For the models by opioid type, we conducted a sensitivity analysis using a consistent subsample of CZs (those who had at least one death involving each opioid type).

To visually depict the relationship between opioid overdose deaths and population, we created three plots (one per period) with the log of population on the x axis and the log of opioid overdose counts on the y axis. We also mapped the residuals (ei above) to explore the geographical distribution of CZs with higher or lower mortality.

For the second aim, we used linear models to estimate the associations between population size, population growth and age-standardised opioid mortality rates, and test an interaction between population size and growth. We used the following equation:

Where Yi is the age-standardised opioid mortality rate for the i-th CZ, Size is the average population for the period, Growth is the relative growth during the period, and Region is a vector of indicator variables for the census region of the CZ. We also examined separate models for prescription opioid, heroin and synthetic opioid mortality.

To assess whether differential quality of coding by area may bias our results, we tested for an association between the proportion of records coded as unspecified drug and population size. This analysis showed that population size was not associated with counts of records coded as unspecified, except for the last period (2015–2017), in which larger CZs were more likely to have fewer counts of records with unspecified codes (online supplemental figure 1). The percentage of unspecified codes also varied across regions, which is likely to be related to differences in drug profile. In the first period, the median per cent varied from about 19% in the Northeast to 25% in the South. In the third period, the median per cent unspecified varied from 3% in the Northeast to 14% in the Midwest (see online supplemental table 1 for median and percentile variation). Finally, we compared the results from the imputed data with the results from the original data to check for potential biases resulting from the imputation.

Detailed results including mortality due to various combinations of opioids, visualisations of the relationships between population metrics and overdose deaths, and per cent of unspecified poisoning deaths by CZs were included in an interactive app available here: https://drexel-uhc.shinyapps.io/Opioid_Scaling/. All analyses were conducted in STATA V.14.39

Results

Figure 1 shows boxplots for age-standardised opioid overdose mortality among CZs grouped by decile of the population in 2005–2009, 2010–2014 and 2015–2017. Median rates were generally higher for CZs in the top deciles of the population compared with smaller CZs. Median rates increased over time, particularly between the second and third period and in deciles 9 and 10, corresponding to CZs with more than 400 000 people.

Figure 1

Boxplots for age-standardised opioid mortality rate in groups of CZs by decile of the population in three periods. Population decile 1 has 75 CZs, 90% of which have less than 10 000 people. The other deciles have 74 CZs each. CZs, commuting zones; NCHS, National Center for Health Statistics.

Opioid overdose deaths scaled superlinearly in the three periods; opioid overdose deaths were disproportionately higher in larger CZs as compared with smaller CZs (table 1 and figure 2). A CZ with 1.00% larger population had 1.10%, 1.10%, and 1.16% higher opioid death count in 2005–2009, 2010–2014, and 2015–2017, respectively (table 1). These coefficients were within the range of those obtained from the original data (ie, before imputation of unspecified drug codes) (online supplemental table 2). Death counts by opioid type show that the superlinear behaviour of opioid overdose deaths was strongly driven by deaths involving heroin in the first and second periods, and deaths involving heroin and synthetic opioids in the third period. Prescription opioid deaths scaled superlinearly but the scaling coefficients were very close to 1, while heroin overdose deaths scaled superlinearly in all periods. Scaling of synthetic opioids changed from sublinear in the first two periods to superlinear in the last period. Figures by opioid types are included in the interactive app. Coefficients obtained from sensitivity analysis using only CZs with at least one opioid death showed similar patterns (online supplemental table 3).

Table 1

Scaling coefficients for opioid deaths and subgroups of opioid in three periods

		2005–2009	2010–2014	2015–2017
All opioids	Scaling coefficient	1.10	1.10	1.16
	95% CI	1.07 to 1.13	1.07 to 1.13	1.13 to 1.19
Prescription opioid	Scaling coefficient	1.03	1.03	1.03
	95% CI	0.99 to 1.07	1.00 to 1.07	1 to 1.06
Heroin	Scaling coefficient	1.07	1.15	1.20
	95% CI	1.00 to 1.13	1.10 to 1.2	1.15 to 1.25
Synthetic opioid	Scaling coefficient	0.88	0.92	1.16
	95% CI	0.84 to 0.91	0.89 to 0.95	1.10 to 1.21

For the analysis by type of opioid, that is, prescription opioid, heroin and synthetic opioid, more than one opioid can be reported in the death certificate. The number of commuting zones included in the models varies by period and by opioid type as commuting zones with zero deaths for any given outcome are dropped from the model when the outcome is log-transformed. A sensitivity analysis, including only commuting zones with at least one death involving each opioid type, resulted in the similar superlinear/sublinear scaling pattern (online supplemental table 3).

Figure 2

Opioid deaths in US commuting zones (CZs) in three periods: 2005–2009, 2010–2014 and 2015–2017. The line represents the regression line of the association between opioid overdose death counts and population counts. Red CZs are those with positive residuals (death counts are larger than predicted by the model) and green CZs are those with negative residuals (death counts are smaller than predicted).

We also mapped the residuals from the scaling regression, finding that heroin deaths were concentrated in Rust Belt and Midatlantic CZs, prescription opioids deaths were concentrated in the Appalachia and some areas in the West, and synthetic opioid deaths were concentrated in the Appalachia, Rust Belt, Midatlantic, and New England (see interactive app).

Table 2 shows the associations of CZ population size, population growth and their interaction with age-standardised opioid overdose mortality rates. There was evidence of a negative interaction between population size, that is, average population, and population growth, that is, change in population over time. For CZs that did not grow in population (stable), greater population size was associated with higher opioid mortality rates: a 1.0% higher population size was associated with an increase in opioid overdose mortality rate of 0.12% (95% CI 0.08% to 0.17%) in 2005–2009, 0.14% (95% CI 0.10% to 0.18%) in 2010–2014 and 0.16% (95% CI 0.12% to 0.20%) in 2015–2017. However, the association became weaker for CZs where the population was growing over time: for every 1.0% increase in population growth, that is, change in population over time, the association between overdose mortality and population size across CZs declined (moved closer to the null) by 0.01%, 0.01%, and 0.02% in 2005–2009, 2010–2014, and 2015–2017, respectively.

Table 2

Associations between opioid overdose mortality rates per 100 000 (log 2), population growth and size in three periods

	2005–2009	2010–2014	2015–2017
Growth (% change for a 1% increase in population for a CZ with population 100 000)	−1.4 (−2.8 to 0.1)	−4.5 (−5.8 to –3.2)	−1.2 (−4.2 to 1.8)
Size (% change for a 1% change in average population for a stable CZ)	0.12 (0.08 to 0.17)	0.14 (0.10 to 0.18)	0.16 (0.12 to 0.20)
Growth×size	−0.01 (−0.01 to 0.00)	−0.01 (−0.02 to to 001)	−0.02 (−0.04 to to 0.01)
Rate in a stable CZ with a population 100 000 in the NE region	5.6 (4.6 to 6.9)	8.5 (7.0 to 10.3)	16.9 (13.7 to 21.1)
				Observations	705	709	699
R2	0.1641	0.1644	0.1756

95% CIs in parentheses.

Models were adjusted for a categorical variable representing region (five categories): NE (reference), MW, S, W and region combinations. Twenty-seven CZs crossed regional lines creating the combinations NE–MW (n=1), NE–S (n=1), MW–S (n=17), MW–W (n=5) and S–W (n=3).

Population size was centred at 100 000 (approximately the median population).

Coefficients for growth were exponentiated and calculated as a per cent change in opioid overdose mortality rate per 1% growth in population using the formula (2β–1)×100. Rates in stable CZ with population 100 000 for each period are shown for reference and are derived from the exponentiated coefficients.

CZ, commuting zone; MW, Midwest; NE, Northeast; S, South; W, West.

In CZs at the median population (about 100 000 residents), greater population growth was associated with a decline in opioid mortality: 1.0% increase in a CZ’s population growth, that is, change in population over time, was associated with a decline in the opioid overdose mortality rate of −1.4% (95% CI −2.8% to 0.1%), −4.5% (95% CI −5.8% to −3.2%), and −1.2% (95% CI −4.2% to 1.8%) in 2005–2009, 2010–2014, and 2015–2017, respectively. The negative interaction between population size and growth demonstrates that the association between higher growth and lower opioid mortality was stronger for larger CZs.

Table 3 shows the associations of CZ population size, growth, and their interactions with age-standardised mortality involving prescription opioid, heroin, and synthetic opioid. Larger population size was associated with higher overdose mortality involving prescription opioid, heroin and synthetic opioid in some periods, but the pattern of association was not consistent across opioid types or periods. Population growth was consistently associated with lower mortality by any opioid type in 2010–2014 but this pattern changed in 2015–2017 when coefficients for all drugs moved closer to the null. The interaction between population size and growth did not have a consistent pattern.

Table 3

Association between opioid overdose mortality rates per 100 000 (log 2) by type and population growth and size

	Prescription opioids			Heroin			Synthetic opioids
	2005–2009	2010–2014	2015–2017	2005–2009	2010–2014	2015–2017	2005–2009	2010–2014	2015–2017
Growth (% change for a 1% increase in population for a CZ with population 100 000)	−0.5 (−2.2 to 1.4)	−3.8 (−5.2 to –2.3)	0.1 (−3.0 to 3.2)	−5.9 (−8.6 to –3.2)	−3.6 (−6.4 to –0.7)	−0.9 (−5.1 to 3.5)	−1.3 (−3.1 to 0.4)	−4.1 (−5.6 to –2.7)	−1.9 (−5.8 to 2.2)
Size (% change for a 1% change in average population for a stable CZ)	0.04 (−0.01 to 0.08)	0.07 (0.03 to 0.11)	0.04 (−0.01 to 0.08)	0.024 (−0.06 to 0.11)	0.14 (0.07 to 0.20)	0.18 (0.12 to 0.25)	−0.12 (−0.17 to –0.07)	−0.07 (−0.12 to –0.03)	0.16 (0.10 to 0.22)
Growth×size	−0.005 (−0.01 to 0.00)	−0.013 (−0.02 to –0.00)	−0.015 (−0.03 to 0.00)	0.024 (0.01 to 0.04)	0.011 (−0.01 to 0.03)	−0.010 (−0.04 to 0.02)	−0.002 (−0.01 to 0.01)	−0.002 (−0.01 to 0.01)	−0.036 (−0.06 to –0.01)
Rate in a stable CZ with a population 100 000 in the NE region	2.3 (1.8 to 2.8)	3.8 (3.1 to 4.8)	4.9 (4.0 to 6.1)	0.7 (0.5 to 0.9)	2.1 (1.6 to 2.8)	5.4 (4.1 to 7.2)	1.2 (1.0 to 1.5)	2.2 (1.8 to 2.7)	9.8 (7.5 to 12.9)
										Observations	660	689	673	447	568	588	610	641	623
R2	0.1527	0.1537	0.0712	0.1388	0.1842	0.1734	0.1040	0.1299	0.2648

95% CIs in parentheses.

Models were adjusted for a categorical variable representing region (five categories): NE (reference), MW, S, W and region combinations. Twenty-seven CZs crossed regional lines creating the combinations NE–MW (n=1), NE–S (n=1), MW–S (n=17), MW–W (n=5) and S–W (n=3).

Population size was centred at 100 000 (approximately the median population).

Coefficients for growth were exponentiated and calculated as a per cent change in opioid overdose mortality rate per 1% growth in population using the formula (2β–1)×100. Rates in stable CZ with population 100 000 for each period are shown for reference and are derived from the exponentiated coefficients.

CZ, commuting zone; MW, Midwest; NE, Northeast; S, South; W, West.

Scatterplots showing the relationship between age-standardised opioid mortality rates and population growth can be found in the interactive app, which also shows the CZs with the highest mortality rates.

Discussion

Opioid overdose deaths increased in US CZs between 2005–2009 and 2015–2017. Opioid deaths were disproportionally concentrated among large CZs, consistent with a superlinear scaling behaviour, and this pattern intensified over time. We also found that population growth was associated with lower opioid mortality and that this association was stronger for larger CZs. Last, we observed differences in the patterns by type of opioid, such that deaths involving heroin scaled superlinearly throughout the study period, deaths involving synthetic opioids shifted from a sublinear to a superlinear scaling behaviour, and deaths involving prescription opioids had a scaling behaviour close to linear.

The diversity of analytical approaches used here, examining opioid mortality rates and counts of opioid overdose deaths and using groups of CZs by deciles of population, points to a pattern in which opioid overdose deaths are disproportionally more common in larger CZs after accounting for population size, particularly in more recent years. These patterns are broadly consistent with the hypothesis that opioid overdose deaths increase disproportionally as a function of population size. Under the urban scaling framework, a disproportionally large number of contacts, experienced in more urbanised areas, increases the potential for successful matches, which explains superlinear behaviours observed for outcomes such as wealth and innovation.40 Similar mechanisms might explain the results observed here. For opioids, these successful matches may involve susceptible populations and sources of prescription and illicit opioids. However, we also found great variability in opioid mortality rates across CZs with similar population sizes, particularly in the group of large CZs (deciles 9 and 10) in the last period. This points to the existence of other potential factors in these large CZs that are likely to be related to overdose deaths that are beyond population metrics, and is part of a divergence in mortality rates across geographies in the USA.41 Future research should aim to identify potential explanations related to differences in exposure to opioids or differences in policies to reduce harm among people who use opioids.

Over and above the relationship between larger CZ population and higher overdose rates, we found that CZs with stagnant or declining population have higher opioid overdose mortality rates. Economically distressed areas of the country may experience population decline over time. In fact, between 2008 and 2017, US metropolitan areas concentrated nearly 99% of all jobs and population growth.42 To the extent that population growth is concentrated in cities with strong economies, our findings are also consistent with literature linking high opioid mortality rates and lack of economic opportunity.29–31 43 In addition, the interaction between population size, across CZs, and population growth shows that positive population growth over time becomes an even stronger predictor of lower opioid mortality rates in larger CZs.

Our findings related to mortality by opioid type show that the superlinear scaling of opioid deaths was largely driven by deaths involving heroin and synthetic opioids, with synthetic opioids strongly influencing the results observed in 2015–2017. This finding is consistent with the rise in the supply of fentanyl, an illicitly manufactured synthetic opioid that entered the US drug markets starting around 2013.44 Fentanyl powder and pill presses are used to produce pills that resemble prescription opioid pills or added to heroin by drug dealers to increase its effects.44 Recent evidence suggests that fentanyl is sought out by people who use drugs in some drug markets.45 Data on fentanyl seizures indicate that these synthetic opioids were initially concentrated in specific parts of the country such as the Appalachian region, the Midatlantic and New England.9 This is consistent with our residual maps (interactive app). Our results highlight the potential for future increases in opioid mortality if fentanyl ‘spreads’ to other regions of the country given the potential for superlinear scaling and the large number of people exposed in highly urbanised areas. In light of these findings, broad public health strategies that increase access to naloxone and medication for opioid use disorders46 47 and allow for safer use of opioids and other substances, for example, supervised injection facilities,48 49 are critical to mitigate harm among people who use opioids.

An important feature of our study is the use of CZs as a geographical unit of analysis.50 Studies examining the geographical patterning of drug fatalities have used county-level or state-level outcomes.30 38 51 State-level and county-level analyses are important to assess the effect of policies adopted by states and local governments. However, these analyses do not account for the complex networks across counties that share interconnected economies. CZs reflect commuting patterns and social interactions that may be at the core of the mechanisms associated with the scaling nature of opioid overdose deaths.

This study has limitations. First, this study does not aim to examine potential causal mechanisms between population metrics and opioid overdose deaths. Several factors, including broader social determinants of health and factors of the local context, have a role in the occurrence of overdose deaths.52 In addition, strategies such as harm reduction programmes have the potential to substantially reduce the burden of opioid overdose mortality.51 If the distribution of these factors is also associated with population size or growth, the associations we found can be a consequence of confounding. We recognise the complexity of the determinants of substance use disorders and the exploratory nature of our work. Nonetheless, we believe that our approach can be a valuable way to understand broad patterns of opioid mortality. In particular, using population growth as a potential proxy for other macro-level determinants can be helpful to predict which areas may be at a greater risk of increases in opioid overdose mortality.

Lastly, there are important limitations related to the ascertainment of the cause of death and the use of contributory cause of death to classify drug poisoning deaths, particularly when examining overdose deaths by type of drug.53 Determining which drugs were involved in an overdose often requires comprehensive toxicology testing, including testing of novel psychoactive substances, but in reality toxicology testing is not always comprehensive or timely.53 We conducted a number of secondary or sensitivity analyses, which showed consistent results. However, we cannot rule out the existence of bias in the ascertainment of opioid overdose deaths, particularly in the analysis by type of opioid (ie, prescription opioids, heroin and synthetic opioids). We assume that the distribution of unspecified drug codes occurs independently of other factors beyond those we adjusted for in the imputation procedure, including the type of death investigation system, local economic and demographic profile, and geography. While the majority of poisoning deaths during our study period were associated with opioids, deaths involving psychostimulants with abuse potential (eg, methamphetamines) and cocaine began to increase around 2013.10 Thus, our results about the last period might have overestimated the number of opioid-associated deaths to some extent, although fentanyl was the drug that showed the greatest increase during that period (either in isolation or associated with other drugs).10 Variation in toxicology testing across jurisdictions may also have affected our results. In this scenario, it is possible that the associations observed were at least partially explained by the existence of more comprehensive toxicology testing in larger metropolitan areas of the country, thus revealing a more accurate picture of the types of drugs associated with overdoses in these jurisdictions.

In summary, our study shows that population dynamics of local areas, specifically size and growth, may have a role in the distribution of opioid deaths in the USA. Our findings expand the understanding of the opioid epidemic by demonstrating that several large and mid-sized CZs across the country have disproportionally high number of opioid overdose deaths. The potential superlinear scaling of opioid overdose deaths due to recent economic and social consequences of the COVID-19 pandemic is a particularly worrisome implication from these results.