Research ArticleResearch Article: New Research, Sensory and Motor Systems

Modeling Effects of Variable preBötzinger Complex Network Topology and Cellular Properties on Opioid-Induced Respiratory Depression and Recovery

Grant M. Chou, Nicholas E. Bush, Ryan S. Phillips, Nathan A. Baertsch and Kameron Decker Harris
eNeuro 22 January 2024, 11 (3) ENEURO.0284-23.2023; https://doi.org/10.1523/ENEURO.0284-23.2023

Abstract

The preBötzinger complex (preBötC), located in the medulla, is the essential rhythm-generating neural network for breathing. The actions of opioids on this network impair its ability to generate robust, rhythmic output, contributing to life-threatening opioid-induced respiratory depression (OIRD). The occurrence of OIRD varies across individuals and internal and external states, increasing the risk of opioid use, yet the mechanisms of this variability are largely unknown. In this study, we utilize a computational model of the preBötC to perform several in silico experiments exploring how differences in network topology and the intrinsic properties of preBötC neurons influence the sensitivity of the network rhythm to opioids. We find that rhythms produced by preBötC networks in silico exhibit variable responses to simulated opioids, similar to the preBötC network in vitro. This variability is primarily due to random differences in network topology and can be manipulated by imposed changes in network connectivity and intrinsic neuronal properties. Our results identify features of the preBötC network that may regulate its susceptibility to opioids.

Significance Statement

The neural network in the brainstem that generates the breathing rhythm is disrupted by opioid drugs. However, this response can be surprisingly unpredictable. By constructing computational models of this rhythm-generating network, we illustrate how random differences in the distribution of biophysical properties and connectivity patterns within individual networks can predict their response to opioids, and we show how modulation of these network features can make breathing more susceptible or resistant to the effects of opioids.

Introduction

Opioid-induced respiratory depression (OIRD) is the primary cause of death associated with opioid overdose. Because both the pain-killing and respiratory depressive effects of opioids require the μ-opioid receptor (MOR) encoded by the Oprm1 gene (Sora et al., 1997; Dahan et al., 2001; Baldo and Rose, 2022; Lynch et al., 2023), there are few effective strategies to protect against OIRD without eliminating the beneficial analgesic effects of opioids. Increasing doses of opioid are often required to maintain analgesia as the neural circuits that mediate pain become desensitized to opioids (Freye and Latasch, 2003; Uniyal et al., 2020), putting patients at greater risk of OIRD. However, a dangerous and less well-understood feature of OIRD is its apparent unpredictability (Dahan et al., 2013). Changes in breathing in response to opioid use can vary substantially between individuals and can be surprisingly inconsistent even within the same individual under different internal and external states or contexts (Cherny et al., 2001; Dahan et al., 2013; Fleming et al., 2015).

Although Oprm1 is expressed in many brain regions (Erbs et al., 2015), including those involved in the regulation of breathing (Varga et al., 2020; Baldo and Rose, 2022), one site of particular importance is the PreBötzinger complex (preBötC), a medullary region that is critical for generating the respiratory rhythm (Smith et al., 1991; Gray et al., 1999; Bachmutsky et al., 2020). This network is composed of interacting excitatory and inhibitory interneurons (Winter et al., 2009; Harris et al., 2017; Baertsch et al., 2018). Although inhibitory neurons are important for regulating the frequency and regularity of breathing (Sherman et al., 2015; Baertsch et al., 2018), GABAergic or glycinergic mechanisms do not seem to play a significant role in OIRD in the preBötC (Gray et al., 1999; Bachmutsky et al., 2020). Instead, glutamatergic neurons are the critical component of the preBötC network needed for both rhythmogenesis and OIRD (Greer et al., 1991; Funk et al., 1993; Sun et al., 2019; Bachmutsky et al., 2020). Among the estimated 40–60% of preBötC neurons that express Oprm1 (Gray et al., 1999; Baertsch et al., 2021; Rousseau et al., 2023), activation of MORs has two primary consequences: excitability is suppressed due to activation of a hyperpolarizing current and the strength of excitatory synaptic interactions is reduced (Baertsch et al., 2021). Together, these mechanisms of opioid action act synergistically to undermine the cellular and network mechanisms that mediate preBötC rhythmogenesis.

Neurons in the preBötC have heterogeneous cellular properties, which are readily observed following pharmacological blockade of synaptic interactions. Under these conditions, the intrinsic activity of preBötC neurons is either silent, bursting, or tonic, which largely depends on their persistent sodium conductance (gNaP) and potassium dominated leak conductance (gleak) (Butera et al., 1999; Del Negro et al., 2002; Koizumi and Smith, 2008; Yamanishi et al., 2018; Phillips and Rubin, 2019). However, gNaP, gleak, and the intrinsic activity of preBötC neurons are not fixed but can be dynamically modulated by conditional factors such as neuromodulation and changes in excitability (Del Negro et al., 2001; Doi and Ramirez, 2008; Ramirez et al., 2011). Thus, unlike the discrete activity states of its constituent neurons, when synaptically coupled, the network collectively produces an inspiratory rhythm that can operate along a continuum of states as the ratios of silent, bursting, and tonic neurons change (Butera and Smith, 1999; Burgraff et al., 2021). As previously demonstrated in rhythmic brainstem slices, the preBötC has an optimal configuration of cellular and network properties that results in a maximally stable inspiratory rhythm. These properties are dynamic, and the state of each individual preBötC network relative to its optimal configuration can predict how susceptible rhythmogenesis is to opioids (Burgraff et al., 2021).

Here, we expand on these findings by utilizing computational modeling to perform preBötC network manipulations and analyses that are experimentally intractable to better understand properties of the network that may contribute to the variation in OIRD and to provide proof of concept for perturbations that may render preBötC rhythmogenesis less vulnerable to opioids. We demonstrate that model preBötC networks exhibit variable responses to simulated opioids. This variation in opioid response is best predicted by differences in “fixed” properties of randomly generated networks, specifically the connectivity between different groups of excitatory and inhibitory neurons as well as which neurons in the network express MOR. In contrast, opioid-induced changes in the intrinsic spiking patterns preBötC neurons (silent, bursting, and tonic) do not predict this variation. In networks with high opioid sensitivity, we find that modulation of either gNaP or gleak can render rhythmogenesis more resistant to opioids. These insights help establish a conceptual framework for understanding how the fixed and dynamic properties of the preBötC shape how this vital network responds when challenged with opioids.

Methods

Computational modeling of OIRD in the PreBötC

We model the preBötC network as a random, directed graph of N = 300 nodes, with each node representing a neuron. The dynamical neuron equations are modified from Butera and Smith (1999), Butera et al. (1999), Harris et al. (2017), and Baertsch et al. (2021). First, we have the total membrane current balance equation−CmdVdt=INa+INaP+Ileak+IK+Isyn, exc+Isyn, inh+Isyn, op+Ihyp, op, (1)where the currents are given byINa=gNam∞3(1−n)(V−ENa),INaP=gNaPDNaP(t)mNaP,∞h(V−ENa),Ileak=gleakDleak(t)(V−Eleak),IK=gKn4(V−EK),Isyn, exc=gexc(V−Eexc),Isyn, inh=ginh(V−Einh),Isyn, op=gsyn, op(t)(V−Eexc),Ihyp, op=DopIhyp, op(t). We implemented the terms DNaP(t) and Dleak(t) to simulate time-dependent “drugs” strengthening or weakening NaP and leak channel conductances by varying between 0 and 1. We also added the opioid-modulated synaptic and hyperpolarizing currents Isyn,op and Ihyp,op as a mechanism for biophysical perturbations through changing Ihyp,op(t) and gsyn,op(t). While many terms in these differential equations are time-dependent, we only explicitly highlight the time-dependence of DNaP, Dleak, gsyn,op, and Ihyp,op because these are set exogenously. The gating equations aredndt=n∞−nτn,dhdt=h∞−hτh with voltage-dependent steady states:m∞=(1+exp(V−Emσm))−1,mNaP,∞=(1+exp(V−Em,NaP,∞σm,NaP,∞))−1,n∞=(1+exp(V−Enσn))−1,h∞=(1+exp(V−Ehσh))−1, and time constants:τn=τnbcosh((V−En)/2σn),τh=τhbcosh((V−Eh)/2σh). Finally, synapses are modeled separately for excitatory, inhibitory, and opioid-sensitive presynaptic neurons. Each synapse conductance s is modeled with first-order dynamics:dsdt=(1−s)msyn,∞−sτsyn,msyn,∞=(1+exp(V\,pre−Esσsyn))−1, and gsyn is given as the sum of gsyn,max · s over all incoming synapses to the neuron. This model was implemented in brian2 (Stimberg et al., 2019). The parameters which are shared across all neurons are given in Table 1; other parameters will be described in the rest of the methods. Our code will be included with the final version of this paper.

