Data-driven modelling of brain activity using neural networks, diffusion maps, and the Koopman operator.

We propose a machine-learning approach to construct reduced-order models (ROMs) to predict the long-term out-of-sample dynamics of brain activity (and in general, high-dimensional time series), focusing mainly on task-dependent high-dimensional fMRI time series. Our approach is a three stage one. First, we exploit manifold learning and, in particular, diffusion maps (DMs) to discover a set of variables that parametrize the latent space on which the emergent high-dimensional fMRI time series evolve. Then, we construct ROMs on the embedded manifold via two techniques: Feedforward Neural Networks (FNNs) and the Koopman operator. Finally, for predicting the out-of-sample long-term dynamics of brain activity in the ambient fMRI space, we solve the pre-image problem, i.e., the construction of a map from the low-dimensional manifold to the original high-dimensional (ambient) space by coupling DMs with Geometric Harmonics (GH) when using FNNs and the Koopman modes per se. For our illustrations, we have assessed the performance of the two proposed schemes using two benchmark fMRI time series: (i) a simplistic five-dimensional model of stochastic discrete-time equations used just for a "transparent" illustration of the approach, thus knowing a priori what one expects to get, and (ii) a real fMRI dataset with recordings during a visuomotor task. We show that the proposed Koopman operator approach provides, for any practical purposes, equivalent results to the FNN-GH approach, thus bypassing the need to train a non-linear map and to use GH to extrapolate predictions in the ambient space; one can use instead the low-frequency truncation of the DMs function space of L2-integrable functions to predict the entire list of coordinate functions in the ambient space and to solve the pre-image problem.

Learning black- and gray-box chemotactic PDEs/closures from agent based Monte Carlo simulation data.

We propose a machine learning framework for the data-driven discovery of macroscopic chemotactic Partial Differential Equations (PDEs)-and the closures that lead to them- from high-fidelity, individual-based stochastic simulations of Escherichia coli bacterial motility. The fine scale, chemomechanical, hybrid (continuum-Monte Carlo) simulation model embodies the underlying biophysics, and its parameters are informed from experimental observations of individual cells. Using a parsimonious set of collective observables, we learn effective, coarse-grained "Keller-Segel class" chemotactic PDEs using machine learning regressors: (a) (shallow) feedforward neural networks and (b) Gaussian Processes. The learned laws can be black-box (when no prior knowledge about the PDE law structure is assumed) or gray-box when parts of the equation (e.g. the pure diffusion part) is known and "hardwired" in the regression process. More importantly, we discuss data-driven corrections (both additive and functional), to analytically known, approximate closures.

Parsimonious physics-informed random projection neural networks for initial value problems of ODEs and index-1 DAEs.

We present a numerical method based on random projections with Gaussian kernels and physics-informed neural networks for the numerical solution of initial value problems (IVPs) of nonlinear stiff ordinary differential equations (ODEs) and index-1 differential algebraic equations (DAEs), which may also arise from spatial discretization of partial differential equations (PDEs). The internal weights are fixed to ones while the unknown weights between the hidden and output layer are computed with Newton's iterations using the Moore-Penrose pseudo-inverse for low to medium scale and sparse QR decomposition with L 2 regularization for medium- to large-scale systems. Building on previous works on random projections, we also prove its approximation accuracy. To deal with stiffness and sharp gradients, we propose an adaptive step-size scheme and address a continuation method for providing good initial guesses for Newton iterations. The "optimal" bounds of the uniform distribution from which the values of the shape parameters of the Gaussian kernels are sampled and the number of basis functions are "parsimoniously" chosen based on bias-variance trade-off decomposition. To assess the performance of the scheme in terms of both numerical approximation accuracy and computational cost, we used eight benchmark problems (three index-1 DAEs problems, and five stiff ODEs problems including the Hindmarsh-Rose neuronal model of chaotic dynamics and the Allen-Cahn phase-field PDE). The efficiency of the scheme was compared against two stiff ODEs/DAEs solvers, namely, ode15s and ode23t solvers of the MATLAB ODE suite as well as against deep learning as implemented in the DeepXDE library for scientific machine learning and physics-informed learning for the solution of the Lotka-Volterra ODEs included in the demos of the library. A software/toolbox in Matlab (that we call RanDiffNet) with demos is also provided.

Time-series forecasting using manifold learning, radial basis function interpolation, and geometric harmonics.

