Putting the theory into ‘burstlet theory’ with a biophysical model of burstlets and bursts in the respiratory preBötzinger complex
Abstract
Inspiratory breathing rhythms arise from synchronized neuronal activity in a bilaterally distributed brainstem structure known as the preBötzinger complex (preBötC). In in vitro slice preparations containing the preBötC, extracellular potassium must be elevated above physiological levels (to 7–9 mM) to observe regular rhythmic respiratory motor output in the hypoglossal nerve to which the preBötC projects. Reexamination of how extracellular K+ affects preBötC neuronal activity has revealed that low-amplitude oscillations persist at physiological levels. These oscillatory events are subthreshold from the standpoint of transmission to motor output and are dubbed burstlets. Burstlets arise from synchronized neural activity in a rhythmogenic neuronal subpopulation within the preBötC that in some instances may fail to recruit the larger network events, or bursts, required to generate motor output. The fraction of subthreshold preBötC oscillatory events (burstlet fraction) decreases sigmoidally with increasing extracellular potassium. These observations underlie the burstlet theory of respiratory rhythm generation. Experimental and computational studies have suggested that recruitment of the non-rhythmogenic component of the preBötC population requires intracellular Ca2+ dynamics and activation of a calcium-activated nonselective cationic current. In this computational study, we show how intracellular calcium dynamics driven by synaptically triggered Ca2+ influx as well as Ca2+ release/uptake by the endoplasmic reticulum in conjunction with a calcium-activated nonselective cationic current can reproduce and offer an explanation for many of the key properties associated with the burstlet theory of respiratory rhythm generation. Altogether, our modeling work provides a mechanistic basis that can unify a wide range of experimental findings on rhythm generation and motor output recruitment in the preBötC.
Editor's evaluation
This article is of significant interest to readers in the field of neural control of breathing and for researchers interested in the generation of neuronal rhythms in general. The study assembles a sophisticated computational modeling approach to test long-standing theories and emerging views in neural control of breathing and more specifically on biophysical mechanisms of burstlet generation in the respiratory network (the preBötzinger complex network). This work is an important contribution to a better understanding of the respiratory rhythm generation, will help validate (or not) running hypotheses and will guide future experiments.
https://doi.org/10.7554/eLife.75713.sa0Introduction
The complex neurological rhythms produced by central pattern generators (CPGs) underlie numerous behaviors in healthy and pathological states. These activity patterns also serve as relatively experimentally accessible instances of the broader class of rhythmic processes associated with brain function. As such, CPGs have been extensively studied using a combination of experimental and computational approaches. The inspiratory CPG located in the preBötzinger complex (preBötC) in the mammalian respiratory brainstem is perhaps one of the most intensively investigated CPGs. Despite decades of research, the mechanisms of rhythm and pattern generation within this circuit remain unresolved and highly controversial; however, it appears that the pieces may now be in place to resolve this controversy.
Much of the debate in contemporary research into the mechanisms of preBötC rhythm and pattern generation revolves around the roles of specific ion currents, such as and (Thoby-Brisson and Ramirez, 2001; Del Negro et al., 2002a; Koizumi and Smith, 2008; Koizumi et al., 2018; Picardo et al., 2019), and whether the observed rhythm is driven by an emergent network process (Rekling and Feldman, 1998; Del Negro et al., 2005; Del Negro et al., 2002b; Del Negro et al., 2002b; Rubin et al., 2009; Sun et al., 2019; Ashhad and Feldman, 2020) and/or by intrinsically rhythmic or pacemaker neurons (Johnson et al., 1994; Koshiya and Smith, 1999; Peña et al., 2004). This debate is fueled by seemingly contradictory pharmacological blocking studies (Del Negro et al., 2002a; Peña et al., 2004; Del Negro et al., 2005; Pace et al., 2007b; Koizumi and Smith, 2008) and by new experimental studies (Kam et al., 2013a; Feldman and Kam, 2015; Kallurkar et al., 2020; Sun et al., 2019; Ashhad and Feldman, 2020) that challenge existing conceptual and computational models about the generation of activity patterns in the preBötC and underlie the so-called burstlet theory of respiratory rhythm generation.
A simple but reasonable hypothesis would be that a group of dedicated preBötC neurons produces a rhythmic output that induces inspiratory movement of the diaphragm, with the strength of that output tuned by some combination of the intensity of firing of these neurons and the number of neurons that become active. The conceptual framework of burstlet theory posits a more complicated two-stage view: first, inspiratory oscillations arise from an emergent, repetitive network process in a specific preBötC subpopulation dedicated to rhythm generation. These oscillations can continue independent of their downstream impact. Second, for inspiration to occur on a particular oscillation cycle, this initial activity must recruit a secondary pattern-generating subpopulation to magnify the oscillation into a full network burst capable of eliciting motor output. This hypothesis is supported by experimental preparations that compared local preBötC neuronal activity and motor output at the hypoglossal (XII) nerve in medullary slices. These studies found that in a low excitability state (controlled by the bath K+ concentration, ), the preBötC generates a regular rhythm featuring a mixture of large and small amplitude network oscillations, dubbed bursts and burstlets, respectively, with only the bursts eliciting XII motor output (Kam et al., 2013a). Moreover, the fraction of low-amplitude preBötC events (burstlet fraction) sigmoidally decreases with increasing and only a subset of preBötC neurons are active during burstlets (Kallurkar et al., 2020). Importantly, preBötC bursts can be blocked by application of cadmium (Cd2+), a calcium channel blocker, without affecting the ongoing burstlet rhythm (Kam et al., 2013a; Sun et al., 2019), supporting the idea that rhythm generation occurs in a distinct preBötC subpopulation from pattern generation and demonstrating that conversion of a burstlet into a burst is a Ca2+-dependent process. Finally, rhythm generation in the burstlet population is hypothesized to result from an emergent network percolation process. This last idea was developed to explain holographic photostimulation experiments, which found that optically stimulating small subsets (4–9) of preBötC inspiratory neurons were sufficient to reliably evoke endogenous-like XII inspiratory bursts with delays averaging (Kam et al., 2013b). The small number of neurons required to evoke a network burst superficially seems to be at odds with reported sparse connectivities among preBötC neurons (Rekling et al., 2000), while models that can capture this effect via fast threshold modulation (Rubin and Terman, 2002) or the presentation of multiple stimulus pulses in a model of network bursting driven by synaptic dynamics (Guerrier et al., 2015) do not produce such extended delay durations. Additionally, these delays are on a similar timescale to the ramping pre-inspiratory neuronal activity that precedes network-wide inspiratory bursts, leading to the hypothesis that stimulation of this small set of preBötC neurons kicks off an endogenous neuronal percolation process underlying rhythm generation, which could be initiated by the near-coincident spontaneous spiking of a small number of preBötC neurons.
The experimental underpinning of burstlet theory challenges current ideas about inspiratory rhythm and pattern generation. However, the proposed mechanisms of burst and burstlet generation remain hypothetical and, to date, there has not been a quantitative model that provides a unified, mechanistic explanation for the key experimental observations or that validates the conceptual basis for this theory. Interestingly, key components of burstlet theory, namely, that inspiratory rhythm and pattern are separable processes and that large amplitude network-wide bursts depend on calcium-dependent mechanisms, are supported by recent experimental and computational studies. Specifically, Koizumi et al., 2018 and Picardo et al., 2019 showed that the amplitude (i.e., pattern) of preBötC and XII bursts is controlled, independently from the ongoing rhythm, by the transient receptor potential channel (TRPM4), a calcium-activated nonselective cation current (). These findings are consistent with burstlet theory as they demonstrate that rhythm and pattern are separable processes at the level of the preBötC. Moreover, these experimental observations are robustly reproduced by a recent computational modeling study (Phillips et al., 2019a), which shows that pattern generation can occur independently of rhythm generation. Consistent with burstlet theory, this model predicts that rhythm generation arises from a small subset of preBötC neurons, which in this model form a persistent sodium ()-dependent rhythmogenic kernel, and that rhythmic synaptic drive from these neurons triggers postsynaptic calcium transients, activation, and amplification of the inspiratory drive potential, which drives bursting in the rest of the network.
These recent results suggest that conversion of burstlets into bursts may be Ca2+ and dependent, occurring when synaptically triggered calcium transients in non-rhythmogenic preBötC neurons are intermittently large enough for activation to occur and to yield recruitment of these neurons into the network oscillation. The biophysical mechanism responsible for periodic amplification of Ca2+ transients is not known, however. In this computational study, we put together and build upon these previous findings to show that periodic amplification of synaptically triggered transients by calcium-induced calcium release (CICR) from intracellular stores provides a plausible mechanism that can produce the observed conversion of burstlets into bursts and can explain diverse experimental findings associated with this process. Altogether, our modeling work suggests a plausible mechanistic basis for the conceptual framework of burstlet theory and the experimental observations that this theory seeks to address.
Results
CICR periodically amplifies intracellular calcium transients
Our first aim in this work was to test whether CICR from endoplasmic reticulum (ER) stores could repetitively amplify synaptically triggered Ca2+ transients. To address this aim, we constructed a cellular model that includes the ER. The model features a Ca2+ pump, which extrudes Ca2+ from the intracellular space, a sarcoendoplasmic reticulum calcium transport ATPase (SERCA) pump, which pumps Ca2+ from the intracellular space into the ER, and the Ca2+-activated inositol trisphosphate (IP3) receptor (Figure 1A). To simulate calcium transients synaptically generated from a rhythmogenic source (i.e., burstlets), we imposed a square wave Ca2+ current into the intracellular space with varied frequency and amplitude but fixed duration (250 ms) and monitored the resulting intracellular Ca2+ transients. Depending on parameter values used, we observed various combinations of low- and high-amplitude Ca2+ responses and characterized how the fraction of Ca2+ transients that have low amplitude depends on values selected within the 2D parameter space parameterized by Ca2+ pulse frequency and amplitude. We found that the fraction of low-amplitude Ca2+ transients decreases as either or both of the Ca2+ pulse frequency and amplitude are increased (Figure 1B and example traces C1–C4).
Bursts and burstlets in a two-neuron preBötC network
Next, we tested whether the CICR mechanism (Figure 1) could drive intermittent recruitment in a reciprocally connected two-neuron network that includes one intrinsically rhythmic and one nonrhythmic neuron as a preliminary step towards considering the rhythm and pattern-generating subpopulations of the preBötC suggested by burstlet theory (Kam et al., 2013a; Cui et al., 2016; Kallurkar et al., 2020; Ashhad and Feldman, 2020) and recent computational investigation (Phillips et al., 2019a). In this network, neuron 1 is an -dependent intrinsically bursting neuron, with a burst frequency that is varied by injecting an applied current, (Figure 2A2–A3). The rhythmic bursting from neuron 1 generates periodic postsynaptic currents () in neuron 2, carried in part by Ca2+ ions, which are analogous to the square wave Ca2+ current in Figure 1. The amplitude of the postsynaptic Ca2+ transient is determined by the number of spikes per burst (Figure 2A4) and by the parameter , which determines the percentage of carried by Ca2+ ions (see ‘Materials and methods’ for a full description of these model components). Conversion of a burstlet (isolated neuron 1 burst) into a network burst (recruitment of neuron 2) is dependent on CICR (see Figure 2—figure supplement 1), which increases intracellular calcium above the threshold for activation.
In the reciprocally connected network, we first quantified the dependence of the burstlet fraction, which was defined as the number of burstlets (neuron 1 bursts without recruitment of neuron 2) divided by the total number of burstlets and network bursts (bursts in neuron 1 with recruitment of neuron 2), on and . Increasing increases the burst frequency in neuron 1 and decreases the number of spikes per neuron 1 burst (Figure 2A3 and A4), consistent with past literature (Butera et al., 1999; Del Negro et al., 2001). These changes do not strongly impact the burstlet fraction until grows enough, at which point the shorter, more rapid bursts of neuron 1 become less effective at recruiting neuron 2 and thus the burstlet fraction increases (Figure 2B2; note that increasing corresponds to a horizontal cut through the panel). In general, increasing decreased the burstlet fraction (i.e., increased the frequency of neuron 2 recruitment) by causing a larger calcium influx with each neuron 1 burst (see Figure 2B2 and C1–C4).
The burst frequency in neuron 2 is determined by the burst frequency of neuron 1 and the burstlet fraction. These effects determine the impact of changes in and on neuron 2 burst frequency (Figure 2B3). As increases, the rise in burstlet frequency implies that neuron 2 bursts in response to a smaller fraction of neuron 1 bursts, yet the rise in neuron 1 burst frequency means that these bursts occur faster. These two effects can balance to yield a relatively constant neuron 2 frequency, although the balance is imperfect and frequency does eventually increase. Increases in more straightforwardly lead to increases in neuron 2 burst frequency as the burstlet fraction drops.
Finally, the number of spikes per burst in neuron 2 is not strongly affected by changes in and (Figure 2B4), suggesting an all-or-none nature of recruitment of bursting in neuron 2. Interestingly, the period between network bursts (i.e., time between neuron 2 recruitment events) can be on the order of hundreds of seconds (e.g., Figure 2C1). This delay is consistent with some of the longer timescales shown in experiments characterizing bursts and burstlets (Kallurkar et al., 2020).
CICR supports bustlets and bursts in a data-constrained preBötC network model
Next, we tested whether the CICR mechanism presented in Figures 1 and 2 could underlie the conversion of burstlets into bursts in a larger preBötC model network including rhythm- and pattern-generating subpopulations (see ‘Data analysis and definitions’ section for details on how these are distinguished in the network setting) and whether this network could capture the -dependent changes in the burstlet fraction characterized in Kallurkar et al., 2020. sets the extracellular K+ concentration, which in turn determines the driving force for any ionic currents that flux K+. In preBötC neurons, these currents include the fast K+ current, which is involved in action potential generation, and the K+-dominated leak conductance, which primarily affects excitability (Figure 3A). In our simulations, we modeled the potassium () and leak () reversal potentials as functions of using the Nernst and Goldman–Hodgkin–Katz equations. The resulting curves were tuned to match existing data from Koizumi and Smith, 2008, as shown in Figure 3B. In our simulations, we found that intrinsic bursting is extremely sensitive to changes in . However, with increasing , intrinsic bursting could be maintained over a wide range of K+ concentrations when accompanied by increases in (Figure 3C). Additionally, the number of spikes per burst in the bursting regime increases with (Figure 3—figure supplement 1). This -dependence of is consistent with experimental data showing that neuronal input resistance decreases with increasing (Okada et al., 2005).
To construct a model preBötC network, we linked rhythm- and pattern-generating subpopulations via excitatory synaptic connections within and between the two populations (Figure 3D). We distinguished the two populations by endowing them with distinct distributions of persistent sodium current conductance (), as documented experimentally (Del Negro et al., 2002a; Koizumi and Smith, 2008). In both populations, we maintained the dependence of on (see Figure 3C and E).
For the full preBötC network model, we first characterized the impact of changes in on network behavior without calcium dynamics by setting . This network condition is analogous to in vitro preparations where all Ca2+ currents are blocked by Cd2+ and the preBötC can only generate burstlets (Kam et al., 2013a; Sun et al., 2019). Not surprisingly, with calcium dynamics blocked, we found that the network can only generate small amplitude network oscillations (burstlets) that first emerge at approximately (Figure 4A). Moreover, under these conditions, increasing results in an increase in the burstlet frequency and amplitude (Figure 4B and C), which is consistent with experimental observations (Kallurkar et al., 2020).
In the full network with calcium dynamics (), burstlets generated by the rhythmogenic subpopulation will trigger postsynaptic calcium transients in the pattern-generating subpopulation. Therefore, in this set of simulations the burstlet activity of the rhythm generating population plays an analogous role to the square wave Ca2+ current in Figure 1 and to bursts of the intrinsically rhythmic neuron in Figure 2. Hence, we characterized the burstlet fraction, burst frequency, and burst amplitude – with a burst defined as an event in which a burstlet from the rhythm generating population recruits a burst in the pattern-generating population – in the full preBötC model network as a function of and (Figure 4D–F). In this case, the frequency of the postsynaptic Ca2+ oscillation is controlled by (Figure 4B). However, because also affects burstlet amplitude (Figure 4C), the postsynaptic Ca2+ amplitude is determined by both and . If is held fixed, then modulating will only affect the amplitude of the postsynaptic Ca2+ transient since burstlet amplitude will not be impacted. The effects of selectively changing the postsynaptic Ca2+ amplitude on the full network can thus be extracted by considering a vertical slice through Figure 4E–F. Note that in the simulations that we show here burstlet generation arises from a mechanism based on ; however, we obtain similar network results if we impose burstlet activity on the burstlet subpopulation and maintain the coupling between populations and Ca2+ dynamics for burst recruitment (Figure 4—figure supplement 1).
We found that increasing or generally decreases the burstlet fraction, increases burst frequency, and slightly increases the burst amplitude (Figure 4D–F and G1–G). The decrease in the burstlet fraction with increasing or is caused by the increase in the burstlet amplitude (Figure 4C) or in Ca2+ influx with each burstlet, respectively, both of which increase the Ca2+ transient in the pattern-generating subpopulation. The increase in burst frequency with increases in or is due to the decreased burstlet fraction (i.e., the burstlet to burst transitions occur on a greater proportion of cycles) and, in the case of , by an increase in the burstlet frequency (Figure 4B). The slight increase in burst amplitude with increasing is largely due to the increase in the burstlet amplitude (Figure 4C). Figure 4H highlights the relative shape of burstlets and bursts as well as the delay between burstlet generation and recruitment of the pattern-generating population and simulated hypoglossal output, which agrees qualitatively with experimental observations (Kallurkar et al., 2020). Experimentally, it is likely that postsynaptic Ca2+ transients will increase with increasing due to the change in the resting in nonrhythmic preBötC neurons (Tryba et al., 2003) relative to the voltage-gated activation dynamics of postsynaptic calcium channels (Elsen and Ramirez, 1998); see ‘Discussion’ for a full analysis of this point. Interestingly, in our simulations, increasing (i.e., the amplitude of the postsynaptic calcium transients) with (Figure 4 traces G1–G4) generated -dependent changes in the burstlet fraction that are consistent with experimental observations (Kallurkar et al., 2020; see Figure 4I).
Note that our model includes synaptic connections from pattern-generating neurons back to rhythm-generating neurons. These connections prolong activity of rhythmic neurons in bursts, relative to burstlets, which in turn yields a longer pause before the next event (e.g., Figure 4G1). This effect can constrain event frequencies somewhat in the fully coupled network relative to the feedforward case (e.g., frequencies in Figure 4B exceed those in Figure 4E for comparable levels).
Calcium and block have distinct effects on the burstlet fraction
Next, we further characterized the calcium dependence of the burstlet to burst transition in our model by simulating calcium blockade or blockade by a progressive reduction of or , respectively. We found that complete block of synaptically triggered Ca2+ transients or block eliminates bursting without affecting the underlying burstlet rhythm (Figure 5A and B). Interestingly, progressive blockades of each of these two mechanisms have distinct effects on the burstlet fraction: blocking postsynaptic Ca2+ transients increases the burstlet fraction by increasing the number of burstlets required to trigger a network burst, whereas block only slightly increases the burstlet fraction (Figure 5C). In both cases, however, progressive blockade smoothly decreases the amplitude of network bursts (Figure 5D). The decrease in amplitude in the case of block is due to derecruitment of neurons from the pattern-forming subpopulation and a decrease in the firing rate of the neurons that remain active, whereas in the case of Ca2+ block the decrease in amplitude results primarily from derecruitment (Figure 5E and F). These simulations provide mechanism-specific predictions that can be experimentally tested.
Dose-dependent effects of opioids on the burstlet fraction
Recent experimental results by Baertsch et al., 2021 showed that opioid application locally within the preBötC decreases burst frequency but also increases the burstlet fraction. In the preBötC, opioids affect neuronal dynamics by binding to the μ-opioid receptor (μOR). The exact number of preBötC neurons expressing μOR is unclear; however, the number appears to be small, with estimates ranging from 8% to 50% (Bachmutsky et al., 2020; Baertsch et al., 2021; Kallurkar et al., 2021). Additionally, μOR is likely to be selectively expressed on neurons involved in rhythm generation, given that opioid application in the preBötC primarily impacts burst frequency rather than amplitude (Sun et al., 2019; Baertsch et al., 2021). Importantly, within the preBötC, opioids ultimately impact network dynamics through two distinct mechanisms: (1) hyperpolarization, presumably via activation of a G protein-gated inwardly rectifying potassium leak (GIRK) current (Kubo et al., 1993; Gray et al., 1999; Montandon et al., 2016), and (2) decreased excitatory synaptic transmission, presumably via decreased presynaptic release (Ballanyi et al., 2009; Wei and Ramirez, 2019; Baertsch et al., 2021).
Taking these considerations into account, we tested if our model could explain the experimental observations. Specifically, we simulated opioids as having a direct impact only on the neurons within the rhythmogenic population and their synaptic outputs (Figure 6A). To understand how preBötC network dynamics are impacted by the two mechanisms through which opioids have been shown to act, we ran separate simulations featuring either activation of GIRK channels or block of the synaptic output from the rhythmogenic subpopulation (Figure 6B–F). We found that both GIRK activation and synaptic block reduced the burst frequency (Figure 6D) and slightly increased burst amplitude (Figure 6E). The decreased frequency with synaptic block comes from an increase in the burstlet fraction, whereas GIRK activation kept the burstlet fraction constant while reducing the burstlet frequency (Figure 6F). Finally, combining these effects, we observed that simultaneously increasing the GIRK channel conductance and blocking the synaptic output of μOR-expressing neurons in our simulations generates slowing of the burst frequency and an increase in the burstlet fraction consistent with in vitro experimental data (Figure 6D–G).
Simultaneous stimulation of subsets of preBötC neurons elicits network bursts with long delays
Simultaneous stimulation of 4–9 preBötC neurons in in vitro slice preparations has been shown to be sufficient to elicit network bursts with similar patterns to those generated endogenously (Kam et al., 2013b). These elicited bursts occur with delays of several hundred milliseconds relative to the stimulation time, which is longer than would be expected from existing models. Interestingly, in the current model, due to the dynamics of CICR, there is a natural delay between the onset of burstlets and the recruitment of the follower population that underlies the transition to a burst. Therefore, we investigated whether our model could match and explain the observations seen in Kam et al., 2013b.
In our model, we first calibrated our stimulation to induce a pattern of spiking that is comparable to the patterns generated in Kam et al., 2013b (10–15 spikes with decrementing frequency, Figure 7A). We found that stimulation of 3–9 randomly selected neurons could elicit network bursts with delays on the order of hundreds of milliseconds (Figure 7B and C). Next, we characterized (1) the probability of eliciting a burst, (2) the delay in the onset of elicited bursts, and (3) the variability in delay, each as a function of the time of stimulation relative to the end of an endogenous burst (i.e., a burst that occurs without stimulation) and of the number of neurons stimulated (Figure 7D–F). In general, we found that increasing the number of stimulated neurons increases the probability of eliciting a burst and decreases the delay between stimulation and burst onset. Moreover, the probability of eliciting a burst increases and the delay decreases as the time after an endogenous burst increases (Figure 7G and H). Additionally, with its baseline parameter tuning, our model had a refractory period of approximately 1 s following an endogenous burst during which stimulation could not evoke a burst (Figure 7). The refractory period in our model is longer than measured experimentally (500 ms) (Kam et al., 2013b).
To determine the mechanisms involved in the refractoriness, we plotted the time courses of key slow variables in the model, namely, persistent sodium inactivation , ER calcium (), and synaptic depression , over one burst cycle in the absence of stimulation (see Figure 7—figure supplement 1). We found that the recovery from synaptic depression and the deinactivation of were the two slow processes with time courses that aligned with the loss of refractoriness. Thus, in our model, it appears that these two factors are crucial to the probability that a stimulus will elicit a sustained response, while calcium-related effects predominantly relate to the recruitment process by which such a response develops into a burst.
To conclude our investigation, we examined how changes in the connection probability within the pattern-forming population () affect the refractory period, probability, and delay of evoked bursts following simultaneous stimulation of 3–9 randomly selected neurons in the preBötC population. We focused on the pattern-forming population because it comprises 75% of the preBötC population, and, therefore, neurons from this population are most likely to be stimulated and the synaptic projections from these neurons are most likely to impact the properties of evoked bursts. To avoid a confound that would arise if increased connection probability led to overall stronger synaptic input, we adjusted to compensate for changes in and keep the network synaptic strength, defined as , at a constant value.
With this scaling, we found that decreasing/increasing decreased/increased the refractory period (Figure 8A–C) by impacting the probability of eliciting a burst in the period immediately after an endogenous burst (Figure 8D and E). More specifically, the change in the probability of evoking a burst, with decreased/increased , is indicated by a leftward/rightward shift in the probability vs. stimulation time curves relative to a control level of () (see Figure 8D and E). That is, relatively small connection probabilities with large connection strengths lead to network dynamics with a shorter refractory period when stimulation cannot elicit a burst and a higher probability that a stimulation at a fixed time since the last burst will evoke a new burst.
It may seem surprising that networks with smaller connection probabilities exhibit a faster emergence of bursting despite their larger connection weights since intuitively, with lower connection probabilities, fewer neurons could be recruited by each action potential, resulting in longer, more time-consuming activation pathways. A key point, however, is that when connection weights are larger, fewer temporally overlapping inputs are needed to recruit each inactive neuron. For example, suppose that we fix and , and we take to scale as . The minimal number of inputs from active neurons needed to activate an inactive neuron depends on the synaptic weight, . Let denote this number for the specific value of that we have selected. We can approximate the expected number of neurons receiving or more inputs from active neurons by computing the expected number receiving exactly inputs, which we denote as , where the brackets indicate an expectation or average. For a network with a random connectivity profile, this expected value is computed from the binomial formula as
Suppose that next we consider another network in which we double and halve , thus keeping their product constant. For this smaller , more inputs will be needed to activate an inactive neuron. Specifically, assume that now at least inputs are needed for activation. The expected number of neurons receiving inputs, , is given by
An elementary calculation shows that for relevant parameter values (such as and small as indicated by the stimulation experiments). Thus, increasing and proportionally scaling down reduces the chance of successful recruitment of inactive neurons by active neurons.
Interestingly, our simulations suggest that the connection probability in the pattern-generating population must be between 1% and 2% to match the approximately 500 ms refractory period measured experimentally (Kam et al., 2013b; Figure 8F). Surprisingly, the mean and distribution of delays from stimulation to burst for all successfully elicited bursts are not strongly affected by changes in (Figure 8F). For a given stimulation time and number of neurons stimulated, however, decreasing decreases the delay of elicited bursts (Figure 8G). Finally, because the neurons in the pattern-generating population appear to play a dominant role in determining if stimulation will elicit a network burst, we characterized how the number of pattern-generating neurons stimulated, out of a total set of nine stimulated neurons, affects the probability of eliciting a network burst as a function of stimulation time (Figure 8H). These simulations were carried out under a baseline condition of . In general, we found that stimulating a relatively larger proportion of pattern-generating neurons increased the probability of eliciting a network burst for all times after the approximately 1 s refractory period, as indicated by the positive slope in Figure 8H. Additionally, eliciting a network burst does not require stimulation of rhythmogenic neurons.
Discussion
Recent experiments have revealed a decoupling of respiratory rhythm generation and output patterning in the preBötC, which has given rise to the conceptual framework of burstlet theory. To date, however, this theory lacks the quantitative basis, grounded in underlying biophysical mechanisms, needed for its objective evaluation. To address this critical gap, in this computational study we developed a data-constrained biophysical model of the preBötC that generates burstlets and bursts as proposed by burstlet theory, with a range of features that match experimental observations. To summarize, we first show that CICR from intracellular stores is a natural mechanism to periodically amplify postsynaptic calcium transients needed for activation and recruitment of pattern-forming neurons into network bursts (Figure 1). Next, we demonstrate that in a two-neuron network CICR can convert baseline rhythmic activity into a mixture of bursts and burstlets, where the burstlet fraction depends largely on the magnitude of postsynaptic calcium transients (Figure 2). In a larger preBötC network containing rhythm- and pattern-forming subpopulations with experimentally constrained intrinsic properties, population sizes, and synaptic connectivity probabilities (Figure 3), similar but more realistic activity patterns arise (Figure 4). Moreover, we show that this model can match a wide range of the key experimental underpinnings of burstlet theory: dependence of the burstlet fraction on extracellular potassium concentration (Figure 4I), the Ca2+ dependence of the burstlet-to-burst transition (Figure 5), the effects of opioids on burst frequency and burstlet fraction (Figure 6), and the long delay and refractory period of bursts evoked by holographic photostimulation of small subsets of preBötC neurons (Figures 7 and 8).
Insights into the mechanisms of burst (pattern) and burstlet (rhythm) generation within the inspiratory preBötC
Burstlet theory to date has largely been an empirical description of the observed features of bursts and burstlets. One idea that has been suggested is that rhythm generation is driven by a stochastic percolation process in which tonic spiking across the rhythm-generating population gradually synchronizes during the inter-burst-interval to generate the burstlet rhythm. Subsequently, a burst (i.e., motor output) only occurs if the burstlet is of sufficient magnitude, resulting from sufficient synchrony, to trigger all-or-none recruitment of the pattern-forming subpopulation (Kam et al., 2013a; Kam et al., 2013b; Feldman and Kam, 2015; Cui et al., 2016; Kallurkar et al., 2020; Ashhad and Feldman, 2020). This theory, however, does not identify or propose specific biophysical mechanisms capable of generating a quantitative explanation of the underlying cellular and network-level dynamics, fails to capture the Ca2+ dependence of the burst-to-burstlet transition, and cannot explain how extracellular potassium concentration impacts the burstlet fraction. Our simulations support an alternative view of the recruitment process associated with this transition that builds directly from previous computational studies (Jasinski et al., 2013; Phillips et al., 2019a; Phillips and Rubin, 2019b; Phillips et al., 2021), which robustly reproduce a wide array of experimental observations. Specifically, in this study we show that amplification of postsynaptic calcium transients in the pattern-generating subpopulation (triggered by burstlets) provides a natural mechanism capable of explaining the Ca2+ dependence of the burstlet-to-burst transition.
Importantly, our model yields the result, and hence the prediction, that the burstlet fraction is determined by the probability that a burstlet will trigger CICR in the pattern-forming subpopulation. In the model, this probability is determined by the magnitude of postsynaptic calcium transients as well as the activation dynamics of the IP3 receptor and the SERCA pump. Therefore, to explain the decrease in the burstlet fraction with increasing extracellular , the magnitude of the burstlet-triggered postsynaptic calcium transients must increase with . Some of this rise can result directly from the increase in burstlet amplitude with increasing (see Kallurkar et al., 2020 and Figure 4C). To fully match the experimentally observed relationship between and the burstlet fraction (Figure 4J), we also explicitly increased the parameter , which sets the proportion of the postsynaptic calcium current carried by Ca2+. Thus, our model predicts that the magnitude of postsynaptic Ca2+ transients triggered by EPSPs should increase as is elevated.
This same prediction arises from considering the voltage-dependent properties of Ca2+ channels characterized in preBötC neurons and the changes in the membrane potential of non-rhythmogenic (i.e., pattern-forming) neurons as a function of . Specifically, it is likely that voltage-gated calcium channels are involved in generating the postsynaptic Ca2+ transients as dendritic Ca2+ transients have been shown to precede inspiratory bursts and to be sensitive to Cd2+, a calcium channel blocker (Del Negro et al., 2011). Consistent with this idea, Cd2+-sensitive Ca2+ channels in preBötC neurons appear to be primarily localized in distal dendritic compartments (Phillips et al., 2018). Voltage-gated calcium channels in the preBötC start to activate at approximately (Elsen and Ramirez, 1998), and importantly, the mean somatic resting membrane potential of non-rhythmogenic preBötC neurons increases from to when extracellular potassium concentration is elevated from to (Tryba et al., 2003). Moreover, at , EPSPs in the preBötC are on the order of 2–5 mV (Kottick and Del Negro, 2015; Morgado-Valle et al., 2015; Baertsch et al., 2021) and the amplitude of EPSCs has been shown to decrease as is lowered (Okada et al., 2005). Putting together these data on resting membrane potential, EPSP sizes, and voltage-dependent activation of Ca2+ channels, we deduce that when , the magnitude of EPSPs may not reach voltages sufficient for significant activation of voltage-gated Ca2+ channels. As is increased, however, both EPSC magnitudes and the membrane potential of pattern-forming neurons increase. Therefore, with increased , the prediction is that EPSCs will result in greater activation of voltage-gated Ca2+ channels and increased postsynaptic calcium transients. This effect is captured in the model by an increase in the parameter , which determines the percentage of the postsynaptic current carried by Ca2+ ions, with .
The idea that dendritic postsynaptic Ca2+ transients and activation play a critical role in regulating the pattern of preBötC dynamics is well supported by experimental and computational studies. Specifically, the dendritic Ca2+ transients that precede inspiratory bursts (Del Negro et al., 2011) have been shown to travel in a wave to the soma, where they activate TRPM4 currents () (Mironov, 2008). Moreover, the rhythmic depolarization of otherwise non-rhythmogenic neurons (inspiratory drive potential) depends on (Pace et al., 2007a), while the inspiratory drive potential is not dependent on Ca2+ transients driven by voltage-gated calcium channels expressed in the soma (Morgado-Valle et al., 2008). Finally, pharmacological blockade of TRPM4 channels, thought to represent the molecular correlates of , reduces the amplitude of preBötC motor output without impacting the rhythm (Koizumi et al., 2018; Picardo et al., 2019). These experimental findings were incorporated into and robustly reproduced in a recent computational model (Phillips et al., 2019a) and are reproduced here (see Figure 5B and D). Consistent with these findings, this previous model suggests that rhythm generation arises from a small subset of preBötC neurons, which form an -dependent rhythmogenic kernel (i.e., burstlet rhythm generator), and that rhythmic synaptic drive from these neurons triggers postsynaptic calcium transients, activation, and amplification of the inspiratory drive potential, which spurs bursting in the rest of the network. This study builds on this previous model by explicitly defining rhythm- and pattern-generating neuronal subpopulations (see Figure 3) and incorporating the mechanisms required for CICR and intermittent amplification of postsynaptic calcium transients.
CICR mediated by the SERCA pump and the IP3 receptor has long been suspected to be involved in the dynamics of preBötC rhythm and/or pattern generation (Pace et al., 2007a; Crowder et al., 2007; Mironov, 2008; Toporikova et al., 2015) and has been explored in individual neurons and network models of the preBötC (Toporikova and Butera, 2011; Jasinski et al., 2013; Rubin et al., 2009; Wang and Rubin, 2020). Experimental studies have not clearly established the role of CICR from ER stores in respiratory circuits, however. For example, Mironov, 2008 showed that application of 1 µM thapsigargin, a SERCA pump inhibitor, abolished rhythmic activity and blocked calcium transients that travel in a wave from the dendrites to the soma. In a separate study, however, block of the SERCA pump by bath application of thapsigargin (2–20 µM) or cyclopiazonic acid (CPA) (30–50 µM) did not significantly affect the amplitude or frequency of hypoglossal motor output in in vitro slice preparations containing the preBötC (Beltran-Parrazal et al., 2012). The explanation for these seemingly contradictory experimental results is unclear, especially since effects of SERCA pump block could be complicated, and will need to be investigated by future studies. It is possible that the role of CICR may be dynamically regulated depending on the state of the preBötC network. For example, the calcium concentration at which the IP3 receptor is activated is dynamically regulated by IP3 (Kaftan et al., 1997), and therefore, activity- or neuromodulatory-dependent changes in the cytoplasmic Ca2+ and/or IP3 concentration may impact ER Ca2+ uptake and release dynamics. Store-operated Ca2+ dynamics are also affected by the transient receptor potential canonical 3 (TRPC3) channels (Salido et al., 2009), which are expressed in the preBötC, and manipulation of TRPC3 has been shown to impact burst amplitude and regularity (Tryba et al., 2003; Koizumi et al., 2018) as would be predicted by this model. It is also possible that calcium release is mediated by the ryanodine receptor, an additional calcium-activated channel located in the ER membrane (Lanner et al., 2010), since bath application of CPA (100 µM) and ryanodine (10 µM) removed large-amplitude oscillations in recordings of preBötC population activity (Toporikova et al., 2015).
Finally, we note that while various markers can be used to define distinct subpopulations of neurons within the preBötC, our model cannot determine which of these ensembles are responsible for rhythm and pattern generation. Past experiments have examined the impact of optogenetic inhibition, applied at various intensities to subpopulations associated with specific markers, on the frequency of inspiratory neural activity, but this activity was measured by motor output, not within the preBötC itself (Tan et al., 2008; Cui et al., 2016; Koizumi et al., 2016). According to burstlet theory and our model, slowed output rhythmicity could derive from inhibition of rhythm-generating neurons, due to a reduced frequency of burstlets, and from inhibition of pattern-generating neurons, due to a reduced success rate of burst recruitment. Thus, measurements within the preBötC will be needed in order to assess the mapping between subpopulations of preBötC neurons and roles in burstlet and burst production.
Additional comparisons to experimental results
In our model (Figure 4), a burstlet rhythm first emerges at a of approximately 5 mM, whereas in the experiments of Kallurkar et al., 2020, the burstlet rhythm continues even down to 3 mM. To explain this discrepancy, we note that our model assumes that the extracellular potassium concentration throughout the network is equal to . Respiratory circuits appear to have some buffering capacity, however, such that for concentrations below approximately 5 mM the extracellular K+ concentration remains elevated above (Okada et al., 2005). The range over which our model generates a rhythm would extend to that seen experimentally if extracellular K+ buffering were accounted for. This buffering effect can also explain why the burstlet fraction remains constant in experimental studies when is lowered from 5 mM to 3 mM (Kallurkar et al., 2020). Our model also does not incorporate short-term extracellular potassium dynamics that depend on and may impact the ramping shape of burstlet onset (Abdulla et al., 2021). Importantly, over the range of values relevant both to experiments and our model, we find clear agreement on the dependence of burstlet fraction on (Figure 4I).
Although our model incorporates various key biological features, it does not include some of the biophysical mechanisms that are known to shape preBötC patterned output or that are hypothesized to contribute to the properties described by burstlet theory. For example, the M-current associated with KCNQ potassium channels has been shown to impact burst duration by contributing to burst termination (Revill et al., 2021). Additionally, intrinsic conductances associated with a hyperpolarization-activated mixed cation current () and a transient potassium current () are hypothesized to be selectively expressed in the pattern- and rhythm-generating preBötC subpopulations (Picardo et al., 2013; Phillips et al., 2018). Thus, our model predicts that while these currents may impact quantitative properties of burstlets and bursts, they are not critical for the presence of burstlets and their transformation into bursts. The current model also does not include a population of inhibitory preBötC neurons. Inhibition is involved in regulating burst amplitude (Baertsch et al., 2018), but it does not have a clear role in burst or burstlet generation, and therefore inhibition was omitted from this work. More globally, it is crucial to recognize that areas outside of the preBötC impact dynamics within the preBötC. These effects, which remain to be fully elucidated, may range from ongoing modulation of the level of excitability of preBötC neurons to timed signaling that contributes to preBötC rhythmicity and patterning (e.g., Mulkey et al., 2004; Dutschmann and Dick, 2012; Phillips et al., 2012; Smith et al., 2013; Dhingra et al., 2019; Richter et al., 2019; Liu et al., 2022). For example, transection studies suggest that pontine regions may make crucial contributions to respiratory circuit excitability and respiratory pattern formation (Jones and Dutschmann, 2016; Smith et al., 2007). Finally, the data on which this study was based comes from a variety of settings, including in vitro and other reduced preparations, and additional factors no doubt complicate the generation and control of respiratory outputs in vivo. Indeed, although experimental results suggest that manipulations to enhance preBötC excitability in slice preparations do not appear to significantly impact the mechanisms of preBötC rhythmicity or the generation of bursts and burstlets, additional investigation of how higher brainstem centers impact preBötC inspiratory rhythm and pattern generation is an important direction for future studies.
Importantly, our model does robustly reproduce all of the range of key experimental observations underlying burstlet theory. Not surprisingly, block of calcium transients or in our model eliminates bursts without affecting the underlying rhythm (Figure 5), which is consistent with experimental observations (Kam et al., 2013b; Sun et al., 2019). Interestingly, our model also provides the experimentally testable predictions that blocking calcium transients will increase the burstlet fraction while block will have no effect on this fraction, whereas both perturbations will smoothly reduce burst amplitude. The calcium-dependent mechanisms that we include in our model pattern-generating population have some common features with a previous model that suggested the possible existence of two distinct preBötC neuronal populations responsible for eupneic burst and sigh generation, respectively, which also included excitatory synaptic transmission from the former to the latter (Toporikova et al., 2015). In the eupnea-sigh model, however, the population responsible for low-frequency, high-amplitude sighs was capable of rhythmic burst generation even without synaptic drive, in contrast to the pattern-generation population as tuned in our model. Also, in contrast to the results on bursts considered in our study, sigh frequency in the earlier model did not vary with extracellular potassium concentration and sigh generation required a hyperpolarization-activated inward current, .
We also considered the effects of opioids in the context of burstlets and bursts, a topic that has not been extensively studied. It is well established that opioids slow the preBötC rhythm in in vitro slice preparations; however, the limited results presented to date on effects of opioids on the burstlet fraction are inconsistent. Specifically, Sun et al., 2019 found that application of the μOR agonist DAMGO at 10 nM and 30 nM progressively decreased the preBötC network frequency but had no impact on the burstlet fraction before the network rhythm was eventually abolished at approximately 100 nM. Similarly, Baertsch et al., 2021 found that DAMGO decreased the preBötC network frequency in a dose-dependent fashion; however, in these experiments the network was less sensitive to DAMGO, maintaining rhythmicity up to approximately 300 nM, and the burstlet fraction increased with increasing DAMGO concentration. The inconsistent effects of DAMGO on the burstlet fraction across these two studies can be explained by our simulations based on the different sensitivities of these two preparations to DAMGO and the two distinct mechanisms of action of DAMGO on neurons that express μOR – decreases in excitability and decreases in synaptic output of neurons – identified by Baertsch et al., 2021. In our simulations, we show that the decreased excitability resulting from activation of a GIRK channel only impacts frequency, whereas decreasing the synaptic output of μOR-expressing neurons results in an increase in the burstlet fraction and a decrease in burst frequency (Figure 6). In experiments, suppression of synaptic output does not appear to occur until DAMGO concentrations are above approximately 50 nM(Baertsch et al., 2021). Therefore, it is not surprising that DAMGO application did not strongly impact the burstlet fraction before the rhythm was ultimately abolished in Sun et al., 2019 due to the higher DAMGO sensitivity of that particular experimental preparation, as indicated by the lower dose needed for rhythm cessation.
Mixed-mode oscillations
Mixed-mode oscillations, in which intrinsic dynamics of a nonlinear system naturally lead to alternations between small- and large-amplitude oscillations (Del Negro et al., 2002c; Bertram and Rubin, 2017), are a mechanism that has been previously proposed to underlie bursts and burstlets under the assumption of differences in intrinsic oscillation frequencies across preBötC neurons (Bacak et al., 2016). This mechanism was not needed to explain the generation of bursts and burstlets in the current model, however. Moreover, systems with mixed-mode oscillations can show a wide range of oscillation amplitudes under small changes in conditions and mixed-mode oscillations only emerge in the preBötC when is elevated above 9 mM (Del Negro et al., 2002c). These properties are not consistent with the burst and burstlet amplitudes or -dependent changes in the burstlet fraction seen experimentally (Kallurkar et al., 2020) and in our model.
Holographic photostimulation, percolation, and rhythm generation
Experimental data supporting burstlet theory has shown that burstlets are the rhythmogenic event in the preBötC. However, although burstlet theory is sometimes referenced as a theory of respiratory rhythm generation, the actual mechanisms of burstlet rhythm generation remain unclear. One idea that has been suggested is that rhythm generation is driven by a stochastic percolation process in which tonic spiking across the rhythm-generating population gradually synchronizes during the inter-burst-interval to generate the burstlet rhythm (Ashhad and Feldman, 2020; Slepukhin et al., 2020). In this framework, a burst (i.e., motor output) only occurs if the burstlet is of sufficient magnitude, resulting from sufficient synchrony, to trigger all-or-none recruitment of the pattern-forming subpopulation (Kam et al., 2013a; Kam et al., 2013b; Feldman and Kam, 2015; Kallurkar et al., 2020; Ashhad and Feldman, 2020; Slepukhin et al., 2020).
The idea that burstlets are the rhythmogenic event within the preBötC is supported by the observation that block of voltage-gated Ca2+ channels by Cd2+ eliminates bursts without affecting the underlying burstlet rhythm (Kam et al., 2013a; Sun et al., 2019). However, the rhythmogenic mechanism based on percolation is speculative and comes from two experimental observations. The first is that the duration and slope (i.e., shape) of the burstlet onset are statistically indistinguishable from the ramping pre-inspiratory activity that immediately precedes inspiratory bursts (Kallurkar et al., 2020). We note, however, that this shape of pre-inspiratory activity can arise through intrinsic mechanisms at the individual neuron level (Abdulla et al., 2021). The second observation evoked in support of the percolation idea is that holographic photostimulation of small subsets (4–9) of preBötC neurons can elicit bursts with delays lasting hundreds of milliseconds (Kam et al., 2013b). These delays are longer than could be explained with existing preBötC models and have approximately the same duration as the pre-inspiratory activity and burstlet onset hypothesized to underlie the rhythm. According to the percolation hypothesis of burstlet rhythm generation, these long delays result from the specific topological architecture of the preBötC, recently proposed to be a heavy-tailed synaptic weight distribution in the rhythmogenic preBötC subpopulation (Slepukhin et al., 2020).
Interestingly, the model presented here naturally captures the long delays characterized by Kam et al., 2013b, and stimulation of small subsets of neurons triggers a growth in population activity in the lead up to a burst that could be described as percolation (Figure 7B). Our model does not require a special synaptic weight distribution to generate the long delays, however. Indeed, our model suggests that the long delays between simulation and burst generation are due in large part to the dynamics of the pattern-forming population, as probabilistically these neurons are the most likely to be stimulated and they appear to play a dominant role in setting the timing of the elicited burst response (Figure 8H). Moreover, the dynamics of this population is strongly impacted by the CICR mechanism proposed here, which is required for burst generation. Interestingly, to match the 500 ms refractory period following an endogenous burst during which holographic stimulation cannot elicit a burst, our model predicts that the connection probability in the pattern-generating preBötC subpopulation must be between 1% and 2% (Figure 8A and B), which is consistent with available experimental data (Ashhad and Feldman, 2020). Experiments applying global, presumably weaker stimulation to the preBötC yield longer (~2 s) refractory periods after endogenous bursts (Baertsch et al., 2018; Kottick and Del Negro, 2015), and our model can also produce similar refractory periods in analogous conditions.
Thus, taken together, previous modeling and our work offer two alternative, seemingly viable hypotheses about the source of the delay between holographic stimulation and burst onset, each related to a proposed mechanism for burstlet and burst generation. Yet additional arguments call into question aspects of the percolation idea. If the burstlet rhythm is driven by a stochastic percolation process, then the period and amplitude of burstlets should be stochastic, irregular, and broadly distributed. Moreover, in the proposed framework of burstlet theory, the pattern of bursts and burstlets for a given burstlet fraction would also be stochastic since the burstlet-to-burst transition is thought to be an all-or-none process that depends on the generation of a burstlet of sufficient magnitude. Example traces illustrating a mixture of bursts and burstlets typically show a pattern of multiple burstlets followed by a burst that appears to regularly repeat (Kam et al., 2013b; Sun et al., 2019; Kallurkar et al., 2020) and hypoglossal output timing has also been found to exhibit high regularity Kam et al., 2013b, however, suggesting that the burstlet-to-burst transition is not dependent on the synchrony and hence magnitudes of individual burstlets but rather on a slow process that gradually evolves over multiple burstlets. The regularity and patterns of burstlets and bursts that arise from such a process in our model match well with those observed experimentally.
We note that the burstlet-to-burst transition mechanism proposed here, based on CICR from ER stores, depends on rhythmic inputs from the rhythm-generating to the pattern-generation population; however, it is independent of the mechanism of rhythm generation. In our simulations, rhythm generation depends on the slowly inactivating persistent sodium current (). The role of in preBötC inspiratory rhythm generation is a contentious topic within the field, largely due to the inconsistent effects of block. We chose to use as the rhythmogenic mechanism in the burstlet population for a number of reasons: (1) consideration of the pharmacological mechanism of action and nonuniform effects of drug penetration can explain the seemingly contradictory experimental findings relating to (Phillips and Rubin, 2019b), (2) -dependent rhythm generation is a well-established and understood idea (Butera et al., 1999), (3) recent computational work on which the current model is based suggests that rhythm generation occurs in a small, -dependent rhythmogenic kernel that is analogous to the burstlet population (Phillips et al., 2019a), and predictions from this model that depend on the specific proposed roles of and in rhythm and pattern formation have been experimentally confirmed in a recent study (Phillips et al., 2021). It is important to note, however, that the findings about burstlets and bursts presented in this work would have been obtained if the burstlet rhythm was imposed (Figure 4—figure supplement 1) or if burstlets were generated by some other means, such as by the percolation mechanism proposed by burstlet theory.
Summary of model predictions
The model presented here is itself a prediction; that is, this work predicts that a CICR-mediated mechanism is critical to the transition of burstlets into bursts. At a more specific level, our model makes the following predictions: (1) the magnitude of postsynaptic calcium transients triggered by EPSCs in preBötC neurons will increase with K+ (see Figure 4 and related text); (2) network-level burstlets and bursts will persist if currents involved in regulating burst shape, such as and , are blocked (see earlier discussion); (3) blocking postsynaptic Ca2+ transients will increase the burstlet fraction and decrease the burst amplitude before network bursts are eventually abolished (see Figure 5); (4) block will not change the burstlet fraction and will decrease burstlet amplitudes (see Figure 5); (5) the synaptic connection probability within the pattern-generating population in the preBötC is low (1–2%, see Figure 8); and (6) selective holographic stimulation of pattern-generating neurons should be more effective than stimulation of rhythm-generating neurons at triggering network bursts (see Figure 8). This could be tested by selectively stimulating Dbx1 preBötC neurons that express Sst (pattern forming) or that do not express Sst (rhythmogenic).
Conclusions
This study has developed the first model-based description of the biophysical mechanism underlying the generation of bursts and burstlets in the inspiratory preBötC. As suggested by burstlet theory and other previous studies, rhythm and pattern generation in this work are represented by two distinct preBötC subpopulations. A key feature of our model is that generation of network bursts (i.e., motor output) requires amplification of postsynaptic Ca2+ transients by CICR in order to activate and drive bursting in the rest of the network. Moreover, the burstlet fraction depends on rate of Ca2+ buildup in intracellular stores, which is impacted by -dependent modulation of preBötC excitability. These ideas complement other recent findings on preBötC rhythm generation (Phillips et al., 2019a; Phillips and Rubin, 2019b; Phillips et al., 2021), together offering a unified explanation for a large body of experimental findings on preBötC inspiratory activity that form a theoretical foundation on which future developments can build.
Materials and methods
Neuron model
Request a detailed protocolModel preBötC neurons include a single compartment and incorporate Hodgkin–Huxley-style conductances adapted from previously described models (Jasinski et al., 2013; Phillips et al., 2019a; Phillips and Rubin, 2019b) and/or experimental data as detailed below. The membrane potential of each neuron is governed by the following differential equation:
where is the membrane capacitance and each represents a current, with i denoting the current’s type. The currents include the action potential generating Na+ and delayed rectifying K+ currents ( and ), persistent Na+ current (), voltage-gated Ca2+ current (), Ca2+-activated nonselective cation (CAN) current (), K+-dominated leak current (), synaptic current (), μ-opioid receptor-activated G protein-coupled inwardly rectifying K+ leak current () (Kubo et al., 1993), and a holographic photostimulation current (). denotes an applied current injected from an electrode. The currents are defined as follows:
where gi is the maximum conductance, is the reversal potential, and mi and hi are gating variables for channel activation and inactivation for each current . The glutamatergic synaptic conductance is dynamic and is defined below. The values used for the gi and are mostly shown in Table 1 , with a few conductances selected from distributions as indicated in Table 2.
Activation (mi) and inactivation (hi) of voltage-dependent channels are described by the following differential equation:
where represents steady-state activation/inactivation and is a time constant. For , , and , the functions and take the forms
The values of the parameters (, , , , and ) corresponding to , and are given in Table 1.
For the voltage-gated potassium channel, the steady-state activation and time constant are given by the expressions
where
The values for the constants , , , , , and are also given in Table 1.
activation depends on the concentration in the cytoplasm () and is given by
The parameters and represent the half-activation Ca2+ concentration and the Hill coefficient, respectively, and are included in Table 1.
The dynamics of is determined in part by the balance of Ca2+ efflux toward a baseline concentration via the Ca2+ pump and Ca2+ influx through voltage-dependent activation of and synaptically triggered Ca2+ transients, with a percentage () of the synaptic current () carried by Ca2+ ions. Additionally, the model includes an intracellular compartment that represents the ER, which impacts . The ER removes Ca2+ from the cytoplasm via a sarcoplasmic/endoplasmic reticulum Ca2+- ATPase (SERCA) pump, which transports Ca2+ from the cytoplasm into the ER () and releases Ca2+ into the cytoplasm via calcium-dependent activation of the inositol triphosphate (IP3) receptor (). Therefore, the dynamics of is described by the following differential equation:
where is a conversion factor relating current to rate of change of , is the time constant for the Ca2+ pump, is a minimal baseline calcium concentration, and is the ratio of free to bound Ca2+ in the ER.
The flux of Ca2+ from the ER to the cytoplasm through the IP3 receptor is modeled as
where represents the leak constant from the ER stores, represents the permeability of the IP3 channel, and are dissociation constants, and is the cytoplasm IP3 concentration. Finally, the Ca2+-dependent IP3 gating variable, , and the Ca2+ concentration in the ER, , are determined by the following equations:
where is a conversion factor, is the dissociation constant for IP3 inactivation, is the total intracellular calcium concentration, and is the ratio of cytosolic to ER volume. The total intracellular calcium concentration is described as
Finally, removal of Ca2+ from the cytoplasm by the SERCA pump is modeled as
where is the maximal flux through the SERCA pump, and is a dissociation constant.
Nondimensionalization of similar models in past work (Wang and Rubin, 2017; Wang and Rubin, 2020) has shown that , and are the slowest variables in the model and evolve on similar timescales, while evolves on a faster timescale that is still significantly slower than that of the voltage dynamics and other current gating variables. Some subtleties arise in that different components of the calcium dynamics evolve on different timescales and their influences depend on the levels of calcium present in various domains within the cell, but these subtleties are not considered in this article.
When we include multiple neurons in the network, we can index them with subscripts. The total synaptic conductance of the ith target neuron is described by the following equation:
where is a fixed or tonic excitatory synaptic conductance (e.g., from respiratory control areas outside of the preBötC) that we assume impinges on all neurons, represents the weight of the synaptic connection from neuron to neuron i, is a scaling factor for short-term synaptic depression in the presynaptic neuron (described in more detail below), is an element of the connectivity matrix ( if neuron makes a synapse with neuron i and otherwise), is the Heaviside step function, and denotes time. is an exponential synaptic decay constant, while is the time at which the nth action potential generated by neuron reaches neuron i.
We included synaptic depression in our model because experiments have revealed that it contributes to termination of inspiratory activity in the preBötC (Kottick and Del Negro, 2015) and past computational models have suggested that it might play an important role in preBötC network oscillations (Rubin et al., 2009; Guerrier et al., 2015). Synaptic depression in the jth neuron () was simulated using an established mean-field model of short-term synaptic dynamics (Abbott et al., 1997; Dayan and Abbott, 2001; Morrison et al., 2008) as follows:
where the parameter sets the maximum value of , sets the rate of recovery from synaptic depression, sets the fractional depression of the synapse each time neuron spikes, and is the Kronecker delta function that equals 1 at the time of each spike in neuron and 0 otherwise. Parameters were chosen to qualitatively match data from Kottick and Del Negro, 2015. Note that with this choice of synaptic depression recovers on a timescale comparable to that of the other slowest variables in the model.
When we consider a two-neuron network (Figure 2), we take and . For the full preBötC population model comprising rhythm- and pattern-generating subpopulations, the weights of excitatory conductances were uniformly distributed such that where is a constant associated with the source and target neurons’ populations; with each such pair, we also associated a connection probability and used this to randomly set the values (see Table 3). Effects of opioids on synaptic transmission for source neurons in the rhythmogenic subpopulation (Figure 6) were simulated by scaling with the parameter , which ranged between 0 and 0.5 and sets the percent synaptic block.
Network construction
Request a detailed protocolThe relative proportions of neurons assigned to the rhythm- and pattern-generating preBötC subpopulations were chosen based on experimental data. For example, Kallurkar et al., 2020 found that of preBötC inspiratory neurons are active during burstlets at . Moreover, the rhythm- and pattern-generating neurons are hypothesized to be represented by the subsets of Dbx1-positive preBötC neurons that are somatostatin-negative () and -positive (), respectively (Cui et al., 2016; Ashhad and Feldman, 2020). Somatostatin-positive neurons are estimated to comprise 72.6% of the preBötC population (Koizumi et al., 2016). Therefore, our preBötC network was constructed such that the rhythm and pattern-forming subpopulations represent 25% and 75% of the neuron preBötC population (i.e., and ). The rhythm- and pattern-generating neurons are distinguished by their , , and conductances. Also, we included the K+ leak current exclusively to the rhythm generating subpopulation, the activation of which we used as one representation of the effects of opioid application (Figure 6).
The synaptic connection probabilities within the rhythm- and pattern-generating neurons, and , were taken from previous experimental findings (Rekling et al., 2000 and Ashhad and Feldman, 2020, respectively). The connection probabilities between the rhythm- and pattern-generating populations are not known and in the model were set at such that the total connection probability in the network is approximately 13% (Rekling et al., 2000).
Heterogeneity was introduced by normally distributing the parameters , , and as well as uniformly distributing the weights () of excitatory synaptic connections (see Table 2 and Table 3). Additionally, was conditionally distributed with in order to achieve a bivariate normal distribution between these two conductances, as suggested by Del Negro et al., 2002a and Koizumi and Smith, 2008. In our simulations, this was achieved by first normally distributing in each neuron according to the values presented in Table 2. Then we used a property of bivariate normal distribution, which says that the conditional distribution of given is itself a normal distribution with mean () and standard deviation () described as follows:
In these equations, and are the mean and and are the standard deviation of the and distributions, while represents the correlation coefficient and represents the persistent sodium current conductance for the ith neuron. All parameters are given in Table 2.
Activation dynamics of
Request a detailed protocolHolographic stimulation was simulated by activating in small sets of randomly selected neurons across the preBötC population. Activation of this current was simulated by the following equation:
where represents the channel activation and ranges between 0 and 1, represents the decay time constant, and is the Kronecker delta function, which represents the instantaneous jump in from 0 to 1 at the time of stimulation (). Parameters were chosen such that the response in stimulated neurons matched those seen in Kam et al., 2013b. All parameters are given in Table 1.
Data analysis and definitions
Request a detailed protocolData generated from simulations was postprocessed in MATLAB (MathWorks, Inc). An action potential was defined to have occurred in a neuron when its membrane potential increased through . Histograms of population activity were calculated as the number of action potentials per bin per neuron, with units of . Network burst and burstlet amplitudes and frequencies were calculated by identifying the peaks and the inverse of the interpeak interval from the population histograms. The thresholds used for burst and burstlet detection were and , respectively. For the simulated holographic stimulation simulations, the start of a network burst was defined as the time at which the integrated preBötC population activity increased through the threshold for burst detection, while the end of a network burst was defined as the time at which the integrated preBötC activity returned to exactly zero.
Integration methods
Request a detailed protocolAll simulations were performed locally on an 8-core Linux-based operating system or on compute nodes at the University of Pittsburgh’s Center for Research Computing. Simulation software was custom written in C++. Numerical integration was performed using the first-order Euler method with a fixed step-size () of .
Data availability
All data generated or analysed during this study are included in the manuscript and supporting file; Source Data files have been provided for Figures 1-8. Simulation code has been published on GitHub at the following link. https://github.com/RyanSeanPhillips/Putting-the-theory-into-burstlet-theory (copy archived at swh:1:rev:fe56e63bc5839b7a6bd60f1fd0e22f8bec2a7669).
References
-
Synaptic depression and cortical gain controlScience (New York, N.Y.) 275:220–224.https://doi.org/10.1126/science.275.5297.221
-
Dynamics of ramping bursts in a respiratory neuron modelJournal of Computational Neuroscience 1:195.https://doi.org/10.1007/s10827-021-00800-w
-
Inhibition of endoplasmic reticulum CaNeuroscience 224:116–124.https://doi.org/10.1016/j.neuroscience.2012.08.016
-
Multi-timescale systems and fast-slow analysisMathematical Biosciences 287:105–121.https://doi.org/10.1016/j.mbs.2016.07.003
-
Models of respiratory rhythm generation in the pre-Bötzinger complex. I. Bursting pacemaker neuronsJournal of Neurophysiology 82:382–397.https://doi.org/10.1152/jn.1999.82.1.382
-
BookTheoretical Neuroscience: Computational and Mathematical Modeling of Neural SystemsComputational Neuroscience SeriesMIT Press.
-
Models of respiratory rhythm generation in the pre-botzinger complexIii. Experimental Tests of Model Predictions. Journal of Neurophysiology 86:59–74.https://doi.org/10.1152/jn.2001.86.1.59
-
Sodium and calcium current-mediated pacemaker neurons and respiratory rhythm generationThe Journal of Neuroscience 25:446–453.https://doi.org/10.1523/JNEUROSCI.2237-04.2005
-
Dendritic calcium activity precedes inspiratory bursts in preBotzinger complex neuronsThe Journal of Neuroscience 31:1017–1022.https://doi.org/10.1523/JNEUROSCI.4731-10.2011
-
Pontine mechanisms of respiratory controlComprehensive Physiology 2:2443–2469.https://doi.org/10.1002/cphy.c100015
-
Calcium currents of rhythmic neurons recorded in the isolated respiratory network of neonatal miceThe Journal of Neuroscience 18:10652–10662.https://doi.org/10.1523/JNEUROSCI.18-24-10652.1998
-
Facing the challenge of mammalian neural microcircuits: taking a few breaths may helpThe Journal of Physiology 593:3–23.https://doi.org/10.1113/jphysiol.2014.277632
-
Sodium and calcium mechanisms of rhythmic bursting in excitatory neural networks of the pre-Bötzinger complex: a computational modelling studyThe European Journal of Neuroscience 37:212–230.https://doi.org/10.1111/ejn.12042
-
Pacemaker behavior of respiratory neurons in medullary slices from neonatal ratJournal of Neurophysiology 72:2598–2608.https://doi.org/10.1152/jn.1994.72.6.2598
-
Inositol 1,4,5-trisphosphate (InsP3) and calcium interact to increase the dynamic range of InsP3 receptor-dependent calcium signalingThe Journal of General Physiology 110:529–538.https://doi.org/10.1085/jgp.110.5.529
-
Distinct Inspiratory Rhythm and Pattern Generating Mechanisms in the preBotzinger ComplexJournal of Neuroscience 33:9235–9245.https://doi.org/10.1523/JNEUROSCI.4143-12.2013
-
Ryanodine receptors: structure, expression, molecular details, and function in calcium releaseCold Spring Harbor Perspectives in Biology 2:a003996.https://doi.org/10.1101/cshperspect.a003996
-
Somatic Ca2+ transients do not contribute to inspiratory drive in preBötzinger Complex neuronsThe Journal of Physiology 586:4531–4540.https://doi.org/10.1113/jphysiol.2008.154765
-
Substitution of extracellular Ca2+ by Sr2+ prolongs inspiratory burst in pre-Bötzinger complex inspiratory neuronsJournal of Neurophysiology 113:1175–1183.https://doi.org/10.1152/jn.00705.2014
-
Phenomenological models of synaptic plasticity based on spike timingBiological Cybernetics 98:459–478.https://doi.org/10.1007/s00422-008-0233-1
-
Respiratory control by ventral surface chemoreceptor neurons in ratsNature Neuroscience 7:1360–1369.https://doi.org/10.1038/nn1357
-
Significance of extracellular potassium in central respiratory control studied in the isolated brainstem–spinal cord preparation of the neonatal ratRespiratory Physiology & Neurobiology 146:21–32.https://doi.org/10.1016/j.resp.2004.10.009
-
Pain-facilitating medullary neurons contribute to opioid-induced respiratory depressionJournal of Neurophysiology 108:2393–2404.https://doi.org/10.1152/jn.00563.2012
-
Dendritic A-Current in Rhythmically Active PreBötzinger Complex Neurons in Organotypic Cultures from Newborn MiceThe Journal of Neuroscience 38:3039–3049.https://doi.org/10.1523/JNEUROSCI.3342-17.2018
-
BookRespiratory rhythm generationIn: Miller AD, Bianchi AL, Bishop BP, editors. Neural Control of the Respiratory Muscles. CRC Press. pp. 119–130.
-
Synchronized Activity and Loss of Synchrony Among Heterogeneous Conditional OscillatorsSIAM Journal on Applied Dynamical Systems 1:146–174.https://doi.org/10.1137/S111111110240323X
-
TRPC channels and store-operated Ca2+ entryBiochimica et Biophysica Acta (BBA) - Molecular Cell Research 1793:223–230.https://doi.org/10.1016/j.bbamcr.2008.11.001
-
Brainstem respiratory networks: building blocks and microcircuitsTrends in Neurosciences 36:152–162.https://doi.org/10.1016/j.tins.2012.11.004
-
Two types of independent bursting mechanisms in inspiratory neurons: an integrative modelJournal of Computational Neuroscience 30:515–528.https://doi.org/10.1007/s10827-010-0274-z
-
Stabilization of bursting in respiratory pacemaker neuronsThe Journal of Neuroscience 23:3538–3546.https://doi.org/10.1523/JNEUROSCI.23-08-03538.2003
-
Timescales and Mechanisms of Sigh-Like Bursting and Spiking in Models of Rhythmic Respiratory NeuronsJournal of Mathematical Neuroscience 7:1–39.https://doi.org/10.1186/s13408-017-0045-5
-
Complex bursting dynamics in an embryonic respiratory neuron modelChaos (Woodbury, N.Y.) 30:043127.https://doi.org/10.1063/1.5138993
Article and author information
Author details
Funding
National Science Foundation (DMS1951095)
- Jonathan E Rubin
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Christopher Del Negro for his comments on a draft of this manuscript.
Copyright
© 2022, Phillips and Rubin
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,019
- views
-
- 177
- downloads
-
- 17
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Computational and Systems Biology
Degree distributions in protein-protein interaction (PPI) networks are believed to follow a power law (PL). However, technical and study bias affect the experimental procedures for detecting PPIs. For instance, cancer-associated proteins have received disproportional attention. Moreover, bait proteins in large-scale experiments tend to have many false-positive interaction partners. Studying the degree distributions of thousands of PPI networks of controlled provenance, we address the question if PL distributions in observed PPI networks could be explained by these biases alone. Our findings are supported by mathematical models and extensive simulations and indicate that study bias and technical bias suffice to produce the observed PL distribution. It is, hence, problematic to derive hypotheses about the topology of the true biological interactome from the PL distributions in observed PPI networks. Our study casts doubt on the use of the PL property of biological networks as a modeling assumption or quality criterion in network biology.
-
- Computational and Systems Biology
- Microbiology and Infectious Disease
Bacterial membranes are complex and dynamic, arising from an array of evolutionary pressures. One enzyme that alters membrane compositions through covalent lipid modification is MprF. We recently identified that Streptococcus agalactiae MprF synthesizes lysyl-phosphatidylglycerol (Lys-PG) from anionic PG, and a novel cationic lipid, lysyl-glucosyl-diacylglycerol (Lys-Glc-DAG), from neutral glycolipid Glc-DAG. This unexpected result prompted us to investigate whether Lys-Glc-DAG occurs in other MprF-containing bacteria, and whether other novel MprF products exist. Here, we studied protein sequence features determining MprF substrate specificity. First, pairwise analyses identified several streptococcal MprFs synthesizing Lys-Glc-DAG. Second, a restricted Boltzmann machine-guided approach led us to discover an entirely new substrate for MprF in Enterococcus, diglucosyl-diacylglycerol (Glc2-DAG), and an expanded set of organisms that modify glycolipid substrates using MprF. Overall, we combined the wealth of available sequence data with machine learning to model evolutionary constraints on MprF sequences across the bacterial domain, thereby identifying a novel cationic lipid.