View this table:
Table 1.

Table of model parameters shared across neurons

Network construction

Our 300 neuron network consists of 60 inhibitory neurons and 240 excitatory neurons. Synapses were randomly constructed, with each neuron having a connection probability of (davg/2)/(N − 1), where davg is the neurons’ average degree (in-degree + out-degree). Our default davg is 6, giving us a connection density of approximately 1%. However, in Figure 3, we increase the connection probability by scaling davg, e.g., davg = 12 results in a 2% connection density.

The intrinsic activity of each neuron is either tonic spiking (T), bursting (B), or quiescent (Q), which is controlled by gleak and gNaP. The gleak value for each neuron was randomly drawn from a mixture of three Gaussians with weights [0.35, 0.1, 0.55], means [0.5, 0.7, 1.2] nS, and standard deviation 0.05 nS. The gNaP values are drawn from a Gaussian with mean 0.8 nS and standard deviation 0.05 nS. Classification of intrinsic activity is done using peak detection on the voltage V recorded with synapses blocked.

MOR targeting

In all simulations, half of the excitatory neurons are opioid-sensitive (MOR+) and can be targeted with opioid (Dop = 1), while the inhibitory neurons and the other half of the excitatory neurons (MOR−) are insensitive (Dop = 0). The opioid-sensitive population’s excitatory synapses follow the Isyn,op equation, whereas the insensitive neurons follow Isyn,exc. Assignment of Dop among excitatory neurons is random except in two cases shown in Figure 5, where opioid is applied to the excitatory neurons with gleak values below or above the median among the excitatory population.

Gradual ramp up of opioid

For opioid ramping simulations (Figs. 1, 3–5), Ihyp,op ramped from 0 to 8 pA, increasing by 0.5% every 3 s, while gsyn,op gradually decreased from 1.0 to 0.0 nS by 0.5% every 3 s. Hence, the total length of the simulation is 10 min. The opioid shutdown dosage was calculated by averaging the Ihyp,op values at the time of the last bursts with amplitudes of 10–15 Hz.

Figure 1.
Figure 1.

Opioid sensitivity varies across model preBötC networks. A, A network graph of the in silico preBötC network. B, Phase diagrams showing intrinsic activities of each neuron (open circles) based on gleak and gNaP conductances. Top-row: MOR− neurons. Bottom-left: control condition, MOR+ neurons. Bottom-right: opioid applied to MOR+ neurons. C, Quantified changes in the number of neurons with silent, bursting, or tonic intrinsic activities in response to opioid (n = 40 networks; two-tailed paired t-tests; ****p < 0.0001). D, Traces of four network simulations where opioid is ramped up. Numbered boxes show the last bursts detected at a given amplitude threshold (10–15 Hz/cell). E, Histogram and kernel density estimation of the distribution of opioid shutdown doses for n = 40 model networks.

Figure 2.
Figure 2.

Intrinsic cellular activities do not predict opioid sensitivity. A, Example rhythms (top) and overlaid burst waveforms (bottom) under control conditions and in the presence of opioid from representative “high-sensitivity” (left) and “low-sensitivity” (right) networks. B, Phase diagrams of high (left) and low (right) sensitivity networks, showing intrinsic activities of MOR− (top) and MOR+ (bottom) neurons (open circles) based on gleak and gNaP conductances. C, Quantified relationship between opioid shutdown and the number of silent, bursting, and tonic neurons under control conditions (top) and in the presence of opioid (bottom) (n = 40 networks; two-tailed paired t-tests; ns = not significant).

Figure 3.
Figure 3.

Increased connection density reduces opioid sensitivity. A, Example traces of four different simulations with varied connection densities where opioid is ramped up (opioid = 0–8) over 10 min. B, Kernel density estimations showing the distribution of shutdown dosages based on connection probabilities. C, Quantified opioid shutdown dose versus connection probability (n = 40 networks; one-way repeated measures analysis of variance (RM ANOVA) with Bonferroni multiple comparisons tests; *p < 0.05, **p < 0.01, ****p < 0.0001).

Figure 4.
Figure 4.

Network structure regulates opioid sensitivity. Correlation analysis of the relationship between opioid shutdown dose and connectivity within and between A, excitatory and inhibitory populations, B, MOR+, MOR−, and inhibitory populations, and C, tonic, bursting, and silent excitatory and inhibitory subpopulations (n = 40, two-tailed correlation analysis; *p < 0.05, **p < 0.01, ***p < 0.001). Numbers in A and B represent the max and min number of each type of connection.

Figure 5.
Figure 5.

Identity of MOR+ neurons regulates opioid sensitivity. A, Example rhythms (top) and burst waveforms (bottom) in response to opioids when MOR is assigned randomly (left) or specifically to low gleak (middle) or high gleak (right) populations. B, Intrinsic activities of MOR− and MOR+ neurons (open circles) in gNaP and gleak space of the example networks shown in A. C, Quantified number of silent, bursting, and tonic neurons under control conditions and in response to opioid when MOR is assigned randomly or to low/high gleak populations (n = 40 each, one-way ANOVA with Bonferroni multiple comparisons tests, ns = not significant, *p < 0.05, **p < 0.01, ****p < 0.0001). D, Example network rhythm during opioid ramp (opioid = 0–8) with MOR assigned randomly or to low/high gleak populations. E, Kernel density estimations showing distributions of opioid shutdown dosages based on the identity of MOR expressing neurons (n = 40, one-way ANOVA with Bonferroni multiple comparisons tests; *p < 0.05, ****p < 0.0001).

Timed all-or-nothing perturbations

In simulations with timed all-or-nothing perturbations (Figs. 2, 5A–C, 6, 7), we allowed for a 10- s transient period before each perturbation. Data from transient periods are not used in our analysis. When the opioid perturbation is turned on, Ihyp,op = 4 pA and gsyn,opioid = 0.5 nS. We varied gNaP (Fig. 6) or gleak (Fig. 7). For each node, gNaP was increased to 110%, 130%, and 150% of control values, whereas gleak was decreased to 90%, 70%, and 50% of control values. In Figures 6 and 7, the 200 s experimental procedure is as follows:

  1. 10 s transient period

  2. 30 s control period

  3. 10 s transient period

  4. 30 s opioid perturbation

  5. 10 s transient period

  6. 30 s control “wash” period

  7. 10 s transient period

  8. 30 s gNaP or gleak perturbation

  9. 10 s transient period

  10. 30 s simultaneous perturbation of opioid and gNaP or gleak

Figure 6.
Figure 6.

Modulation of gNaP renders the network resistant to opioids. A, Example rhythm and burst waveforms from a network (MOR randomly assigned) in response to opioid and during concurrent modulation of gNaP to 110% , 130% , and 150% of control values. B, Quantified effects on frequency (top) and burst amplitude (bottom) during opioid and gNaP modulation (n = 40, one-way RM ANOVA with Bonferroni multiple comparisons tests, **p < 0.01, ***p < 0.001, ****p < 0.0001). C, Changes in the intrinsic activities in gNaP and gleak space of MOR− and MOR+ neurons from the example network shown in A. D, Quantified changes in the number of silent, bursting, and tonic neurons in response to opioid and subsequent modulation of gNaP (n = 40, one-way RM ANOVA with Bonferroni multiple comparisons tests, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).

Figure 7.
Figure 7.