We address a three-tier numerical framework based on nonlinear manifold learning for the forecasting of high-dimensional time series, relaxing the "curse of dimensionality" related to the training phase of surrogate/machine learning models. At the first step, we embed the high-dimensional time series into a reduced low-dimensional space using nonlinear manifold learning (local linear embedding and parsimonious diffusion maps). Then, we construct reduced-order surrogate models on the manifold (here, for our illustrations, we used multivariate autoregressive and Gaussian process regression models) to forecast the embedded dynamics. Finally, we solve the pre-image problem, thus lifting the embedded time series back to the original high-dimensional space using radial basis function interpolation and geometric harmonics. The proposed numerical data-driven scheme can also be applied as a reduced-order model procedure for the numerical solution/propagation of the (transient) dynamics of partial differential equations (PDEs). We assess the performance of the proposed scheme via three different families of problems: (a) the forecasting of synthetic time series generated by three simplistic linear and weakly nonlinear stochastic models resembling electroencephalography signals, (b) the prediction/propagation of the solution profiles of a linear parabolic PDE and the Brusselator model (a set of two nonlinear parabolic PDEs), and (c) the forecasting of a real-world data set containing daily time series of ten key foreign exchange rates spanning the time period 3 September 2001-29 October 2020.

Coarse-scale PDEs from fine-scale observations via machine learning.

Complex spatiotemporal dynamics of physicochemical processes are often modeled at a microscopic level (through, e.g., atomistic, agent-based, or lattice models) based on first principles. Some of these processes can also be successfully modeled at the macroscopic level using, e.g., partial differential equations (PDEs) describing the evolution of the right few macroscopic observables (e.g., concentration and momentum fields). Deriving good macroscopic descriptions (the so-called "closure problem") is often a time-consuming process requiring deep understanding/intuition about the system of interest. Recent developments in data science provide alternative ways to effectively extract/learn accurate macroscopic descriptions approximating the underlying microscopic observations. In this paper, we introduce a data-driven framework for the identification of unavailable coarse-scale PDEs from microscopic observations via machine-learning algorithms. Specifically, using Gaussian processes, artificial neural networks, and/or diffusion maps, the proposed framework uncovers the relation between the relevant macroscopic space fields and their time evolution (the right-hand side of the explicitly unavailable macroscopic PDE). Interestingly, several choices equally representative of the data can be discovered. The framework will be illustrated through the data-driven discovery of macroscopic, concentration-level PDEs resulting from a fine-scale, lattice Boltzmann level model of a reaction/transport process. Once the coarse evolution law is identified, it can be simulated to produce long-term macroscopic predictions. Different features (pros as well as cons) of alternative machine-learning algorithms for performing this task (Gaussian processes and artificial neural networks) are presented and discussed.

Amplitude spectrum EEG signal evidence for the dissociation of motor and perceptual spatial working memory in the human brain.

This study investigated the question whether spatial working memory related to movement plans (motor working memory) and spatial working memory related to spatial attention and perceptual processes (perceptual spatial working memory) share the same neurophysiological substrate or there is evidence for separate motor and perceptual working memory streams of processing. Towards this aim, ten healthy human subjects performed delayed responses to visual targets presented at different spatial locations. Two tasks were attained, one in which the spatial location of the target was the goal for a pointing movement and one in which the spatial location of the target was used for a perceptual (yes or no) change detection. Each task involved two conditions: a memory condition in which the target remained visible only for the first 250 ms of the delay period and a delay condition in which the target location remained visible throughout the delay period. The amplitude spectrum analysis of the EEG revealed that the alpha (8-12 Hz) band signal was smaller, while the beta (13-30 Hz) and gamma (30-45 Hz) band signals were larger in the memory compared to the non-memory condition. The alpha band signal difference was confined to the frontal midline area; the beta band signal difference extended over the right hemisphere and midline central area, and the gamma band signal difference was confined to the right occipitoparietal area. Importantly, both in beta and gamma bands, we observed a significant increase in the movement-related compared to the perceptual-related memory-specific amplitude spectrum signal in the central midline area. This result provides clear evidence for the dissociation of motor and perceptual spatial working memory.

Task-oriented machine learning surrogates for tipping points of agent-based models.

