Study setting and design

We conducted a population-based repeated cross-sectional study of all accidental toxicity deaths where alcohol, stimulants, benzodiazepines or opioids directly contributed to death in Ontario, Canada from January 1^st 2018 to June 30^th 2022. We defined a substance-related death as an accidental acute toxicity death which resulted from the direct contribution of the consumed substance, regardless of how that substance was obtained (i.e., both pharmaceutical and non-pharmaceutical use). We excluded deaths that occurred among individuals without a valid Ontario health card for whom data linkage was not possible. The study is reported as per the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) guidelines [19] with a preliminary analysis of this data previously made available online [20,21].

Data sources

We obtained all data for this analysis from ICES, formerly known as the Institute for Clinical Evaluative Science, an independent, non-profit research institute whose legal status under Ontario health information privacy law allows it to collect and analyse healthcare and demographic data without consent, for health system evaluation and improvement. We first accessed this data on July 13, 2023. Coroners in Ontario investigate all sudden and unexpected deaths. For this study, we used two datasets from the Office of the Chief Coroner of Ontario/Ontario Forensic Pathology Service. We used the Drug and Drug/Alcohol Related Death database to identify confirmed opioid toxicity deaths (with or without other substance involvement) and a newly derived database that includes data on alcohol, stimulant and benzodiazepine related-deaths where opioids did not directly contribute to death to support complete capture of all substance-related toxicity deaths in Ontario involving at least one of these four substances. Sociodemographic characteristics and population denominators were captured using the Registered Persons Database. We used the Canadian Institute for Health Information’s National Ambulatory Care Reporting System and Discharge Abstract Database, the Ontario Mental Health Reporting System and the Ontario Health Insurance Plan to identify emergency department (ED) visits, acute hospital admissions, mental health related hospital admissions and outpatient visits for substance use disorders. All the above datasets were linked using anonymized unique encoded identifiers and analysed at ICES. The use of the data in this study is authorized under Section 45 of Ontario’s Personal Health Information Protection Act (PHIPA) and does not require review by a Research Ethics Board.

Trends and characteristics of individuals who died from a substance-related toxicity

Over the study period, we reported monthly population-adjusted death rates (per 100,000 population) directly attributable to each of alcohol, stimulants, benzodiazepines and opioids (these groups were not mutually exclusive as substance toxicity death may be attributed to more than one substance group). We stratified monthly rates of overall substance-related toxicity deaths by the number of substances directly contributing to death (1, 2, or 3+). We identified all combinations of substances directly contributing to death over the study period and determined the most common substance combinations.

We summarized demographic characteristics of all individuals who died of a substance-related toxicity over the study period, including age (grouped as: < 25, 25–44, 45–64 and ≥65 years old),[22] sex, neighborhood income quintile and location of residence (urban/rural and Northern/Southern Ontario). We also evaluated circumstances surrounding deaths, including the location of incident, and the type and number of substances directly contributing to death. Finally, we summarized prior healthcare encounters before death (S1 Table), including hospital encounters for non-fatal substance-related toxicities in the past year and substance use disorder diagnoses identified during hospitalizations (both any diagnosis type and most responsible diagnosis type) or outpatient physician visits (S2 Table).

Statistical analyses

To examine the impact of the COVID-19 pandemic-related state of emergency on the trends of substance-related toxicity deaths over time, we used a combination of Joinpoint regression and ARIMA models [23]. Joinpoint regression models are recommended in situations where there is no pre-identified notion of where breakpoints in trend data may occur and can be used in the analysis of any outcome assessed over equally distributed time points [24]. Prior research has shown that the impact of the COVID-19 pandemic on rising rates of substance-related deaths began to plateau in late 2020 and early 2021 [25,26], however it was not possible for us to determine a priori the time at which this change in trend would have started. The Joinpoint software identifies points where a trend significantly changes, using the simplest model the data will allow [27]. We used the recommended Grid Search method [28] to select the best fitting Joinpoint model for our data. Using this model, we identified one monthly time-point (August 2020) which represented a significant breakpoint following the onset of the pandemic.

We considered this time-point as an intervention and conducted an ARIMA analysis to simultaneously model the onset of the COVID-19 pandemic-related state of emergency in March 2020 and this additional intervention point. We first plotted our monthly data to visually examine trends with our two intervention points to check for gradual or immediate changes. We determined that the onset of COVID-19 and the additional intervention point in September 2020 should both be fit using a ramp transfer function due to the gradual impact on death rates. As such, we did not consider pulse and step functions. Our ramp transfer functions began in March and September 2020, respectively and increased by a value of one each month until the end of the study period [29]. Our parameter estimates therefore represent the average change in the rate of substance-related toxicities each month. We determined that seasonal differencing was required to achieve stationarity, confirmed with the Dickey-Fuller test to ensure the unit root was removed from our time series. Although seasonal ARIMA (SARIMA) models have been shown to improve forecasting accuracy in the presence of strong seasonality [30–32], this was not the focus of our analysis. To minimize model complexity and reduce the risk of overfitting [33,34], we selected an ARIMA model with seasonal differencing. We used the residual ACF, partial autocorrelation function, and inverse autocorrelation function plots to select appropriate model parameters. Model diagnostics including Akaike’s Information Criteria (AIC), Bayesian Information Criteria (BIC) and the auto-correlation plots were used to select our final model parameters (p, q, d) (S3 Table). We used the Ljung-Box chi-square test for white noise to ensure that there was no significant autocorrelation in model residuals.