Modulation of gleak renders the network resistant to opioids. A, Example rhythm and burst waveforms from a network in response to opioid and during concurrent modulation of gleak to 90% , 70% , and 50% of control values. B, Quantified effects on frequency (top) and burst amplitude (bottom) during opioid and gleak modulation (n = 40, one-way RM ANOVA with Bonferroni multiple comparisons tests, **p < 0.01, ****p < 0.0001). C, Changes in the intrinsic activities in gNaP and gleak space of MOR− and MOR+ neurons from the example network shown in A. D, Quantified changes in the number of silent, bursting, and tonic neurons in response to opioid and subsequent modulation of gleak (n = 40, one-way RM ANOVA with Bonferroni multiple comparisons tests, ****p < 0.0001).

Analysis
Burst detection and opioid shutdown dosage

Bursts were detected using a basic peak-finding algorithm [find_peaks function in scipy (Virtanen et al., 2020)] where each peak must have a minimum height of 4 Hz/cell and minimum prominence of 10 Hz/cell. We then compute the opioid shutdown dosage by finding the averaging Iopioid values at the time of the last bursts meeting an amplitude threshold of 10, 11, 12, 13, 14, and 15 Hz/cell. Statistical analysis of measured variables was performed using GraphPad Prism 10 software, and data were visualized using a combination of python, GraphPad Prism, and Powerpoint.

Phase diagrams

Each phase boundary was computed by simulating the network under synaptic block, sweeping across a grid of conductances gleak ∈ {0.2, 0.3, …, 1.5} nS and gNaP = {0.6, 0.7, …, 1.5} nS. The points plotted within the phase boundaries represent the neurons in the two-population network simulated under synaptic block. Population neurons and phase sections are both colored by intrinsic activity classified as described in “Burst detection and opioid shutdown dosage”.

Code availability

Code is available in our github repository or upon request to the corresponding author. Code was run using Linux systems with python simulation stack.

Code files

Download Code files, ZIP file.

Statistics Summary

Download Statistics Summary, XLS file.

Results

Opioid sensitivity varies across model preBötC networks

In sparsely connected (1% ) preBötC networks, connections were drawn randomly between excitatory and inhibitory neurons with different intrinsic activities as determined by varied gNaP and gleak values. We implemented a two-population distribution of gleak and gNaP (Baertsch et al., 2021) as described in Methods, Network construction to reduce the number of neurons in the network that exhibit intrinsic bursting to 5–10% (Fig. 1), as estimated by experimental recordings in vitro. With 50% of excitatory neurons randomly designated as opioid-responsive (MOR+) and the remaining 50% non-responsive (MOR−), we performed simulations on 40 different preBötC networks where the effects of opioid (i.e., presynaptic suppression and hyperpolarization) were gradually increased over the course of 10  min. An example of how simulated opioid affected the intrinsic activities of MOR+ and MOR− neurons in gleak and gNaP parameter space is shown in Figure 1B. Across all networks, opioids transitioned the intrinsic activity of preBötC neurons from tonic or bursting to silent (Fig. 1C), similar to observations in preBötC slices (Burgraff et al., 2021). Some networks were highly sensitive to opioid, as shown by a quicker decline in the respiratory rhythm (e.g., Fig. 1D, traces 1 and 2), while other networks were quite resistant (e.g., Fig. 1D, traces 3 and 4). This variability is reflected in the distribution of shutdown dosages, which ranged from 3.73 to 7.51 with a mean of 5.26 and was approximately Gaussian (Fig. 1E). Thus, despite all simulated networks having the same number of neurons designated as excitatory, inhibitory, and opioid sensitive (MOR+), there was considerable variation in how individual networks responded to opioids.

Changes in intrinsic cell activities do not predict opioid sensitivity

To explore how random differences in the intrinsic cellular activities of the networks may predict the varied responses of the network rhythm to opioids, we compared networks with high and low opioid sensitivity. “High-sensitivity” networks were defined as those with an above-median opioid shutdown dosage, while networks with a below-median shutdown dosage were considered “low-sensitivity”. Rather than the gradual opioid ramping as shown in Figure 1, in Figure 2 we instead simulated a 30-s control period followed by a 30-s period with a moderate dose of opioid applied (opioid = 4). The variation in opioid sensitivity is exemplified in Figure 2A, where we see clear differences in how the rhythm responded to opioid. In the high-sensitivity case, the rhythm became weak and irregular, whereas rhythms produced by the most resistant networks were able to maintain consistent frequencies and burst amplitudes close to baseline. Changes in the intrinsic cellular activities of these representative high- and low-sensitivity networks are shown in gleak and gNaP parameter space in Figure 2B. Under control conditions and in the presence of opioid, the proportions of neurons with silent, bursting, or tonic intrinsic activity were similar between high- and low-sensitivity networks (Fig. 2C). Indeed, regardless of opioid sensitivity, a similar number of MOR+ neurons that were tonic or bursting in control conditions became silent in the presence of opioid, which was consistent across all 40 networks. Thus, differences in how opioids affect the intrinsic activities of neurons in our model networks are unlikely to explain their variable responses to opioids.

Connection density and network structure regulate opioid sensitivity

To test how the total amount of connectivity with the preBötC model networks affects how they respond to opioids, we ran simulations where the connection probability of each neuron was increased to 2% , 4% , 8% , or 16% for 40 networks each (the default for all other experiments is 1% connection density), while maintaining total synaptic strength in the network constant. The results are shown in Figure 3. For each trace in Figure 3A, we can see that networks with higher connectivity are able to maintain a network rhythm at higher doses of opioid. The distribution of dosages that effectively shut down each network also tends to be slightly less variable at higher connection probabilities (Fig. 3B and C). Thus, preBötC networks with higher total connection densities are more resistant to opioids.

Next, we examined how random differences in connection topology may contribute to the variation in opioid responses observed across our 40 randomly drawn model networks. To do so, we first tested whether the total number of excitatory and inhibitory connections (excitation/inhibition balance) within each model network was related to its sensitivity to opioids (Fig. 4A). Correlation analysis revealed that, in general, networks with a more highly connected excitatory population and fewer inhibitory inputs to these excitatory neurons were more resistant to opioids (i.e., higher opioid shutdown dose). In contrast, overall connectivity within the inhibitory population or from excitatory to inhibitory neurons was not correlated with the sensitivity of the network rhythm to opioids. Next, we tested more specifically whether the number of connections within and between excitatory MOR+, excitatory MOR−, and inhibitory neurons was correlated with the opioid dose that shutdown rhythm generation (Fig. 4B). We found that when the population of excitatory MOR− neurons was more interconnected and received less inhibitory input, the network was more likely to be resistant to opioids.

In a third analysis, we broke the network connectivity down even further by computing correlations between opioid shutdown dose and the number of connections among intrinsically tonic, bursting, and silent excitatory MOR+, excitatory MOR−, and inhibitory neuron subpopulations (Fig. 4C). This revealed three primary observations. First, the number of connections from silent to tonic excitatory MOR− neurons was the strongest driver of opioid resistance among this MOR− population. Second, although the total number of connections within the MOR+ population was not predictive of opioid sensitivity, networks with more connections between intrinsically tonic MOR+ neurons and fewer connections between intrinsically silent MOR+ neurons were more resistant to opioids. And third, networks were also more likely to be resistant to opioids if they had more connections from tonic MOR+ neurons to tonic or silent MOR− neurons and fewer connections from bursting MOR+ neurons to silent MOR− neurons. Overall, these correlation analyses suggest that differences in network topology as a result of randomness in the assignment of network connections contribute to the variable responses of preBötC networks to opioids.

Identity of MOR+ neurons regulates opioid sensitivity

Because 50% of the excitatory neurons in our model networks are randomly designated as MOR+, we next wondered how the opioid sensitivity of the model networks may be altered if the identity of MOR+ neurons is non-random. To address this question, we performed simulations to compare opioid responses in networks where the intrinsically silent neurons (high gleak) or the tonic/bursting neurons (low gleak) were designated as MOR+, as described in “MOR targeting”. Example network activity during these experiments is shown in Figure 5A. Compared to random assignment of MOR as described above, assigning MOR to the low gleak population made the rhythms more resistant to opioids, whereas assigning MOR to the high gleak population made them more sensitive. When low gleak (primarily intrinsically tonic/bursting) neurons are MOR+, 92.3% of the population became intrinsically silent in response to opioids (Fig. 5B). However, this only reflects the intrinsic activity of those neurons with synapses blocked. With synaptic interactions intact, the network remained rhythmic, with a slower frequency than control conditions but a similar amplitude. On the other hand, when the high gleak (primarily intrinsically silent) population is MOR+, the rhythm collapsed under only a moderate dose of opioid (Ihyp,op = 4 pA) (Fig. 5A). In this case, changes in the intrinsic activities of neurons in the network in response to opioid were minimal (Fig. 5B).