We present a machine learning framework bridging manifold learning, neural networks, Gaussian processes, and Equation-Free multiscale approach, for the construction of different types of effective reduced order models from detailed agent-based simulators and the systematic multiscale numerical analysis of their emergent dynamics. The specific tasks of interest here include the detection of tipping points, and the uncertainty quantification of rare events near them. Our illustrative examples are an event-driven, stochastic financial market model describing the mimetic behavior of traders, and a compartmental stochastic epidemic model on an Erdös-Rényi network. We contrast the pros and cons of the different types of surrogate models and the effort involved in learning them. Importantly, the proposed framework reveals that, around the tipping points, the emergent dynamics of both benchmark examples can be effectively described by a one-dimensional stochastic differential equation, thus revealing the intrinsic dimensionality of the normal form of the specific type of the tipping point. This allows a significant reduction in the computational cost of the tasks of interest.

Single-trial magnetoencephalography signals encoded as an unfolding decision process.

The model of a stochastic decision process unfolding in motor and premotor regions of the brain was encoded in single-trial magnetoencephalographic (MEG) recordings while ten healthy subjects performed a sensorimotor Reaction Time (RT) task. The duration of single-trial MEG signals preceding the motor response, recorded over the motor cortex contralateral to the responding hand, co-varied with RT across trials according to the model's prediction. Furthermore, these signals displayed the same properties of a "rising-to-a-fixed-threshold" decision process as posited by the model and observed in the activity of single neurons in the primate cortex. The present findings demonstrate that non-averaged, single-trial MEG recordings can be used to test models of cognitive processes, like decision-making, in humans.

Designing social distancing policies for the COVID-19 pandemic: A probabilistic model predictive control approach.

The effective control of the COVID-19 pandemic is one the most challenging issues of recent years. The design of optimal control policies is challenging due to a variety of social, political, economical and epidemiological factors. Here, based on epidemiological data reported in recent studies for the Italian region of Lombardy, which experienced one of the largest and most devastating outbreaks in Europe during the first wave of the pandemic, we present a probabilistic model predictive control (PMPC) approach for the systematic study of what if scenarios of social distancing in a retrospective analysis for the first wave of the pandemic in Lombardy. The performance of the proposed PMPC was assessed based on simulations of a compartmental model that was developed to quantify the uncertainty in the level of the asymptomatic cases in the population, and the synergistic effect of social distancing during various activities, and public awareness campaign prompting people to adopt cautious behaviors to reduce the risk of disease transmission. The PMPC takes into account the social mixing effect, i.e. the effect of the various activities in the potential transmission of the disease. The proposed approach demonstrates the utility of a PMPC approach in addressing COVID-19 transmission and implementing public relaxation policies.

Construction of embedded fMRI resting-state functional connectivity networks using manifold learning.

We construct embedded functional connectivity networks (FCN) from benchmark resting-state functional magnetic resonance imaging (rsfMRI) data acquired from patients with schizophrenia and healthy controls based on linear and nonlinear manifold learning algorithms, namely, Multidimensional Scaling, Isometric Feature Mapping, Diffusion Maps, Locally Linear Embedding and kernel PCA. Furthermore, based on key global graph-theoretic properties of the embedded FCN, we compare their classification potential using machine learning. We also assess the performance of two metrics that are widely used for the construction of FCN from fMRI, namely the Euclidean distance and the cross correlation metric. We show that diffusion maps with the cross correlation metric outperform the other combinations.

A bulletin from Greece: a health system under the pressure of the second COVID-19 wave.

In this article, we analyze the cascade of events since the beginning of the coronavirus disease 2019 (COVID-19) pandemic in Greece, with emphasis on the crisis' management so as to preserve the functionality of the national health system, which remains vulnerable due to the financial recession of the previous decade and chronic shortcomings . We compare and contrast the situation during the first and second epidemic wave. Understanding what possibly went wrong and when, is crucial . Such knowledge provides valuable guidance for the confrontation of the strong second wave that we are currently facing in Europe and other regions around the globe, as well as for the future waves that may follow.

Lessons from the devastating impact of the first COVID-19 wave in Italy.

Herein, we are critically examining the chain of events and discussing previously unrecognized factors that led to the 'perfect COVID-19 storm' in northern Italy during the first epidemic wave in spring 2020. SARS-CoV-2 was circulating uncontrollably at least for five weeks before the adoption of containment measures, and the role of exponential growth in the spread of the virus, conveyed by a high R0, was likely underestimated. An understanding of this failure's causes and contexts will help us to control the strong second wave of the pandemic we are now facing in Europe, and to be better prepared for future outbreaks.

ISOMAP and machine learning algorithms for the construction of embedded functional connectivity networks of anatomically separated brain regions from resting state fMRI data of patients with Schizophrenia.

