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

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.


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 (g NaP ) and potassium dominated leak conductance (g leak ) (Butera et al., 1999;Del Negro et al., 2002;Koizumi and Smith, 2008;Yamanishi et al., 2018;Phillips and Rubin, 2019).However, g NaP , g leak , 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 g NaP or g leak 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), andBaertsch et al. (2021).First, we have the total membrane current balance equation where the currents are given by We implemented the terms D NaP (t) and D leak (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 I syn,op and I hyp,op as a mechanism for biophysical perturbations through changing I hyp,op (t) and g syn,op (t).While many terms in these differential equations are time-dependent, we only explicitly highlight the time-dependence of D NaP , D leak , g syn,op , and I hyp,op because these are set exogenously.The gating equations are with voltage-dependent steady states: , , and time constants: Finally, synapses are modeled separately for excitatory, inhibitory, and opioid-sensitive presynaptic neurons.Each synapse conductance s is modeled with first-order dynamics: , and g syn is given as the sum of g syn,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.
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 (d avg /2)/(N − 1), where d avg is the neurons' average degree (in-degree + out-degree).Our default d avg is 6, giving us a connection density of approximately 1%.However, in Figure 3, we increase the connection probability by scaling d avg , e.g., d avg = 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 g leak and g NaP .The g leak 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 g NaP 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 (D op = 1), while the inhibitory neurons and the other half of the excitatory neurons (MOR−) are insensitive (D op = 0).The opioid-sensitive population's excitatory synapses follow the I syn,op equation, whereas the insensitive neurons follow I syn,exc .Assignment of D op among excitatory neurons is random except in two cases shown in Figure 5, where opioid is applied to the excitatory neurons with g leak values below or above the median among the excitatory population.
Gradual ramp up of opioid.For opioid ramping simulations (Figs. 1, 3-5), I hyp,op ramped from 0 to 8 pA, increasing by 0.5% every 3 s, while g syn,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 I hyp,op values at the time of the last bursts with amplitudes of 10-15 Hz.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, I hyp,op = 4 pA and g syn,opioid = 0.5 nS.We varied g NaP (Fig. 6) or g leak (Fig. 7).For each node, g NaP was increased to 110%, 130%, and 150% of control values, whereas g leak 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 g NaP or g leak perturbation 9. 10 s transient period 10. 30 s simultaneous perturbation of opioid and g NaP or g leak 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 I opioid 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 g leak ∈ {0.2, 0.3, …, 1.5} nS and g NaP = {0.6,0.7, …, 1.5} nS.The points plotted within the 15 ms τ nb 10 ms τ hb 10 s g Na 28 nS g K 11.2 nS 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".

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 g NaP and g leak values.We implemented a two-population distribution of g leak and g NaP (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 g leak and g NaP 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 g leak and g NaP 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 g leak ) or the tonic/bursting neurons (low g leak ) 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 g leak population made the rhythms more resistant to opioids, whereas assigning MOR to the high g leak population made them more sensitive.When low g leak (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 g leak (primarily intrinsically silent) population is MOR+, the rhythm collapsed under only a moderate dose of opioid (I hyp,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 g leak MOR+ neurons to transition from tonic/bursting activity to silent, whereas high g leak MOR+ neurons were largely unaffected.Under non-random conditions, when all low g leak 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 g leak 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 g leak populations (Fig. 5D and E).Notably, despite a much larger proportion of the network becoming intrinsically silent, networks with low g leak 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 g leak 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 g NaP or g leak 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 g NaP 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 g leak 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 g NaP is shown in Figure 6A.Increasing g NaP by 30% in the model networks reversed the effects of opioids on burst frequency and amplitude (Fig. 6B).However, recovery of the rhythm by g NaP 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 g NaP was increased, the MOR+ neurons that were originally tonic under control conditions transitioned to bursting.Specifically, when g NaP was increased by 30%, 55% of the population entered a g NaP and g leak parameter space that supports intrinsic bursting.Thus, despite recovery of a rhythm with similar frequency and amplitude characteristics following g NaP 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 g leak (Fig. 7).The rhythmic activity of a representative network under control conditions, in opioid, and following a subsequent 10%, 30%, and 50% reduction in g leak 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 g leak was reduced to 70% of control, there was a large increase in the number of bursting neurons, and upon further reduction of g leak 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., 2005Dahan et al., , 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 (g leak ) 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 g leak (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 g leak , 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 g leak 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 g leak ), 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 g leak 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 g leak and g NaP determines intrinsic activity in our model and is sufficient to capture the silent, bursting, or tonic phenotypes of preBötC neurons.Both g leak and g NaP contribute to cellular excitability (resting membrane potential), and the voltage-dependent properties of g NaP allow some neurons with appropriate g leak 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 g leak and g NaP from set distributions were not a significant factor driving variable opioid responses, we found that scaling the distribution of g NaP or g leak across the whole population did alter the sensitivity of model networks to opioids.Interestingly, manipulation of g NaP was more effective since a 30% increase in g NaP was sufficient to restore both frequency and amplitude of the rhythm, whereas effects were more specific to burst amplitude following a 30% decrease in g leak .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.

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 g leak and g NaP 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 g leak and g NaP 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 g leak (middle) or high g leak (right) populations.B, Intrinsic activities of MOR− and MOR+ neurons (open circles) in g NaP and g leak 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 g leak 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 g leak 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).

Figure 6 .
Figure 6.Modulation of g NaP 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 g NaP to 110%, 130%, and 150% of control values.B, Quantified effects on frequency (top) and burst amplitude (bottom) during opioid and g NaP 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 g NaP and g leak 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 g NaP (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 g leak renders the network resistant to opioids.A, Example rhythm and burst waveforms from a network in response to opioid and during concurrent modulation of g leak to 90%, 70%, and 50% of control values.B, Quantified effects on frequency (top) and burst amplitude (bottom) during opioid and g leak 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 g NaP and g leak 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 g leak (n = 40, one-way RM ANOVA with Bonferroni multiple comparisons tests, ****p < 0.0001).

Table 1 .
Table of model parameters shared across neurons