The above results were for a single exemplar network. In Figure 5C, for each condition, we compared the number of intrinsically silent, bursting, and tonic neurons and how the distributions of these intrinsic activities change in response to opioids across 40 different model networks. As expected (see Fig. 1B), when the identity of MOR+ neurons was randomly assigned, opioids caused many of the low gleak MOR+ neurons to transition from tonic/bursting activity to silent, whereas high gleak MOR+ neurons were largely unaffected. Under non-random conditions, when all low gleak neurons were designated as MOR+, changes in the intrinsic activities within the network were exaggerated such that nearly all intrinsically tonic activity was lost as 92% of the network became intrinsically silent. In contrast, when high gleak neurons were designated as MOR+, there were minimal changes in the distribution of intrinsic activities within the networks (Fig. 5C). To further test how the identity of MOR+ neurons may alter how the preBötC network rhythm responds to opioids, we performed simulations ramping up the opioid effect to compare shutdown dosages for networks with MOR identity assigned randomly, or selectively to low or high gleak populations (Fig. 5D and E). Notably, despite a much larger proportion of the network becoming intrinsically silent, networks with low gleak neurons designated as MOR+ were more resistant to opioids, than when MOR identity was randomly assigned. On the other hand, the average shutdown dosage was lower when high gleak neurons were designated as MOR+, indicating that, despite the minimal effects on the intrinsic activities of the neurons, the network rhythm was substantially more sensitive to opioids under these conditions. These findings support the conclusion that changes in intrinsic cellular activities within the network are not predictive of its sensitivity to opioids (see Fig. 2), but that the distribution of MOR+ expression among preBötC neurons may be an important determinant of how the network responds to opioids.

Modulation of gNaP or gleak can render the preBötC resistant to opioids

Considering these results, we tested whether manipulations of the intrinsic properties of preBötC neurons may represent a viable strategy to protect the preBötC rhythm from the effects of opioids. Specifically, we tested whether increasing gNaP would allow for sustained rhythmogenesis in the presence of relatively high opioid doses as previously hypothesized based on pharmacological experiments in vitro (Burgraff et al., 2021). We also tested whether decreasing the leak conductance gleak would have a similar protective effect on rhythmogenesis. Rhythmic activity of a representative network under control conditions, in opioid, and during a subsequent 10% , 30% , and 50% increase in gNaP is shown in Figure 6A. Increasing gNaP by 30% in the model networks reversed the effects of opioids on burst frequency and amplitude (Fig. 6B). However, recovery of the rhythm by gNaP modulation did not restore intrinsic cellular activities to near control. Instead, it was associated with a change in the intrinsic activity of both MOR+ and MOR− neurons from silent to bursting, with little effect on the number of tonic neurons (Fig. 6C and D). Under control conditions, the network was composed of mostly intrinsically tonic and silent neurons (52.3% and 40.7% , respectively). In response to opioid, the proportion of silent neurons increased to 73.7% as MOR+ neurons transitioned from tonic to silent. As gNaP was increased, the MOR+ neurons that were originally tonic under control conditions transitioned to bursting. Specifically, when gNaP was increased by 30% , 55% of the population entered a gNaP and gleak parameter space that supports intrinsic bursting. Thus, despite recovery of a rhythm with similar frequency and amplitude characteristics following gNaP modulation, the number of intrinsically tonic neurons remained reduced, whereas the number of bursting neurons was increased relative to control conditions.

We next performed similar simulations during manipulation of gleak (Fig. 7). The rhythmic activity of a representative network under control conditions, in opioid, and following a subsequent 10% , 30% , and 50% reduction in gleak is shown in Figure 7A. In this case, burst amplitude but not frequency could be significantly recovered toward control values (Fig. 7B). This was associated with changes in the intrinsic activities of primarily MOR− neurons (Fig. 7C). When gleak was reduced to 70% of control, there was a large increase in the number of bursting neurons, and upon further reduction of gleak to 50% of control, these neurons became tonic, leaving only 2.7% of the population as bursting (Fig. 7D). Thus, our model predicts that manipulations that directly or indirectly affect persistent sodium and/or potassium leak conductances may be effective for increasing the resistance of preBötC function to opioids.

Discussion

The effect of opioids on respiratory function is variable in brain slices in vitro, animal models in vivo, and in individual humans (Cherny et al., 2001; Dahan et al., 2005, 2013; Burgraff et al., 2021). Here we adopt a computational model of the respiratory rhythm generator to dissect plausible network topology and cellular properties that contribute to variable respiratory responses to opioids. We leverage computational models that allow us to instantiate networks of the preBötC with connectivity patterns and conductances drawn from random distributions. These networks are statistically indistinguishable on the “macro”-scale; they have the same overall numbers of excitatory and inhibitory neurons, the same numbers of MOR+ and MOR− neurons, the same probability of connections per neuron, and conductance values are drawn from the same distributions. Yet, due to the random assignment of some of these properties, each network differs on the level of individual neurons (nodes), which vary in their exact connectivity patterns and conductance strengths. Surprisingly, this “micro”-level randomness is sufficient to create quite variable responses at the network level to the same stimulus—in this case simulated opioids. We suspect that these differences may contribute to the observed variable responses to opioids seen in experimental preparations (Burgraff et al., 2021). Furthermore, this micro-level variability could, for example, explain how individuals may respond differently to network perturbations despite the preBötC network developing with the same general set of instructions (e.g., genome, transcriptome, axonal targeting mechanisms, etc.). While OIRD arises from the effects of opioids on multiple central and peripheral sites (Ramirez et al., 2021), our simulations illustrate how variation in the architecture of the inspiratory rhythm generator could be an important factor underlying the unpredictability of opioid overdose.

The computational approach here allows for directed manipulations that are experimentally intractable. For instance, we are able to ask if the response of the preBötC to opioids depends on MOR being expressed in populations with particular conductance profiles. More concretely, we target the opioid effect directly to neurons that have a particular leak conductance. This leak conductance (gleak) is an important determinant of whether a neuron is intrinsically “tonic”, “bursting”, or “silent” (Butera et al., 1999; Del Negro et al., 2002; Koizumi and Smith, 2008; Yamanishi et al., 2018). Surprisingly, introducing MOR selectively to low gleak (intrinsically excited neurons with tonic/bursting activity) decreased the response of the network to opioids making the rhythm more resilient. Conversely, introducing MOR selectively to the less excitable population (the high gleak, quiescent cells) increased the susceptibility of the network rhythm to opioids. We speculate that a robust preBötC rhythm relies on the existence of a population of “recruitable” neurons that are not strongly intrinsically active, but are capable of becoming active with a small amount of synaptic input. When opioids affect neurons in the low gleak population, their intrinsic activity is reduced but they remain in the recruitable pool and therefore can continue to participate in the network, allowing the rhythm to continue at higher opioid doses. Conversely, we expect that when opioids further suppress neurons that already have low intrinsic excitability (high gleak), they are removed from the recruitable pool and unable to participate in network bursts, making coordinated network activity more vulnerable to opioids. When the effect of opioids is randomly targeted to 50% of neurons, the proportion that remains recruitable in the presence of opioid depends on how MOR expression is randomly assigned within the high and low gleak populations, contributing to variable opioid responses at the network level.