We construct Functional Connectivity Networks (FCN) from resting state fMRI (rsfMRI) recordings towards the classification of brain activity between healthy and schizophrenic subjects using a publicly available dataset (the COBRE dataset) of 145 subjects (74 healthy controls and 71 schizophrenic subjects). First, we match the anatomy of the brain of each individual to the Desikan-Killiany brain atlas. Then, we use the conventional approach of correlating the parcellated time series to construct FCN and ISOMAP, a nonlinear manifold learning algorithm to produce low-dimensional embeddings of the correlation matrices. For the classification analysis, we computed five key local graph-theoretic measures of the FCN and used the LASSO and Random Forest (RF) algorithms for feature selection. For the classification we used standard linear Support Vector Machines. The classification performance is tested by a double cross-validation scheme (consisting of an outer and an inner loop of "Leave one out" cross-validation (LOOCV)). The standard cross-correlation methodology produced a classification rate of 73.1%, while ISOMAP resulted in 79.3%, thus providing a simpler model with a smaller number of features as chosen from LASSO and RF, namely the participation coefficient of the right thalamus and the strength of the right lingual gyrus.

Construction of functional brain connectivity networks from fMRI data with driving and modulatory inputs: an extended conditional Granger causality approach.

We propose a numerical-based approach extending the conditional MVAR Granger causality (MVGC) analysis for the construction of directed connectivity networks in the presence of both exogenous/stimuli and modulatory inputs. The performance of the proposed scheme is validated using both synthetic stochastic data considering also the influence of haemodynamics latencies and a benchmark fMRI dataset related to the role of attention in the perception of visual motion. The particular fMRI dataset has been used in many studies to evaluate alternative model hypotheses using the Dynamic Causal Modelling (DCM) approach. Based on the use of the Bayes factor, we show that the obtained GC connectivity network compares well to a reference model that has been selected through DCM analysis among other candidate models. Thus, our findings suggest that the proposed scheme can be successfully used as a stand-alone or complementary to DCM approach to find directed causal connectivity patterns in task-related fMRI studies.

Data-based analysis, modelling and forecasting of the COVID-19 outbreak.

Since the first suspected case of coronavirus disease-2019 (COVID-19) on December 1st, 2019, in Wuhan, Hubei Province, China, a total of 40,235 confirmed cases and 909 deaths have been reported in China up to February 10, 2020, evoking fear locally and internationally. Here, based on the publicly available epidemiological data for Hubei, China from January 11 to February 10, 2020, we provide estimates of the main epidemiological parameters. In particular, we provide an estimation of the case fatality and case recovery ratios, along with their 90% confidence intervals as the outbreak evolves. On the basis of a Susceptible-Infectious-Recovered-Dead (SIDR) model, we provide estimations of the basic reproduction number (R0), and the per day infection mortality and recovery rates. By calibrating the parameters of the SIRD model to the reported data, we also attempt to forecast the evolution of the outbreak at the epicenter three weeks ahead, i.e. until February 29. As the number of infected individuals, especially of those with asymptomatic or mild courses, is suspected to be much higher than the official numbers, which can be considered only as a subset of the actual numbers of infected and recovered cases in the total population, we have repeated the calculations under a second scenario that considers twenty times the number of confirmed infected cases and forty times the number of recovered, leaving the number of deaths unchanged. Based on the reported data, the expected value of R0 as computed considering the period from the 11th of January until the 18th of January, using the official counts of confirmed cases was found to be ∼4.6, while the one computed under the second scenario was found to be ∼3.2. Thus, based on the SIRD simulations, the estimated average value of R0 was found to be ∼2.6 based on confirmed cases and ∼2 based on the second scenario. Our forecasting flashes a note of caution for the presently unfolding outbreak in China. Based on the official counts for confirmed cases, the simulations suggest that the cumulative number of infected could reach 180,000 (with a lower bound of 45,000) by February 29. Regarding the number of deaths, simulations forecast that on the basis of the up to the 10th of February reported data, the death toll might exceed 2,700 (as a lower bound) by February 29. Our analysis further reveals a significant decline of the case fatality ratio from January 26 to which various factors may have contributed, such as the severe control measures taken in Hubei, China (e.g. quarantine and hospitalization of infected individuals), but mainly because of the fact that the actual cumulative numbers of infected and recovered cases in the population most likely are much higher than the reported ones. Thus, in a scenario where we have taken twenty times the confirmed number of infected and forty times the confirmed number of recovered cases, the case fatality ratio is around ∼0.15% in the total population. Importantly, based on this scenario, simulations suggest a slow down of the outbreak in Hubei at the end of February.

