Our first aim in this work was to test whether CICR from endoplasmic reticulum (ER) stores could repetitively amplify synaptically triggered Ca²⁺ transients. To address this aim, we constructed a cellular model that includes the ER. The model features a Ca²⁺ pump, which extrudes Ca²⁺ from the intracellular space, a sarcoendoplasmic reticulum calcium transport ATPase (SERCA) pump, which pumps Ca²⁺ from the intracellular space into the ER, and the Ca²⁺-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 Ca²⁺ current into the intracellular space with varied frequency and amplitude but fixed duration (250 ms) and monitored the resulting intracellular Ca²⁺ transients. Depending on parameter values used, we observed various combinations of low- and high-amplitude Ca²⁺ responses and characterized how the fraction of Ca²⁺ transients that have low amplitude depends on values selected within the 2D parameter space parameterized by Ca²⁺ pulse frequency and amplitude. We found that the fraction of low-amplitude Ca²⁺ transients decreases as either or both of the Ca²⁺ pulse frequency and amplitude are increased (Figure 1B and example traces C1–C4).

Figure 1

Download asset Open asset

Figure 1—source data 1

Calcium-induced calcium release.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig1-data1-v2.xlsx

Download elife-75713-fig1-data1-v2.xlsx

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 $I_{NaP}$-dependent intrinsically bursting neuron, with a burst frequency that is varied by injecting an applied current, $I_{APP}$ (Figure 2A2–A3). The rhythmic bursting from neuron 1 generates periodic postsynaptic currents ($I_{Syn}$) in neuron 2, carried in part by Ca²⁺ ions, which are analogous to the square wave Ca²⁺ current in Figure 1. The amplitude of the postsynaptic Ca²⁺ transient is determined by the number of spikes per burst (Figure 2A4) and by the parameter $P_{SynCa}$, which determines the percentage of $I_{Syn}$ carried by Ca²⁺ 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 $I_{CAN}$ activation.

Figure 2 with 1 supplement see all

Download asset Open asset

Figure 2—source data 1

Burstlets and bursts in a two-neuron network.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig2-data1-v2.xlsx

Download elife-75713-fig2-data1-v2.xlsx

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 $I_{APP}$ and $P_{SynCa}$. Increasing $I_{APP}$ 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 $I_{APP}$ 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 $I_{APP}$ corresponds to a horizontal cut through the panel). In general, increasing $P_{SynCa}$ 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 $P_{SynCa}$ and $I_{APP}$ on neuron 2 burst frequency (Figure 2B3). As $I_{APP}$ 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 $P_{SynCa}$ 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 $I_{APP}$ and $P_{SynCa}$ (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).

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 $K_{bath}$-dependent changes in the burstlet fraction characterized in Kallurkar et al., 2020. $K_{bath}$ 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 ($E_K$) and leak ($E_{Leak}$) reversal potentials as functions of $K_{bath}$ 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 $K_{bath}$. However, with increasing $K_{bath}$, intrinsic bursting could be maintained over a wide range of K⁺ concentrations when accompanied by increases in $g_{Leak}$ (Figure 3C). Additionally, the number of spikes per burst in the bursting regime increases with $K_{bath}$ (Figure 3—figure supplement 1). This $K_{bath}$-dependence of $g_{Leak}$ is consistent with experimental data showing that neuronal input resistance decreases with increasing $K_{bath}$ (Okada et al., 2005).

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 3—source data 1

Bath potassium concentration dependence of cellular and network dynamics.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig3-data1-v2.xlsx

Download elife-75713-fig3-data1-v2.xlsx

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 ($g_{NaP}$), as documented experimentally (Del Negro et al., 2002a; Koizumi and Smith, 2008). In both populations, we maintained the dependence of $g_{Leak}$ on $K_{bath}$ (see Figure 3C and E).