Network connectivity is difficult to study and manipulate experimentally. Thus, computational models, where the number and strength of all connections between every neuron are known, can be an important tool to provide “proof of concept” insights into how network topology can influence network function and determine its response to perturbations. We took advantage of this by performing correlation analysis to better understand how the number of connections between certain subgroups of preBötC neurons may predict how susceptible the network is to opioids. These analyses revealed that, in general, when neurons that do not respond to opioid (MOR−) are more interconnected and receive less inhibitory input, the network is more resistant to opioids. We suspect that this connectivity configuration may allow the network of MOR− neurons to remain rhythmogenic even when very few opioid-sensitive (MOR+) neurons are able to contribute to network function. In another analysis, we scaled the number of connections within the network without altering total synaptic strength, which consistently increased the robustness of the network to opioids. Because opioids weaken the presynaptic strength of excitatory interactions (Baertsch et al., 2021), we anticipate that networks with lower numbers of connections become “fractured” into isolated sub-networks when opioid-induced weakening of synapses impairs the network’s ability to effectively recruit portions of the population. Indeed, the preBötC rhythm in vitro has a higher proportion of failed bursts with low amplitude in response to opioids (Baertsch et al., 2021; Phillips and Rubin, 2022). In networks with more connections, activity more consistently propagates to all neurons (Kam et al., 2013), efficiently recruiting the whole population despite the effect of opioids on synaptic transmission. This could also contribute to the variable opioid responses observed in in vitro experiments since both within and across labs where the creation of rhythmic brain stem slices invariably samples slightly different portions of the preBötC population that may be more or less densely connected (Ruangkittisakul et al., 2014; Baertsch et al., 2019). Although these simulations illustrate that network topology could be an important determinant of opioid sensitivity, because connection density and patterns are considered “fixed” properties of the network, at least on short time scales, manipulation of network topology is an unlikely avenue for therapeutic interventions. In contrast, the strength of existing excitatory synaptic connections can be pharmacologically altered acutely via, e.g., ampakines, which may render the preBötC less vulnerable to opioids and shows promise as an intervention for OIRD (Ren et al., 2006; Xiao et al., 2020; Sunshine and Fuller, 2021).

The intrinsic activity of preBötC neurons is determined by multiple interacting cellular properties (Ramirez et al., 2012). Not all are known and not all can be incorporated into our simplified model network. Yet, like many other computational studies (Lindsey et al., 2012), the interaction between gleak and gNaP determines intrinsic activity in our model and is sufficient to capture the silent, bursting, or tonic phenotypes of preBötC neurons. Both gleak and gNaP contribute to cellular excitability (resting membrane potential), and the voltage-dependent properties of gNaP allow some neurons with appropriate gleak to exhibit intrinsic bursting or “pacemaker” activity (Koizumi and Smith, 2008). Whether such neurons with intrinsic bursting capabilities have a specialized role in network rhythmogenesis is a matter of ongoing debate (Smith et al., 2000; Feldman and Del Negro, 2006; Ramirez and Baertsch, 2018a,b; da Silva et al., 2023) that we do not address here. Instead, we aimed to understand how opioids alter the intrinsic activities of preBötC neurons. In the model network, opioids reduce the number of neurons with intrinsic bursting or tonic activity and increase the number of silent neurons. To our surprise, the extent of these changes was not a significant predictor of the network response to opioids. This suggests that the intrinsic activity of a given neuron may not be representative of its contribution to network function, and that other factors, such as those discussed above, play more substantial roles in determining how the preBötC responds to opioids. Although network differences due to random sampling of gleak and gNaP from set distributions were not a significant factor driving variable opioid responses, we found that scaling the distribution of gNaP or gleak across the whole population did alter the sensitivity of model networks to opioids. Interestingly, manipulation of gNaP was more effective since a 30% increase in gNaP was sufficient to restore both frequency and amplitude of the rhythm, whereas effects were more specific to burst amplitude following a 30% decrease in gleak. Unlike network topology, intrinsic conductances that regulate cellular excitability and activity are not “fixed” but are dynamic and can be modified by conditional changes in, e.g., neuromodulators and ion concentrations (Rybak et al., 2007; Ramirez et al., 2012) and are also more amenable to pharmacological manipulations (Bedoya et al., 2019; Verneuil et al., 2020; Burgraff et al., 2021). Thus, further experimental investigation of these approaches is warranted as they may hold promise as potential therapeutic strategies to protect against opioid-induced failure of preBötC network function.

Footnotes

  • The authors declare no competing financial interests.

  • We thank our funding sources NIH R00-HL145004 (N.A.B.), NIH R01-HL166317 (N.A.B.), K01-DA058543 (R.S.P.), F32-HL159904 (N.E.B.), and Western Washington University (K.D.H., G.M.C.) for supporting this work.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.