Electroencephalography source localization analysis in epileptic children during a visual working-memory task.

We localize the sources of brain activity of children with epilepsy based on electroencephalograph (EEG) recordings acquired during a visual discrimination working memory task. For the numerical solution of the inverse problem, with the aid of age-specific MRI scans processed from a publicly available database, we use and compare three regularization numerical methods, namely the standardized low resolution brain electromagnetic tomography (sLORETA), the weighted minimum norm estimation (wMNE) and the dynamic statistical parametric mapping (dSPM). We show that all three methods provide the same spatio-temporal patterns of differences between the groups of epileptic and control children. In particular, our analysis reveals statistically significant differences between the two groups in regions of the parietal cortex indicating that these may serve as "biomarkers" for diagnostic purposes and ultimately localized treatment.

Tracing day-zero and forecasting the COVID-19 outbreak in Lombardy, Italy: A compartmental modelling and numerical optimization approach.

INTRODUCTION: Italy became the second epicenter of the novel coronavirus disease 2019 (COVID-19) pandemic after China, surpassing by far China's death toll. The disease swept through Lombardy, which remained in lockdown for about two months, starting from the 8th of March. As of that day, the isolation measures taken in Lombardy were extended to the entire country. Here, assuming that effectively there was one case "zero" that introduced the virus to the region, we provide estimates for: (a) the day-zero of the outbreak in Lombardy, Italy; (b) the actual number of asymptomatic infected cases in the total population until March 8; (c) the basic (R0)and the effective reproduction number (Re) based on the estimation of the actual number of infected cases. To demonstrate the efficiency of the model and approach, we also provide a tentative forecast two months ahead of time, i.e. until May 4, the date on which relaxation of the measures commenced, on the basis of the COVID-19 Community Mobility Reports released by Google on March 29. METHODS: To deal with the uncertainty in the number of the actual asymptomatic infected cases in the total population Volpert et al. (2020), we address a modified compartmental Susceptible/ Exposed/ Infectious Asymptomatic/ Infected Symptomatic/ Recovered/ Dead (SEIIRD) model with two compartments of infectious persons: one modelling the cases in the population that are asymptomatic or experience very mild symptoms and another modelling the infected cases with mild to severe symptoms. The parameters of the model corresponding to the recovery period, the time from the onset of symptoms to death and the time from exposure to the time that an individual starts to be infectious, have been set as reported from clinical studies on COVID-19. For the estimation of the day-zero of the outbreak in Lombardy, as well as of the "effective" per-day transmission rate for which no clinical data are available, we have used the proposed SEIIRD simulator to fit the numbers of new daily cases from February 21 to the 8th of March. This was accomplished by solving a mixed-integer optimization problem. Based on the computed parameters, we also provide an estimation of the basic reproduction number R0 and the evolution of the effective reproduction number Re. To examine the efficiency of the model and approach, we ran the simulator to "forecast" the epidemic two months ahead of time, i.e. from March 8 to May 4. For this purpose, we considered the reduction in mobility in Lombardy as released on March 29 by Google COVID-19 Community Mobility Reports, and the effects of social distancing and of the very strict measures taken by the government on March 20 and March 21, 2020. RESULTS: Based on the proposed methodological procedure, we estimated that the expected day-zero was January 14 (min-max rage: January 5 to January 23, interquartile range: January 11 to January 18). The actual cumulative number of asymptomatic infected cases in the total population in Lombardy on March 8 was of the order of 15 times the confirmed cumulative number of infected cases, while the expected value of the basic reproduction number R0 was found to be 4.53 (min-max range: 4.40- 4.65). On May 4, the date on which relaxation of the measures commenced the effective reproduction number was found to be 0.987 (interquartiles: 0.857, 1.133). The model approximated adequately two months ahead of time the evolution of reported cases of infected until May 4, the day on which the phase I of the relaxation of measures was implemented over all of Italy. Furthermore the model predicted that until May 4, around 20% of the population in Lombardy has recovered (interquartile range: ∼10% to ∼30%).

Bautin bifurcations in a forest-grassland ecosystem with human-environment interactions.