For the full preBötC network model, we first characterized the impact of changes in $K_{bath}$ on network behavior without calcium dynamics by setting $P_{SynCa}=0$. This network condition is analogous to in vitro preparations where all Ca²⁺ currents are blocked by Cd²⁺ 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 ${\displaystyle K_{bath}=5\mathrm{m}\mathrm{M}}$ (Figure 4A). Moreover, under these conditions, increasing $K_{bath}$ results in an increase in the burstlet frequency and amplitude (Figure 4B and C), which is consistent with experimental observations (Kallurkar et al., 2020).

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

Burstlets and Bursts in a larger network.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig4-data1-v2.xlsx

Download elife-75713-fig4-data1-v2.xlsx

In the full network with calcium dynamics (${\displaystyle P_{SynCa}>0}$), 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 Ca²⁺ 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 $K_{bath}$ and $P_{SynCa}$ (Figure 4D–F). In this case, the frequency of the postsynaptic Ca²⁺ oscillation is controlled by $K_{bath}$ (Figure 4B). However, because $K_{bath}$ also affects burstlet amplitude (Figure 4C), the postsynaptic Ca²⁺ amplitude is determined by both $K_{bath}$ and $P_{SynCa}$. If $K_{bath}$ is held fixed, then modulating $P_{SynCa}$ will only affect the amplitude of the postsynaptic Ca²⁺ transient since burstlet amplitude will not be impacted. The effects of selectively changing the postsynaptic Ca²⁺ 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 $I_{NaP}$; however, we obtain similar network results if we impose burstlet activity on the burstlet subpopulation and maintain the coupling between populations and Ca²⁺ dynamics for burst recruitment (Figure 4—figure supplement 1).

We found that increasing $P_{SynCa}$ or $K_{bath}$ 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 $K_{bath}$ or $P_{SynCa}$ is caused by the increase in the burstlet amplitude (Figure 4C) or in Ca²⁺ influx with each burstlet, respectively, both of which increase the Ca²⁺ transient in the pattern-generating subpopulation. The increase in burst frequency with increases in $K_{bath}$ or $P_{SynCa}$ 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 $K_{bath}$, by an increase in the burstlet frequency (Figure 4B). The slight increase in burst amplitude with increasing $K_{bath}$ 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 Ca²⁺ transients will increase with increasing $K_{bath}$ due to the change in the resting $V_m$ 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 $P_{SynCa}$ (i.e., the amplitude of the postsynaptic calcium transients) with $K_{bath}$ (Figure 4 traces G1–G4) generated $K_{bath}$-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 $K_{bath}$ levels).

Next, we further characterized the calcium dependence of the burstlet to burst transition in our model by simulating calcium blockade or $I_{CAN}$ blockade by a progressive reduction of $P_{SynCa}$ or $g_{CAN}$, respectively. We found that complete block of synaptically triggered Ca²⁺ transients or $I_{CAN}$ 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 Ca²⁺ transients increases the burstlet fraction by increasing the number of burstlets required to trigger a network burst, whereas $I_{CAN}$ 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 $I_{CAN}$ 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 Ca²⁺ block the decrease in amplitude results primarily from derecruitment (Figure 5E and F). These simulations provide mechanism-specific predictions that can be experimentally tested.

Figure 5

Download asset Open asset

Figure 5—source data 1

Simulated calcium or CAN current blockade.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig5-data1-v2.xlsx

Download elife-75713-fig5-data1-v2.xlsx

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).

Figure 6

Download asset Open asset

Figure 6—source data 1

Effects of simulated opioids on burstlets and bursts.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig6-data1-v2.xlsx

Download elife-75713-fig6-data1-v2.xlsx

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).

Figure 7 with 1 supplement see all

Download asset Open asset

Figure 7—source data 1

Simulated holographic stimulation.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig7-data1-v2.xlsx

Download elife-75713-fig7-data1-v2.xlsx