References

  1. ↵
    1. Bachmutsky I,
    2. Wei XP,
    3. Kish E,
    4. Yackle K
    (2020) Opioids depress breathing through two small brainstem sites. Elife 9:e52694.
    OpenUrlCrossRef
  2. ↵
    1. Baertsch NA,
    2. Baertsch HC,
    3. Ramirez JM
    (2018) The interdependence of excitation and inhibition for the control of dynamic breathing rhythms. Nat Commun 9:843–843.
    OpenUrl
  3. ↵
    1. Baertsch NA,
    2. Bush NE,
    3. Burgraff NJ,
    4. Ramirez J-M
    (2021) Dual mechanisms of opioid-induced respiratory depression in the inspiratory rhythm-generating network. Elife 10:e67523.
    OpenUrl
  4. ↵
    1. Baertsch NA,
    2. Severs LJ,
    3. Anderson TM,
    4. Ramirez JM
    (2019) A spatially dynamic network underlies the generation of inspiratory behaviors. Proc Natl Acad Sci U S A 116:7493–7502.
    OpenUrlAbstract/FREE Full Text
  5. ↵
    1. Baldo BA,
    2. Rose MA
    (2022) Mechanisms of opioid-induced respiratory depression. Arch Toxicol 96:2247–2260.
    OpenUrl
  6. ↵
    1. Bedoya M,
    2. Rinné S,
    3. Kiper AK,
    4. Decher N,
    5. Gonzaález W,
    6. Ramírez D
    (2019) Task channels pharmacology: new challenges in drug design. J Med Chem 62:10044–10058.
    OpenUrl
  7. ↵
    1. Burgraff NJ,
    2. Bush NE,
    3. Ramirez JM,
    4. Baertsch NA
    (2021) Dynamic rhythmogenic network states drive differential opioid responses in the in vitro respiratory network. J Neurosci 41:9919–9931.
    OpenUrlAbstract/FREE Full Text
  8. ↵
    1. Butera RJ,
    2. Rinzel J,
    3. Smith JC
    (1999) Models of respiratory rhythm generation in the pre-Botzinger complex. I. Bursting pacemaker neurons. J Neurophysiol 82:382–397.
    OpenUrlCrossRefPubMed
  9. ↵
    1. Butera RJ,
    2. Smith JC
    (1999) Models of respiratory rhythm generation in the pre-Botzinger complex. II. Populations of coupled pacemaker neurons. J Neurophysiol 82:398–415.
    OpenUrlCrossRefPubMed
  10. ↵
    1. Cherny N,
    2. Ripamonti C,
    3. Pereira J,
    4. Davis C,
    5. Fallon M,
    6. McQuay H,
    7. Mercadante S,
    8. Pasternak G,
    9. Ventafridda V, Expert Working Group of the European Association of Palliative Care Network
    (2001) Strategies to manage the adverse effects of oral morphine: an evidence-based report. J Clin Oncol 19:2542–2554.
    OpenUrlAbstract/FREE Full Text
  11. ↵
    1. Dahan A,
    2. Overdyk F,
    3. Smith T,
    4. Aarts L,
    5. Niesters M
    (2013) Pharmacovigilance: a review of opioid-induced respiratory depression in chronic pain patients. Pain Phys 16:E85.
    OpenUrl
  12. ↵
    1. Dahan A,
    2. Sarton E,
    3. Teppema L,
    4. Olievier C,
    5. Nieuwenhuijs D,
    6. Matthes HW,
    7. Kieffer BL
    (2001) Anesthetic potency and influence of morphine and sevoflurane on respiration in μ-opioid receptor knockout mice. J Am Soc Anesthesiol 94:824–832.
    OpenUrl
  13. ↵
    1. Dahan A,
    2. Yassen A,
    3. Bijl H,
    4. Romberg R,
    5. Sarton E,
    6. Teppema L,
    7. Olofsen E,
    8. Danhof M
    (2005) Comparison of the respiratory effects of intravenous buprenorphine and fentanyl in humans and rats. Br J Anaesth 94:825–834.
    OpenUrlCrossRefPubMed
  14. ↵
    1. da Silva CA,
    2. Grover CJ,
    3. Picardo MCD,
    4. Del Negro CA
    (2023) Role of NaV1.6-mediated persistent sodium current and bursting-pacemaker properties in breathing rhythm generation. Cell Rep 42:113000.
    OpenUrl
  15. ↵
    1. Del Negro CA,
    2. Johnson SM,
    3. Butera RJ,
    4. Smith JC
    (2001) Models of respiratory rhythm generation in the pre-Bötzinger complex. III. Experimental tests of model predictions. J Neurophysiol 86:59–74.
    OpenUrlCrossRefPubMed
  16. ↵
    1. Del Negro CA,
    2. Koshiya N,
    3. Butera Jr RJ,
    4. Smith JC
    (2002) Persistent sodium current, membrane properties and bursting behavior of pre-Botzinger complex inspiratory neurons in vitro. J Neurophysiol 88:2242–2250.
    OpenUrlCrossRefPubMed
  17. ↵
    1. Doi A,
    2. Ramirez JM
    (2008) Neuromodulation and the orchestration of the respiratory rhythm. Respir Physiol Neurobiol 164:96–104.
    OpenUrlCrossRefPubMed
  18. ↵
    1. Erbs E
    , et al. (2015) A mu–delta opioid receptor brain atlas reveals neuronal co-occurrence in subcortical networks. Brain Struct Funct 220:677–702.
    OpenUrlCrossRefPubMed
  19. ↵
    1. Feldman JL,
    2. Del Negro CA
    (2006) Looking for inspiration: new perspectives on respiratory rhythm. Nat Rev Neurosci 7:232–242.
    OpenUrlCrossRefPubMed
  20. ↵
    1. Fleming E,
    2. Voscopoulos C,
    3. George E
    (2015) Non-invasive respiratory volume monitoring identifies opioid-induced respiratory depression in an orthopedic surgery patient with diagnosed obstructive sleep apnea: a case report. J Med Case Rep 9:1–5.
    OpenUrl
  21. ↵
    1. Freye E,
    2. Latasch L
    (2003) Development of opioid tolerance–molecular mechanisms and clinical consequences. Anasthesiol Intensivmed Notfallmed Schmerzther AINS 38:14–26.
    OpenUrl
  22. ↵
    1. Funk GD,
    2. Smith JC,
    3. Feldman JL
    (1993) Generation and transmission of respiratory oscillations in medullary slices: role of excitatory amino acids. J Neurophysiol 70:1497–1515.
    OpenUrlCrossRefPubMed
  23. ↵
    1. Gray PA,
    2. Rekling JC,
    3. Bocchiaro CM,
    4. Feldman JL
    (1999) Modulation of respiratory frequency by peptidergic input to rhythmogenic neurons in the preBotzinger complex. Science 286:1566–1568.
    OpenUrlAbstract/FREE Full Text
  24. ↵
    1. Greer J,
    2. Smith J,
    3. Feldman J
    (1991) Role of excitatory amino acids in the generation and transmission of respiratory drive in neonatal rat. J Physiol 437:727–749.
    OpenUrlCrossRefPubMed
  25. ↵
    1. Harris KD,
    2. Dashevskiy T,
    3. Mendoza J,
    4. Garcia AJ,
    5. Shea-Brown E
    (2017) Different roles for inhibition in the rhythm-generating respiratory network. J Neurophysiol 118:2070–2088.
    OpenUrlCrossRefPubMed
  26. ↵
    1. Kam K,
    2. Worrell JW,
    3. Janczewski WA,
    4. Cui Y,
    5. Feldman JL
    (2013) Distinct inspiratory rhythm and pattern generating mechanisms in the preBotzinger complex. J Neurosci 33:9235–9245.
    OpenUrlAbstract/FREE Full Text
  27. ↵
    1. Koizumi H,
    2. Smith JC
    (2008) Persistent Na+ and K+-dominated leak currents contribute to respiratory rhythm generation in the pre-Botzinger complex in vitro. J Neurosci 28:1773–1785.
    OpenUrlAbstract/FREE Full Text
  28. ↵
    1. Lindsey BG,
    2. Rybak IA,
    3. Smith JC
    (2012) Computational models and emergent properties of respiratory neural networks. Compr Physiol 2:1619.
    OpenUrl
  29. ↵
    1. Lynch N,
    2. Lima JD,
    3. Spinieli RL,
    4. Kaur S
    (2023) Opioids, sleep, analgesia and respiratory depression: their convergence on mu (μ)-opioid receptors in the parabrachial area. Front Neurosci 17:1134842.
    OpenUrl
  30. ↵
    1. Phillips RS,
    2. Rubin JE
    (2019) Effects of persistent sodium current blockade in respiratory circuits depend on the pharmacological mechanism of action and network dynamics. PLoS Comput Biol 15:e1006938.
    OpenUrlCrossRefPubMed
  31. ↵
    1. Phillips RS,
    2. Rubin JE
    (2022) Putting the theory into “burstlet theory” with a biophysical model of burstlets and bursts in the respiratory preBötzinger complex. eLife 11:e75713. doi:10.7554/eLife.75713
    OpenUrlCrossRef
  32. ↵
    1. Ramirez J-M,
    2. Baertsch NA
    (2018a) The dynamic basis of respiratory rhythm generation: one breath at a time. Annu Rev Neurosci 41:475–499.
    OpenUrlCrossRefPubMed
  33. ↵
    1. Ramirez JM,
    2. Baertsch NA
    (2018b) The dynamic basis of respiratory rhythm generation: one breath at a time. Annu Rev Neurosci 41:475–499.
    OpenUrlCrossRefPubMed
  34. ↵
    1. Ramirez J-M,
    2. Burgraff NJ,
    3. Wei AD,
    4. Baertsch NA,
    5. Varga AG,
    6. Baghdoyan HA,
    7. Lydic R,
    8. Morris KF,
    9. Bolser DC,
    10. Levitt ES
    (2021) Neuronal mechanisms underlying opioid-induced respiratory depression: our current understanding. J Neurophysiol 125:1899–1919.
    OpenUrlCrossRefPubMed
  35. ↵
    1. Ramirez J,
    2. Doi A,
    3. Garcia III A,
    4. Elsen F,
    5. Koch H,
    6. Wei A
    (2012) The cellular building blocks of breathing. Compr Physiol 2:2683.
    OpenUrl
  36. ↵
    1. Ramirez J-M,
    2. Koch H,
    3. Garcia AJ,
    4. Doi A,
    5. Zanella S
    (2011) The role of spiking and bursting pacemakers in the neuronal control of breathing. J Biol Phys 37:241–261.
    OpenUrlCrossRefPubMed
  37. ↵
    1. Ren J,
    2. Poon BY,
    3. Tang Y,
    4. Funk GD,
    5. Greer JJ
    (2006) Ampakines alleviate respiratory depression in rats. Am J Respir Crit Care Med 174:1384–1391.
    OpenUrlCrossRefPubMed
  38. ↵
    1. Rousseau J-P,
    2. Furdui A,
    3. da Silveira Scarpellini C,
    4. Horner RL,
    5. Montandon G
    (2023) Medullary tachykinin precursor 1 neurons promote rhythmic breathing. Elife 12:e85575.
    OpenUrl
  39. ↵
    1. Ruangkittisakul A,
    2. Kottick A,
    3. Picardo MC,
    4. Ballanyi K,
    5. Del Negro CA
    (2014) Identification of the pre-bötzinger complex inspiratory center in calibrated “sandwich” slices from newborn mice with fluorescent dbx1 interneurons. Physiol Rep 2:e12111.
    OpenUrlAbstract/FREE Full Text
  40. ↵
    1. Rybak IA,
    2. Abdala AP,
    3. Markin SN,
    4. Paton JF,
    5. Smith JC
    (2007) Spatial organization and state-dependent mechanisms for respiratory rhythm and pattern generation. Prog Brain Res 165:201–220.
    OpenUrlCrossRefPubMed
  41. ↵
    1. Sherman D,
    2. Worrell JW,
    3. Cui Y,
    4. Feldman JL
    (2015) Optogenetic perturbation of preBötzinger complex inhibitory neurons modulates respiratory pattern. Nat Neurosci 18:408–416.
    OpenUrlCrossRefPubMed
  42. ↵
    1. Smith JC,
    2. Butera RJ,
    3. Koshiya N,
    4. Del Negro C,
    5. Wilson CG,
    6. Johnson SM
    (2000) Respiratory rhythm generation in neonatal and adult mammals: the hybrid pacemaker-network model. Respir Physiol 122:131–147.
    OpenUrlCrossRefPubMed
  43. ↵
    1. Smith JC,
    2. Ellenberger HH,
    3. Ballanyi K,
    4. Richter DW,
    5. Feldman JL
    (1991) Pre-Bötzinger complex: a brainstem region that may generate respiratory rhythm in mammals. Science 254:726–729.
    OpenUrlAbstract/FREE Full Text
  44. ↵
    1. Sora I,
    2. Funada M,
    3. Uhl GR
    (1997) The μ-opioid receptor is necessary for [D-Pen2, D-Pen5] enkephalin-induced analgesia. Eur J Pharmacol 324:R1–R2.
    OpenUrlCrossRefPubMed
  45. ↵
    1. Stimberg M,
    2. Brette R,
    3. Goodman DF
    (2019) Brian 2, an intuitive and efficient neural simulator. eLife 8:e47314.
    OpenUrlCrossRef
  46. ↵
    1. Sun X,
    2. Pérez CT,
    3. Halemani DN,
    4. Shao XM,
    5. Greenwood M,
    6. Heath S,
    7. Feldman JL,
    8. Kam K
    (2019) Opioids modulate an emergent rhythmogenic process to depress breathing. Elife 8:e50613.
    OpenUrlCrossRef
  47. ↵
    1. Sunshine MD,
    2. Fuller DD
    (2021) Automated classification of whole body plethysmography waveforms to quantify breathing patterns. Front Physiol 12:690265. doi:10.3389/fphys.2021.690265
    OpenUrlCrossRef
  48. ↵
    1. Uniyal A,
    2. Gadepalli A,
    3. Akhilesh ,
    4. Tiwari V
    (2020) Underpinning the neurobiological intricacies associated with opioid tolerance. ACS Chem Neurosci 11:830–839.
    OpenUrl
  49. ↵
    1. Varga AG,
    2. Reid BT,
    3. Kieffer BL,
    4. Levitt ES
    (2020) Differential impact of two critical respiratory centres in opioid-induced respiratory depression in awake mice. J Physiol 598:189–205.
    OpenUrlCrossRefPubMed
  50. ↵
    1. Verneuil J,
    2. Brocard C,
    3. Trouplin V,
    4. Villard L,
    5. Peyronnet-Roux J,
    6. Brocard F
    (2020) The m-current works in tandem with the persistent sodium current to set the speed of locomotion. PLoS Biol 18:e3000738.
    OpenUrlCrossRef
  51. ↵
    1. Virtanen P
    , et al. (2020) SciPy 1.0: fundamental algorithms for scientific computing in python. Nat Methods 17:261–272.
    OpenUrlCrossRefPubMed
  52. ↵
    1. Winter SM,
    2. Fresemann J,
    3. Schnell C,
    4. Oku Y,
    5. Hirrlinger J,
    6. Hülsmann S
    (2009) Glycinergic interneurons are functionally integrated into the inspiratory network of mouse medullary slices. Pflugers Archiv Eur J Physiol 458:459–469.
    OpenUrlCrossRefPubMed
  53. ↵
    1. Xiao D,
    2. Xie F,
    3. Xu X,
    4. Zhou X
    (2020) The impact and mechanism of ampakine CX1739 on protection against respiratory depression in rats. Future Med Chem 12:2093–2104.
    OpenUrl
  54. ↵
    1. Yamanishi T,
    2. Koizumi H,
    3. Navarro MA,
    4. Milescu LS,
    5. Smith JC
    (2018) Kinetic properties of persistent Na+ current orchestrate oscillatory bursting in respiratory neurons. J Gen Physiol 150:1523–1540.
    OpenUrlAbstract/FREE Full Text