Ecosystems may be characterized by a complex dynamical behaviour where external disturbances and/or internal perturbations may trigger sudden/irreversible changes, called catastrophic shifts. Simple mathematical models in the form of ordinary and/or partial differential equations have been proposed to approximate in a qualitatively manner the observed complex phenomena, where catastrophic shifts are determined by bifurcation points. In this work, we show that in ecosystems, gradual/smooth changes may be transformed in sudden/catastrophic shifts as a consequence of codimension-2 bifurcations. We stress the importance of using the full arsenal of numerical bifurcation theory to systematically identify and characterize criticalities in ecological models in the 2D parameter space. For our demonstrations, we revisit the analysis of a simple model of a forest-grassland mosaic ecosystem constructing the 2D bifurcation diagram with respect to the impact of human influence and that of natural causes. Our numerical analysis reveals that this simple model is able to approximate both abrupt (catastrophic) and smooth transitions as the system undergoes Bautin bifurcations.

Varicella zoster virus transmission dynamics in Vojvodina, Serbia.

This study aimed at establishing baseline key epidemiological parameters for varicella zoster virus (VZV) infection in Vojvodina, Serbia, with the ultimate goal to quantify the VZV transmission potential in the population. Seroprevalence data generated during the first large cross-sectional VZV serosurvey were modelled, using a two-tiered modelling approach to calculate age-specific forces of infection (FOI), the basic reproduction number (R0) and herd immunity threshold (H). Seroprevalence and modelling data were compared with corresponding pre-vaccination epidemiological parameters from 11 countries participating in the European Sero-Epidemiology Network 2 (ESEN2) project. Serbia fits into the general dynamic VZV transmission patterns in Europe in the pre-vaccine era, with estimated R0 = 4.12, (95% CI: 2.69-7.07) and H = 0.76 (95% CI: 0.63-0.86). The highest VZV transmission occurs among preschool children, as evidenced by the estimation of the highest FOI (0.22, 95% CI: 0.11-0.34) in the 0.5-4 age group, with a peak FOI of 0.25 at 2.23 years. Seroprevalence was consistently lower in 5-14 year-olds, resulting in considerable shares of VZV-susceptible adolescents (7.3%), and young adults (6%), resembling the situation in a minority of European countries. The obtained key epidemiological parameters showed most intense VZV transmission in preschool children aged <4 years, justifying the consideration of universal childhood immunization in the future. National immunization strategy should consider programs for VZV serologic screening and immunization of susceptible groups, including adolescents and women of reproductive age. This work is an important milestone towards the evaluation of varicella immunization policy options in Serbia.

Correlation of Opioid Mortality with Prescriptions and Social Determinants: A Cross-sectional Study of Medicare Enrollees.

BACKGROUND: The opioid epidemic is an escalating health crisis. We evaluated the impact of opioid prescription rates and socioeconomic determinants on opioid mortality rates, and identified potential differences in prescription patterns by categories of practitioners. METHODS: We combined the 2013 and 2014 Medicare Part D data and quantified the opioid prescription rate in a county level cross-sectional study with data from 2710 counties, 468,614 unique prescribers and 46,665,037 beneficiaries. We used the CDC WONDER database to obtain opioid-related mortality data. Socioeconomic characteristics for each county were acquired from the US Census Bureau. RESULTS: The average national opioid prescription rate was 3.86 claims per beneficiary that received a prescription for opioids (95% CI 3.86-3.86). At a county level, overall opioid prescription rates (p < 0.001, Coeff = 0.27) and especially those provided by emergency medicine (p < 0.001, Coeff = 0.21), family medicine physicians (p = 0.11, Coeff = 0.008), internal medicine (p = 0.018, Coeff = 0.1) and physician assistants (p = 0.021, Coeff = 0.08) were associated with opioid-related mortality. Demographic factors, such as proportion of white (p white < 0.001, Coeff = 0.22), black (p black < 0.001, Coeff = - 0.19) and male population (p male < 0.001, Coeff = 0.13) were associated with opioid prescription rates, while poverty (p < 0.001, Coeff = 0.41) and proportion of white population (p white < 0.001, Coeff = 0.27) were risk factors for opioid-related mortality (p model < 0.001, R 2 = 0.35). Notably, the impact of prescribers in the upper quartile was associated with opioid mortality (p < 0.001, Coeff = 0.14) and was twice that of the remaining 75% of prescribers together (p < 0.001, Coeff = 0.07) (p model = 0.03, R 2 = 0.03). CONCLUSIONS: The prescription opioid rate, and especially that by certain categories of prescribers, correlated with opioid-related mortality. Interventions should prioritize providers that have a disproportionate impact and those that care for populations with socioeconomic factors that place them at higher risk.