Finally, to test for changes in the prevalence of combinations of substances directly contributing to toxicity deaths between the first and last 12 months of the study period, we conducted chi-square tests using a type 1 error rate of 0.05 to define statistical significance. All analyses were conducted using SAS version 9.4 (SAS institute, Cary, North Carolina, USA) and Joinpoint Trend Analysis Software from the National Cancer Institute website (Version 4.8.0.1).

Involvement of people with lived experience

The Ontario Drug Policy Research Network’s opioid-related research is informed by a Lived Experience Advisory Group, comprised of individuals with lived and living experience using opioids who were consulted before, during, and after the study design process. Subsequently, we identified several additional people with lived and living experience with a range of substances including opioids, stimulants, benzodiazepines, and alcohol. Throughout the conduct of this study, these individuals participated as study team members, providing feedback on the study scope, methodology, contextualization, and reporting of results. All people with lived and living experience were compensated for their time spent on this study and were offered co-authorship (or when preferred by the individual, acknowledgment) in this manuscript.

Characteristics of individuals

Over the study period, we identified 10,785 substance toxicity deaths of which 10,022 (92.9%) met our inclusion criteria. Among these 10,022 substance-related toxicity deaths, opioids directly contributed to the majority of deaths (84.1%; N = 8,431), followed by stimulants (60.9%; N = 6,108), alcohol (13.4%; N = 1,346) and benzodiazepines (7.8%; N = 782; Table 1). Overall, 74.9% (N = 7,505) were among males, approximately half were aged 25–44 (51.9%; N = 5,197), and 37.6% were aged 45–64 (N = 3,768). The location of incident was recorded for 9,892 individuals, with the majority occurring at a private residence (72.8%; N = 7,298), followed by other indoor spaces [excluding hospitals] (8.2%; N = 820), outdoors (7.2%; N = 720), and rooming houses or collective dwellings (5.9%; N = 596). Overall, deaths were skewed towards neighbourhoods with the lowest income quintile (40.9%; N = 4,101). While the majority of deaths occurred in urban areas (89.3%) and Southern Ontario (86.9%), the rate of substance-toxicity deaths in Northern Ontario was almost three times higher compared to Southern Ontario (47.9 per 100,000 compared to 16.9 per 100,000, in 2021), and rates were similar between urban and rural areas (18.5 per 100,000 and 17.0 per 100,000, respectively).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Descriptive characteristics, overall substance-related toxicity deaths.

https://doi.org/10.1371/journal.pone.0324732.t001

Combinations of substances involved in toxicity deaths

We observed a 72.1% increase in the number of substance toxicity deaths from the first 12 months of the study period (N = 1,570) compared to the last 12 months (N = 2,702), and the combination of substances directly contributing to death changed considerably over this time (Fig 1). Notably, in the first 12 months of the study period, deaths involving opioids alone occurred most frequently (31.3% of all substance toxicity deaths; N = 492), followed by deaths involving both opioids and stimulants (30.7%, N = 482). However, the proportion of deaths involving both opioids and stimulants grew significantly, reaching 45.3% of deaths by the last 12 months of the study period (N = 1,225; p < 0.001). In parallel, the proportion of deaths involving opioids alone fell significantly over this period, reaching 26.2% by the last 12 months of the study (N = 708; p < 0.001). Although much less common, over the study period there were also significant declines in the prevalence of deaths attributed to alcohol only (3.3% to 1.5%; p < 0.001) and where benzodiazepine directly contributed to death in combination with either opioids (5.2% to 1.3%; p < 0.001), or both stimulants and opioids (2.7% to 1.8%; p = 0.04).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Proportion of substance-related toxicity deaths directly attributable to each substance combination, compared during the first 12 months and last 12 months of the study period. * indicates values of p < 0.05. Note: Categories <1% are not reported here due to space constraints.

https://doi.org/10.1371/journal.pone.0324732.g001

Trends in rates of substance-related toxicity deaths