Synthesis

Reviewing Editor: Leandro Vendruscolo, National Institute on Drug Abuse

Decisions are customarily a result of the Reviewing Editor and the peer reviewers coming together and discussing their recommendations until a consensus is reached. When revisions are invited, a fact-based synthesis statement explaining their decision and outlining what is needed to prepare a revision will be listed below. The following reviewer(s) agreed to reveal their identity: NONE.

Reviewer 1:

The step-by-step increments on the complexity of the manipulation, building on the previous findings made it easy to follow and a good read.

Reviewer 2:

This manuscript takes an in silico approach to try to understand the complexities underlying variability in responses of breathing networks, particularly in the preBotC, to opioids. Models with random sparsely connected networks of 300 neurons with varying levels of intrinsic activity were generated and found to have low and high levels of sensitivity to simulated opioid. Opioid sensitivity was not predicted by intrinsic activity or the number of MOR+ neurons silenced by opioids. Higher network connectivity, especially between MOR- excitatory neurons, was protective from opioid. Surprisingly, high gLeak in MOR+ neurons made networks more susceptible to simulated opioid, even though this rendered most MOR+ neurons silent in the absence of simulated opioid. Opioid effects on network output were reversed by increasing gNaP and decreasing gLeak. This study is an appropriate use of in silico modeling to answer questions that are not accessible by other techniques. The simulated data are generally well presented and accessible, and the findings are of interest to the control of breathing field. The Discussion nicely highlights the advantages and limitations of the study, and provides physiological relevance without overstating.

I have a few questions regarding the design of the model, the parameters that were used to simulate opioid exposure and some minor comments.

1) I thought there are equal numbers of inhibitory and excitatory neurons in the preBotC. So, why does the model have only 60 inhibitory neurons compared to 240 excitatory neurons?

2) If I am reading the figures correctly, the distribution of intrinsic properties (gLeak and gNaP) is not random along a continuum, but clustered into two groups. What is the basis of these clusters?

3) There is no attempt to classify non-silent model neurons as inspiratory or expiratory, since MORs are expressed in both populations. Why? In Baertsch, 2021 all in silico neurons were inspiratory related. Assuming the same is true here, it is not really reflective of the physiology. Please comment.

4) Figure 1B and 2B: There's no change in the MOR- neurons after application of the opioid condition? This seems non physiological. Wouldn't some of these MOR- neurons receive input from MOR+ synapses, or at least the condition of the network under the opioid perturbation would also alter the activity of MOR- neurons? In Baertsch, 2021, the activity of some MOR- neurons is also affected by DAMGO and simulated opioid. Is the readout of "silent" too binary in this study to detect such changes? Please comment.

5) In the case where MOR+ neurons are assigned low gLeak and MOR- neurons are therefore high gLeak, nearly all nodes are silent during opioid. Where is the evidence that these nodes are "recruitable", as speculated by the authors? Wouldn't this be detected by the model and the nodes would not be classified as "silent"? Rather it seems a very small number of MOR- nodes are bursting and able to generate opioid resistant output. Please comment.

6) Page 13-15, lines 218-220. I don't see data to support the statement on lines 218-220 that keeping synaptic interactions intact changes the network response to opioids.

Minor:

5) Please describe what the open circles on intrinsic activity figures (for example in Fig 1B) represent in the legend.

Fig 4. What do the numbers inside the boxes indicate? The legend shows stats from "correlation analysis". What test was used?

Fig. 5E. Label for green line (high gLeak) is missing.

Fig. 6D. Y axis on bursting neurons is illegible.

Add DaSilva, Cell Reports, 2023 to references listed in line 347.

Author Response

Response to Reviewers:

Reviewer 1:

The step-by-step increments on the complexity of the manipulation, building on the previous findings made it easy to follow and a good read.

Thank you very much for your time spent reviewing our manuscript. We are thrilled you found it a good read! Reviewer 2:

This manuscript takes an in silico approach to try to understand the complexities underlying variability in responses of breathing networks, particularly in the preBotC, to opioids. Models with random sparsely connected networks of 300 neurons with varying levels of intrinsic activity were generated and found to have low and high levels of sensitivity to simulated opioid. Opioid sensitivity was not predicted by intrinsic activity or the number of MOR+ neurons silenced by opioids. Higher network connectivity, especially between MOR- excitatory neurons, was protective from opioid. Surprisingly, high gLeak in MOR+ neurons made networks more susceptible to simulated opioid, even though this rendered most MOR+ neurons silent in the absence of simulated opioid. Opioid effects on network output were reversed by increasing gNaP and decreasing gLeak. This study is an appropriate use of in silico modeling to answer questions that are not accessible by other techniques. The simulated data are generally well presented and accessible, and the findings are of interest to the control of breathing field. The Discussion nicely highlights the advantages and limitations of the study, and provides physiological relevance without overstating.