To determine the mechanisms involved in the refractoriness, we plotted the time courses of key slow variables in the model, namely, persistent sodium inactivation $h_{NaP}$, ER calcium (${[Ca]}_{ER}$), and synaptic depression $D$, 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 $h_{NaP}$ 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 ($P_{PP}$) 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 $W_{PP}$ to compensate for changes in $P_{PP}$ and keep the network synaptic strength, defined as $S=N_P\cdot P_{PP}\cdot W_{PP}$, at a constant value.

With this scaling, we found that decreasing/increasing $P_{PP}$ 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 $P_{PP}$, is indicated by a leftward/rightward shift in the probability vs. stimulation time curves relative to a control level of $P_{PP}$ ($P_{PP}=2\%$) (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.

Figure 8

Download asset Open asset

Figure 8—source data 1

Refractory period of evoked bursts following holographic stimulation.

https://cdn.elifesciences.org/articles/75713/elife-75713-fig8-data1-v2.xlsx

Download elife-75713-fig8-data1-v2.xlsx

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 $N_P$ and $W_{PP}$, and we take $P_{PP}$ to scale as $1/N_P$. The minimal number of inputs from active neurons needed to activate an inactive neuron depends on the synaptic weight, $W_{PP}$. Let $r$ denote this number for the specific value of $W_{PP}$ that we have selected. We can approximate the expected number of neurons receiving $r$ or more inputs from $A$ active neurons by computing the expected number receiving exactly $r$ inputs, which we denote as $[I_r]$, 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 $P_{PP}$ and halve $W_{PP}$, thus keeping their product constant. For this smaller $W_{PP}$, more inputs will be needed to activate an inactive neuron. Specifically, assume that now at least $2r$ inputs are needed for activation. The expected number of neurons receiving $2r$ inputs, $[I_{2r}]$, is given by

An elementary calculation shows that ${\displaystyle [I_{2r}]<[I_r]}$ for relevant parameter values (such as $N_P=300$ and small $r$ as indicated by the stimulation experiments). Thus, increasing $P_{PP}$ and proportionally scaling down $W_{PP}$ 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 $P_{PP}$ (Figure 8F). For a given stimulation time and number of neurons stimulated, however, decreasing $P_{PP}$ 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 $P_{PP}=2\%$. 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.

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 Ca²⁺ 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 Ca²⁺ 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 $K_{bath}$, the magnitude of the burstlet-triggered postsynaptic calcium transients must increase with $K_{bath}$. Some of this rise can result directly from the increase in burstlet amplitude with increasing $K_{bath}$ (see Kallurkar et al., 2020 and Figure 4C). To fully match the experimentally observed relationship between $K_{bath}$ and the burstlet fraction (Figure 4J), we also explicitly increased the parameter $P_{SynCa}$, which sets the proportion of the postsynaptic calcium current carried by Ca²⁺. Thus, our model predicts that the magnitude of postsynaptic Ca²⁺ transients triggered by EPSPs should increase as $K_{bath}$ is elevated.

This same prediction arises from considering the voltage-dependent properties of Ca²⁺ 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 $K_{bath}$. Specifically, it is likely that voltage-gated calcium channels are involved in generating the postsynaptic Ca²⁺ transients as dendritic Ca²⁺ transients have been shown to precede inspiratory bursts and to be sensitive to Cd²⁺, a calcium channel blocker (Del Negro et al., 2011). Consistent with this idea, Cd²⁺-sensitive Ca²⁺ 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 ${\displaystyle -65\mathrm{m}\mathrm{V}}$(Elsen and Ramirez, 1998), and importantly, the mean somatic resting membrane potential of non-rhythmogenic preBötC neurons increases from ${\displaystyle -67.034\mathrm{m}\mathrm{V}}$ to ${\displaystyle -61.78\mathrm{m}\mathrm{V}}$ when extracellular potassium concentration is elevated from ${\displaystyle 3\mathrm{m}\mathrm{M}}$ to ${\displaystyle 8\mathrm{m}\mathrm{M}}$ (Tryba et al., 2003). Moreover, at ${\displaystyle K_{bath}=9\mathrm{m}\mathrm{M}}$, 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 $K_{bath}$ is lowered (Okada et al., 2005). Putting together these data on resting membrane potential, EPSP sizes, and voltage-dependent activation of Ca²⁺ channels, we deduce that when ${\displaystyle K_{bath}=3\mathrm{m}\mathrm{M}}$, the magnitude of EPSPs may not reach voltages sufficient for significant activation of voltage-gated Ca²⁺ channels. As $K_{bath}$ is increased, however, both EPSC magnitudes and the membrane potential of pattern-forming neurons increase. Therefore, with increased $K_{bath}$, the prediction is that EPSCs will result in greater activation of voltage-gated Ca²⁺ channels and increased postsynaptic calcium transients. This effect is captured in the model by an increase in the parameter $P_{SynCa}$, which determines the percentage of the postsynaptic current carried by Ca²⁺ ions, with $K_{bath}$.

The idea that dendritic postsynaptic Ca²⁺ transients and $I_{CAN}$ activation play a critical role in regulating the pattern of preBötC dynamics is well supported by experimental and computational studies. Specifically, the dendritic Ca²⁺ 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 ($I_{CAN}$) (Mironov, 2008). Moreover, the rhythmic depolarization of otherwise non-rhythmogenic neurons (inspiratory drive potential) depends on $I_{CAN}$ (Pace et al., 2007a), while the inspiratory drive potential is not dependent on Ca²⁺ 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 $I_{CAN}$, 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 $I_{NaP}$-dependent rhythmogenic kernel (i.e., burstlet rhythm generator), and that rhythmic synaptic drive from these neurons triggers postsynaptic calcium transients, $I_{CAN}$ 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 Ca²⁺ and/or IP3 concentration may impact ER Ca²⁺ uptake and release dynamics. Store-operated Ca²⁺ 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.

In our model (Figure 4), a burstlet rhythm first emerges at a $K_{bath}$ 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 $K_{bath}$. Respiratory circuits appear to have some buffering capacity, however, such that for $K_{bath}$ concentrations below approximately 5 mM the extracellular K⁺ concentration remains elevated above $K_{bath}$ (Okada et al., 2005). The $K_{bath}$ 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 $K_{bath}$ 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 $K_{bath}$ and may impact the ramping shape of burstlet onset (Abdulla et al., 2021). Importantly, over the range of $K_{bath}$ values relevant both to experiments and our model, we find clear agreement on the dependence of burstlet fraction on $K_{bath}$ (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 ($I_h$) and a transient potassium current ($I_A$) 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 $I_{CAN}$ 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 $I_{CAN}$ 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, $I_h$.

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, 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 $K_{bath}$ is elevated above 9 mM (Del Negro et al., 2002c). These properties are not consistent with the burst and burstlet amplitudes or $K_{bath}$-dependent changes in the burstlet fraction seen experimentally (Kallurkar et al., 2020) and in our model.

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 Ca²⁺ channels by Cd²⁺ 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 ($I_{NaP}$). The role of $I_{NaP}$ in preBötC inspiratory rhythm generation is a contentious topic within the field, largely due to the inconsistent effects of $I_{NaP}$ block. We chose to use $I_{NaP}$ 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 $I_{NaP}$ (Phillips and Rubin, 2019b), (2) $I_{NaP}$-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, $I_{NaP}$-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 $I_{NaP}$ and $I_{CAN}$ 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.

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 $I_h$ and $I_A$, are blocked (see earlier discussion); (3) blocking postsynaptic Ca²⁺ transients will increase the burstlet fraction and decrease the burst amplitude before network bursts are eventually abolished (see Figure 5); (4) $I_{CAN}$ 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).

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 Ca²⁺ transients by CICR in order to activate $I_{CAN}$ and drive bursting in the rest of the network. Moreover, the burstlet fraction depends on rate of Ca²⁺ buildup in intracellular stores, which is impacted by $K_{bath}$-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.

Model 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 ${\displaystyle C=36\mathrm{p}\mathrm{F}}$ is the membrane capacitance and each $I_i$ represents a current, with i denoting the current’s type. The currents include the action potential generating Na⁺ and delayed rectifying K⁺ currents ($I_{Na}$ and $I_K$), persistent Na⁺ current ($I_{NaP}$), voltage-gated Ca²⁺ current ($I_{Ca}$), Ca²⁺-activated nonselective cation (CAN) current ($I_{CAN}$), K⁺-dominated leak current ($I_{Leak}$), synaptic current ($I_{Syn}$), μ-opioid receptor-activated G protein-coupled inwardly rectifying K⁺ leak current ($I_{GIRK}$) (Kubo et al., 1993), and a holographic photostimulation current ($I_{Holo}$). $I_{APP}$ denotes an applied current injected from an electrode. The currents are defined as follows:

where g_i is the maximum conductance, $E_i$ is the reversal potential, and m_i and h_i are gating variables for channel activation and inactivation for each current $I_i$. The glutamatergic synaptic conductance $g_{Syn}$ is dynamic and is defined below. The values used for the g_i and $E_i$ are mostly shown in Table 1 , with a few conductances selected from distributions as indicated in Table 2.

Activation (m_i) and inactivation (h_i) of voltage-dependent channels are described by the following differential equation:

where $X_{\mathrm{\infty }}$ represents steady-state activation/inactivation and ${\tau }_X$ is a time constant. For $I_{Na}$, $I_{NaP}$, and $I_{Ca}$, the functions $X_{\mathrm{\infty }}$ and ${\tau }_X$ take the forms

The values of the parameters ($X_{1/2}$, $k_X$, ${\tau }_{max}^X$, ${\tau }_{1/2}^X$, and $k_{\tau }^X$) corresponding to $I_{Na}{I}_{NaP}$, ${\displaystyle I_{NaP}}$ and $I_{Ca}$ are given in Table 1.

For the voltage-gated potassium channel, the steady-state activation $m_{\mathrm{\infty }}^K(V)$ and time constant ${\tau }_m^K(V)$ are given by the expressions

where

The values for the constants $A_{\alpha }$, $A_{\beta }$, $B_{\alpha }$, $B_{\beta }$, $k_{\alpha }$, and $k_{\beta }$ are also given in Table 1.

$I_{CAN}$ activation depends on the $C{a}^{2+}$ concentration in the cytoplasm (${[Ca]}_{Cyto}$) and is given by

The parameters $C{a}_{1/2}$ and $n$ represent the half-activation Ca²⁺ concentration and the Hill coefficient, respectively, and are included in Table 1.

The dynamics of ${[Ca]}_{Cyto}$ is determined in part by the balance of Ca²⁺ efflux toward a baseline concentration via the Ca²⁺ pump and Ca²⁺ influx through voltage-dependent activation of $I_{Ca}$ and synaptically triggered Ca²⁺ transients, with a percentage ($P_{SynCa}$) of the synaptic current ($I_{Syn}$) carried by Ca²⁺ ions. Additionally, the model includes an intracellular compartment that represents the ER, which impacts ${[Ca]}_{Cyto}$. The ER removes Ca²⁺ from the cytoplasm via a sarcoplasmic/endoplasmic reticulum Ca²⁺- ATPase (SERCA) pump, which transports Ca²⁺ from the cytoplasm into the ER ($J_{SERCA}$) and releases Ca²⁺ into the cytoplasm via calcium-dependent activation of the inositol triphosphate (IP3) receptor ($J_{IP3}$). Therefore, the dynamics of ${[Ca]}_{Cyto}$ is described by the following differential equation:

where ${\displaystyle {\alpha }_{Ca}=2.5\cdot {10}^{-5}\mathrm{m}\mathrm{M}/\mathrm{f}\mathrm{C}}$ is a conversion factor relating current to rate of change of ${[Ca]}_{Cyto}$, ${\displaystyle {\tau }_{pump}=500\mathrm{m}\mathrm{s}}$ is the time constant for the Ca²⁺ pump, ${\displaystyle C{a}_{min}=5.0\cdot {10}^{-6}\mathrm{m}\mathrm{M}}$ is a minimal baseline calcium concentration, and ${\alpha }_{ER}=2.5\cdot {10}^{-5}$ is the ratio of free to bound Ca²⁺ in the ER.

The flux of Ca²⁺ from the ER to the cytoplasm through the IP3 receptor is modeled as

where ${\displaystyle E{R}_{leak}=0.1/\mathrm{m}\mathrm{s}}$ represents the leak constant from the ER stores, ${\displaystyle G_{IP3}=77,500/\mathrm{m}\mathrm{s}}$ represents the permeability of the IP3 channel, ${\displaystyle K_a=1.0\cdot {10}^{-4}\mathrm{m}\mathrm{M}}$ and ${\displaystyle K_l=1.0\cdot {10}^{-3}\mathrm{m}\mathrm{M}}$ are dissociation constants, and ${\displaystyle [IP3{]}_i=1.5\cdot {10}^{-3}\mathrm{m}\mathrm{M}}$ is the cytoplasm IP3 concentration. Finally, the Ca²⁺-dependent IP3 gating variable, $l$, and the Ca²⁺ concentration in the ER, ${[Ca]}_{ER}$, are determined by the following equations:

where ${\displaystyle A=0.1\mathrm{m}\mathrm{M}/\mathrm{m}\mathrm{s}}$ is a conversion factor, ${\displaystyle K_d=0.2\cdot {10}^{-3}\mathrm{m}\mathrm{M}}$ is the dissociation constant for IP3 inactivation, ${[Ca]}_{total}$ is the total intracellular calcium concentration, and ${\sigma }_{Ca}=0.185$ is the ratio of cytosolic to ER volume. The total intracellular calcium concentration is described as

Finally, removal of Ca²⁺ from the cytoplasm by the SERCA pump is modeled as

where ${\displaystyle G_{SERCA}=0.45\mathrm{m}\mathrm{M}/\mathrm{m}\mathrm{s}}$ is the maximal flux through the SERCA pump, and ${\displaystyle K_{SERCA}=7.5\cdot {10}^{-5}\mathrm{m}\mathrm{M}}$ is a dissociation constant.

Nondimensionalization of similar models in past work (Wang and Rubin, 2017; Wang and Rubin, 2020) has shown that $h_{NaP},l$, and ${[Ca]}_{ER}$ are the slowest variables in the model and evolve on similar timescales, while ${[Ca]}_{Cyto}$ 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 ${(g_{Syn})}_i$ of the ith target neuron is described by the following equation:

where $g_{Tonic}$ 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, $W_{j,i}$ represents the weight of the synaptic connection from neuron $j$ to neuron i, $D_j$ is a scaling factor for short-term synaptic depression in the presynaptic neuron $j$ (described in more detail below), $C_{j,i}$ is an element of the connectivity matrix ($C_{j,i}=1$ if neuron $j$ makes a synapse with neuron i and $C_{j,i}=0$ otherwise), $H(.)$ is the Heaviside step function, and $t$ denotes time. ${\tau }_{Syn}$ is an exponential synaptic decay constant, while $t_{j,n}$ is the time at which the nth action potential generated by neuron $j$ 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 ($D_j$) 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 $D_0=1$ sets the maximum value of $D_j$, ${\displaystyle {\tau }_D=1000\mathrm{m}\mathrm{s}}$ sets the rate of recovery from synaptic depression, ${\alpha }_D=0.2$ sets the fractional depression of the synapse each time neuron $j$ spikes, and $\delta (.)$ is the Kronecker delta function that equals 1 at the time of each spike in neuron $j$ and 0 otherwise. Parameters were chosen to qualitatively match data from Kottick and Del Negro, 2015. Note that with this choice of ${\tau }_D$ 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 $W_{1,2}=W_{2,1}=0.006$ and $C_{1,2}=C_{2,1}=1$. For the full preBötC population model comprising rhythm- and pattern-generating subpopulations, the weights of excitatory conductances were uniformly distributed such that $W_{j,i}=U(0,W_{Max})$ where $W_{Max}$ 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 $C_{j,i}$ values (see Table 3). Effects of opioids on synaptic transmission for source neurons in the rhythmogenic subpopulation (Figure 6) were simulated by scaling $W_{j,i}$ with the parameter ${\gamma }_{\mu OR}$, which ranged between 0 and 0.5 and sets the percent synaptic block.

The 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 $20\pm 9\%$ of preBötC inspiratory neurons are active during burstlets at ${\displaystyle K_{Bath}=9\mathrm{m}\mathrm{M}}$. 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 ($SS{T}^{-}$) and -positive ($SS{T}^{+}$), respectively (Cui et al., 2016; Ashhad and Feldman, 2020). Somatostatin-positive neurons are estimated to comprise 72.6% of the $Dbx{1}^{+}$ 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 $N=400$ neuron preBötC population (i.e., $N_R=100$ and $N_P=300$). The rhythm- and pattern-generating neurons are distinguished by their $I_{NaP}$, $I_{Leak}$, and $I_{CAN}$ conductances. Also, we included the K⁺ leak current $I_{GIRK}$ 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, ${\displaystyle P_{RR}=13\mathrm{\%}}$ and ${\displaystyle P_{PP}=2\mathrm{\%}}$, 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 $P_{RP}=P_{PR}=30\%$ such that the total connection probability in the network is approximately 13% (Rekling et al., 2000).

Heterogeneity was introduced by normally distributing the parameters $g_{leak}$, $g_{NaP}$, and $g_{CAN}$ as well as uniformly distributing the weights ($W_{j,i}$) of excitatory synaptic connections (see Table 2 and Table 3). Additionally, $g_{leak}$ was conditionally distributed with $g_{NaP}$ 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 $g_{NaP}$ 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 $g_{leak}$ given $g_{NaP}$ is itself a normal distribution with mean (${\mu }_{Leak}^{*}$) and standard deviation (${\sigma }_{Leak}^{*}$) described as follows:

In these equations, ${\mu }_{Leak}$ and ${\mu }_{NaP}$ are the mean and ${\sigma }_{Leak}$ and ${\sigma }_{NaP}$ are the standard deviation of the $g_{Leak}$ and $g_{NaP}$ distributions, while ${\displaystyle \rho =0.8}$ represents the correlation coefficient and $g_{NaP}^i$ represents the persistent sodium current conductance for the ith neuron. All parameters are given in Table 2.

Holographic stimulation was simulated by activating $I_{Holo}$ in small sets of randomly selected neurons across the preBötC population. Activation of this current was simulated by the following equation:

where $m_{Holo}$ represents the channel activation and ranges between 0 and 1, ${\tau }_{Holo}$ represents the decay time constant, and $\delta (.)$ is the Kronecker delta function, which represents the instantaneous jump in $m_{Holo}$ from 0 to 1 at the time of stimulation ($t_{stim}$). 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 generated from simulations was postprocessed in MATLAB (MathWorks, Inc). An action potential was defined to have occurred in a neuron when its membrane potential $V_m$ increased through ${\displaystyle -35\mathrm{m}\mathrm{V}}$. Histograms of population activity were calculated as the number of action potentials per ${\displaystyle 20\mathrm{m}\mathrm{s}}$ bin per neuron, with units of $APs/(s\cdot neuron)$. 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 ${\displaystyle 30\mathrm{s}\mathrm{p}\mathrm{k}/\mathrm{s}/\mathrm{N}}$ and ${\displaystyle 2.5\mathrm{s}\mathrm{p}\mathrm{k}/\mathrm{s}/\mathrm{N}}$, 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.

All 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 ($\mathrm{\Delta }t$) of ${\displaystyle 0.025\mathrm{m}\mathrm{s}}$.