The monthly rate of substance-related toxicity deaths in Ontario rose throughout the study period, increasing 2-fold from 0.7 per 100,000 in January 2018 up to 1.4 per 100,000 in June 2022 (Fig 2). Using our selected ARIMA model (AIC = -14.60, BIC = -7.75), we observed a significant ramp increase in the rate of substance-related toxicities in March 2020 (+0.11 per 100,000 monthly, 95% CI: 0.01, 0.21 p = 0.03), followed by significant decrease in the trend beginning in September 2020 (-0.14 per 100,000 monthly, 95% CI: -0.23, -0.04, p = 0.01; Table 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Rate of substance-related toxicities: results of the interventional ARIMA analyses.

https://doi.org/10.1371/journal.pone.0324732.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Monthly substance-related toxicity death rates (per 100,000 population), overall and stratified by the number of substances directly involved in death.

^*Vertical bars represents the declaration of a pandemic-related state of emergency for COVID-19 and the additional intervention point identified using joinpoint models.

https://doi.org/10.1371/journal.pone.0324732.g002

In total, 4,324 people (43.1%) who died from a substance-related toxicity death had one substance contributing to death, while 4,790 people had two substances contributing (47.8%) and 908 people (9.1%) had three or more substances directly involved in death. When stratified by the number of substances involved, three or more substances remained the least common throughout the study period, but monthly rates still rose over the four year period (from 0.06 up to 0.07 per 100,000 from January 2018 to June 2022) (Fig 2). In contrast, deaths involving one or two substances were more common. While the rate of single-substance deaths increased 25% (from 0.4 to 0.5 per 100,000 from January 2018 to June 2022), the rate of deaths involving two substances increased nearly 2.5 fold (from 0.3 to 0.7 per 100,000 from January 2018 to June 2022) and was the most common category by the end of study period.

Prior healthcare encounters for substance-related toxicity and substance use disorder

Approximately one-fifth of individuals (19.9%; N = 1,994) were treated in a hospital setting for a substance-related toxicity involving opioids, stimulants, benzodiazepines and/or alcohol in the year prior to death (Table 3). The majority of these toxicities were opioid-related (17.4%; N = 1,748), followed by stimulant-related toxicities (5.0%; N = 499), benzodiazepine-related toxicities (2.0%; N = 204) and alcohol-related toxicities (1.0%; N = 103). Additionally, 60.9% of people who died from a substance-related toxicity (N = 6,098) had a substance use disorder diagnosis (any type) at time of death. Specifically, almost one-quarter of people had a hospital visit relating to opioid use disorder (23.2%; N = 2,323), stimulant use disorder (22.1%; N = 2,219) and/or alcohol use disorder (24.2%; N = 2,429) in the past five years; whereas hospital visits related to benzodiazepine use disorder diagnoses were much less common (2.5%; N = 252). Similarly, over one-third of people had an outpatient physician visit related to a substance use disorder in the past year (excluding alcohol use disorder) (36.6%; N = 3,671), with 5.4% (N = 537) having an outpatient visit relating to alcohol use disorder.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Previous healthcare encounters among individuals who died from a substance-related toxicity.

https://doi.org/10.1371/journal.pone.0324732.t003

Limitations

A strength of this study includes the use of a novel database to provide insights into the trends of substance-related toxicity deaths across Ontario. However, a few limitations warrant discussion. First, we excluded individuals who did not have a valid Ontario health card (approximately 7% of all deaths) to allow for linkage to healthcare records. Therefore, these analyses are an underestimate of the full burden of substance toxicity deaths in Ontario over the study period and the generalizability of our findings may be limited to people who are disconnected from the healthcare system. Second, we only included confirmed deaths where alcohol, stimulants, benzodiazepines, or opioids directly contributed to death. It is possible that suspected substance toxicity deaths may be later determined by the investigating coroner to be related to the above substances, and these were not included in this study. Data from the Chief Coroner of Ontario on drug toxicity deaths are received only for death investigations that are nearly complete and as such we anticipate the potential underestimation of death rates in our study is minimal. Third, there is no validated definition of substance use disorders using health administrative data. Therefore, we relied on hospital records over a five year period, and physician visits in the one year prior to death to identify people with a likely active substance use disorder at time of death. However, it is possible that misclassification exists in this definition. Fourth, we only include benzodiazepine toxicity deaths where a benzodiazepine was identified as a direct contributor to death. Because thresholds are not well established for benzodiazepine toxicity in post-mortem toxicology, our findings may not accurately reflect the presence of benzodiazepine detection in substance toxicity deaths, which has reached over 45% in other studies of opioid toxicity deaths in Ontario during the COVID-19 pandemic [54]. Fifth, the use of seasonal differencing in our ARIMA model to achieve stationarity may complicate the interpretation of trend components, as it removes the seasonal pattern from the data. This can make it more difficult to distinguish between long-term trends and short-term fluctuations. Finally, our reliance on healthcare administrative data limited our ability to capture some pandemic-related factors such as stress, job loss, and housing instability, while the joinpoint model, though effective for detecting trend shifts, may miss gradual or complex changes and does not account for confounding variables. These limitations may affect the generalizability of our findings, as the observed trends may not reflect the experiences of people disproportionately impacted by the pandemic. Future research could benefit from integrating complementary data sources such as surveys or interviews to help capture broader social impacts.