Thank you very much for the careful assessment of our work and the kind comments.

I have a few questions regarding the design of the model, the parameters that were used to simulate opioid exposure and some minor comments.

1) I thought there are equal numbers of inhibitory and excitatory neurons in the preBotC. So, why does the model have only 60 inhibitory neurons compared to 240 excitatory neurons? The reviewer is correct that roughly 50% of inspiratory neurons in preBotC are inhibitory. However, the fraction of those neurons that project locally to inhibit other preBotC neurons vs. more distally to areas outside the preBotC isn't well-understood. Our model is meant to capture the "core" preBotC network and its local connections. Modeling 20% within-group inhibitory cells is consistent with our past experience modeling inhibition in the preBotC (Harris et al, 2017). In fact, in that study it was shown that 50% within-group inhibition was not able to generate a reliable rhythm due to desynchronization of rhythmic excitatory cells. In addition the to fraction of neurons that are inhibitory vs. excitatory, the relative strength and number of excitatory vs. inhibitory synapses in the preBotC is also not well-understood and could play a substantial role in determining the overall excitation/inhibition balance in the network that regulates its function.

2) If I am reading the figures correctly, the distribution of intrinsic properties (gLeak and gNaP) is not random along a continuum, but clustered into two groups. What is the basis of these clusters? You are correct, the distributions of gleak and gnap were split into high gleak and low gleak populations as done in Baertsch et al., 2021 eLife. This split distribution was done to reduce the number of intrinsic bursting neurons to be more representative of experimental observations, and also to mimic the somewhat distinct rhythm- and pattern-generating populations that are hypothesized to be present in the preBotC - the rhythm-generating population comprised of tonic and bursting neurons, and the "follower" pattern-generating populations comprised of intrinsically silent "recruitable" neurons. This also allowed us to subsequently test the effects of simulated opioids on these two preBotC populations.

3) There is no attempt to classify non-silent model neurons as inspiratory or expiratory, since MORs are expressed in both populations. Why? In Baertsch, 2021 all in silico neurons were inspiratory related. Assuming the same is true here, it is not really reflective of the physiology. Please comment.

In a previous modeling study examining the role of local and distal inhibition of the preBotC (Harris et al., 2017), we found that local inhibition naturally segmented some neurons to fire during the expiratory phase. Increasing the fraction of local "within group" inhibitory neurons increased the number of expiratory neurons, but this also made the network rhythm considerably weaker, and once ~15-20% of the network became expiratory, the rhythm fell apart. Thus a main conclusion of this study was, to achieve the number of expiratory neurons expected while maintaining a robust inspiratory rhythm, connectivity among excitatory and inhibitory neurons cannot be 100% random. Two populations with within- and between-group inhibition can contain larger numbers of inhibitory and expiratory neurons while maintaining robust rhythmogenesis. These two populations could be locally intermixed within the preBotC, anatomically distributed (e.g. BotC, PiCo, pFRG), or both.

In this study and in Baertsch et al., 2021, we used a simplified configuration of a single population with all connections drawn randomly. Because this population contains some inhibitory neurons, these models do produce a small fraction of expiratory neurons (~1-5 out of 300), but because the number of expiratory neurons is so minimal in this case, we did not specifically analyze the effects of MOR+ or MOR- expiratory neurons. Moreover, blocking inhibition, as done in Baertsch et al., 2021, transforms all expiratory neurons into inspiratory neurons, but the effects of simulated opioids remain largely unchanged, suggesting the small number of expiratory neurons in these model networks play a minimal role in their opioid response.

4) Figure 1B and 2B: There's no change in the MOR- neurons after application of the opioid condition? This seems non physiological. Wouldn't some of these MOR- neurons receive input from MOR+ synapses, or at least the condition of the network under the opioid perturbation would also alter the activity of MOR- neurons? In Baertsch, 2021, the activity of some MOR- neurons is also affected by DAMGO and simulated opioid. Is the readout of "silent" too binary in this study to detect such changes? Please comment.

Thanks for this comment. The reviewer is correct that in the synaptically coupled network, changes in the activity of MOR+ neurons would be expected to indirectly affect the activity of MOR- neurons as seen in Baertsch et al., 2021. However, the silent, bursting, tonic activities in the gnap,gleak plots shown in Figures 1B, 2B, 5B, 6C, and 7C represent the intrinsic activity of the neurons, i.e. the activity of the neurons in the absence of synaptic interactions. Because of this, the silent, bursting, or tonic activities of the MOR- population are not affected by opioid-induced changes in the intrinsic activities of the MOR+ population. We have made edits to the text in an attempt to clarify this distinction.

5) In the case where MOR+ neurons are assigned low gLeak and MOR- neurons are therefore high gLeak, nearly all nodes are silent during opioid. Where is the evidence that these nodes are "recruitable", as speculated by the authors? Wouldn't this be detected by the model and the nodes would not be classified as "silent"? Rather it seems a very small number of MOR- nodes are bursting and able to generate opioid resistant output. Please comment.

Similar to the point above, these plots depict the intrinsic activity patterns of MOR+ and MOR- neurons. The reviewer is correct that these plots show that nearly all nodes become silent in response to opioid. However, these plots show that these nodes become intrinsically silent, but this does not mean that they remain silent in the synaptically coupled network. To the contrary, what we are suggesting in the Discussion is that, although these nodes become intrinsically silent in response to opioids, they remain "recruitable" meaning that in the synaptically coupled network, they can still become rhythmic when synaptically driven by the few remaining nodes that remain intrinsically tonic/bursting (in this case, likely MOR-).

6) Page 13-15, lines 218-220. I don't see data to support the statement on lines 218-220 that keeping synaptic interactions intact changes the network response to opioids.

This may also relate to the clarifications above? Fig 5A shows network activity (with synapses intact). Without synaptic interactions, there is no network activity (rhythm) in any case under any configuration. What remains is the intrinsic activity of each neuron (shown in Fig 5B). So what we are showing in the case where opioid is targeted to the low-gleak population (middle panels) is that despite most neurons becoming intrinsically silent in response to opioid (5B), the network remains able to produce a robust rhythm when synaptic interactions are present (5A). We suggest that this is possible because many of the neurons that become intrinsically silent remain "recruitable" and become rhythmic due to synaptic inputs from the few remaining intrinsically bursting/tonic neurons in the network.

Minor:

5) Please describe what the open circles on intrinsic activity figures (for example in Fig 1B) represent in the legend.

The open circles in these plots represent each neuron within the network (all circles are open). This has been clarified in the legends.

Fig 4. What do the numbers inside the boxes indicate? The legend shows stats from "correlation analysis". What test was used? Thanks for noting this omission. The numbers in A and B represent the max and min number of each type of connection. This information is now included in the legend.

As far as we know, the name of this test is "correlation analysis", sometimes also referred to as bivariate analysis, which is used to determine whether a relationship exists between variables and then determining the magnitude and direction of that relationship. This is reflected in the Pearson correlation coefficient (r), which ranges from -1 to 1 and is represented by the red and blue shading in Fig 4.

Fig. 5E. Label for green line (high gLeak) is missing.

Thank you. This has been corrected.

Fig. 6D. Y axis on bursting neurons is illegible.

Thank you. This has been corrected.

Add DaSilva, Cell Reports, 2023 to references listed in line 347.

We have added this reference as requested.

Back to top

In this issue

Email
Print
View Full Page PDF
Citation Tools
Respond to this article
Modeling Effects of Variable preBötzinger Complex Network Topology and Cellular Properties on Opioid-Induced Respiratory Depression and Recovery
Grant M. Chou, Nicholas E. Bush, Ryan S. Phillips, Nathan A. Baertsch, Kameron Decker Harris
eNeuro 22 January 2024, 11 (3) ENEURO.0284-23.2023; DOI: 10.1523/ENEURO.0284-23.2023
Twitter logo Facebook logo Mendeley logo

Jump to section

Keywords

Responses to this article

Respond to this article

Jump to comment:

No eLetters have been published for this article.