Introduction

The neural control of respiration is a complex physiological process that requires precise coordination among multiple brainstem nuclei. One such nucleus, the preBötzinger Complex (preBötC), plays a crucial role in generating the inspiratory rhythm and pattern (Gray, Janczewski, Mellen, McCrimmon, & Feldman, 2001; Smith, Ellenberger, Ballanyi, Richter, & Feldman, 1991; Tan et al., 2008; Wang et al., 2014). In contrast, the lateral parafacial region (pFL), also referred to as pFRG in the initial studies, has emerged as an important neuronal structure containing neurons that are responsible for the generation of active expiration via the recruitment of expiratory musculature, such as the abdominal (ABD) muscles (de Britto & Moraes, 2017; Huckstepp, Cardoza, Henderson, & Feldman, 2015; Pagliardini et al., 2011; Pisanski & Pagliardini, 2019). Although the involvement of preBötC in respiration, as well as its origin, anatomical location, markers and anatomical projections has been extensively studied (Biancardi et al., 2023; Del Negro, Funk, & Feldman, 2018; Gray et al., 2001; Hayes et al., 2017; Tan et al., 2008; Wang et al., 2014; Yackle et al., 2017; Yang & Feldman, 2018; Yang, Kim, Callaway, & Feldman, 2020), the expiratory oscillator remains a challenging structure to locate accurately due to the absence of a definite anatomical marker and the lack of its activation during experimental resting conditions (de Britto & Moraes, 2017; Magalhães et al., 2021; Pagliardini et al., 2011). Understanding the precise location and functional properties of the expiratory oscillator is crucial for unraveling the complete neural circuitry underlying respiratory control.

Previous research has employed a variety of tools to study and elucidate this expiratory oscillator that include pharmacological disinhibition and excitation (Boutin, Alsahafi, & Pagliardini, 2017; de Britto, Magalhães, da Silva, Paton, & Moraes, 2020; de Britto & Moraes, 2017; Korsak, Sheikhbahaei, Machhada, Gourine, & Huckstepp, 2018; Magalhães et al., 2021; Pagliardini et al., 2011; Zoccal et al., 2018), as well as chemogenetic and optogenetic approaches (Huckstepp et al., 2015; Pagliardini et al., 2011; Pisanski, Ding, Koch, & Pagliardini, 2019). The coordinates used for the latter studies considered the core of the expiratory oscillator proximal to the caudal tip of the facial nucleus (VIIc) (-0.2 to +0.5 mm from VIIc). Interestingly, the use of chemogenetics to inhibit pFL at the level of the caudal tip of the facial nucleus (-0.3 to +0.3 mm from VIIc) partially inhibited the putative expiratory muscle output, with ABD recruitment still occurring naturally during sleep, albeit at a diminished rate (Pisanski et al., 2019). Similarly, chemogenetic inhibition of pFL at +0.5 mm from VIIc decreased the intensity of the ABD signals obtained in response to bicuculline/strychnine injections, but did not silence it entirely (Huckstepp et al., 2015). Although both of these elegant chemogenetic studies have extensively contributed to our understanding of the pFL, the existing evidence suggests that the expiratory oscillator may expand beyond the limits of the viral expression achieved in said studies, as proposed by Huckstepp et al., (2015). Intriguingly, other research (Silva, Oliveira, Souza, Moreira, & Takakura, 2019) located its core at more rostral coordinates (+0.3 to +1.0 mm from VIIc). This dichotomy in the coordinates used to study the functional properties of the expiratory oscillator demonstrates that the characterization of the anatomical and functional boundaries of pFL are incomplete.

To address this gap, we used focal injections of bicuculline along different rostrocaudal locations of the brainstem as a tool to create a functional map of the pFL. Bicuculline, a competitive antagonist of γ-aminobutyric acid type A (GABA-A) receptors, has been widely employed to study neuronal excitability and functional circuitry in various brain regions, including the pFL (Pagliardini et al., 2011). By strategically administering localized volumes of bicuculline at multiple rostrocaudal levels of the ventral brainstem, we aimed to selectively enhance the excitability of neurons driving active expiration, thereby revealing the extension of the pharmacological response and the most efficient site in generating active expiration. Because all the locations studied have been previously shown to produce active expiration (Pagliardini et al., 2011; Boutin et al., 2017; de Britto et al., 2017; Huckstepp et al., 2015; Silva et al., 2019) we conducted a more holistic analysis of all the respiratory changes induced by bicuculline injections. This analysis included traditional respiratory measures such as: tidal volume, minute ventilation, respiratory rate and oxygen metabolism on top of the abdominal changes induced by active expiration. Additionally, using a novel multidimensional cycle-by-cycle analysis specifically developed for this study, we characterized the effect that bicuculline elicited on the different phases of the respiratory cycle at each injection site, as the differences in area under the curve of the airflow, ᶴDIA EMG and ᶴABD EMG, as well as the combined differences of these three respiratory signals using a novel phase-plane analysis. This approach enabled us to construct a comprehensive functional map of the pFL in conjunction with anatomical immunostaining techniques.

Our results indicate that the injection of bicuculline produced active expiration at all locations studied (-0.2 mm to +0.8 mm from VIIc). However, the strongest effects in terms of changes in tidal volume (V_T), minute ventilation (V_E) and differences in combined respiratory signals were observed at the two most rostral locations (+0.6 mm and +0.8 mm). Similarly, ABD activation lasted longer at these sites and the swiftest onset of the ABD response was observed at the +0.6 mm location. Interestingly, our multivariate analysis of the respiratory cycle permitted further differentiation of the response elicited in the two rostral locations, with the strongest deformations of the respiratory loop during late-Expiration (late-E) and post-Inspiration (post-I) with injections at the +0.8 mm locations, whereas injections at the level of +0.6 mm producing more pronounced deformations of the respiratory cycle during the inspiratory phase. The use of the multi dimensional trajectory map will certainly enhance our understanding of the role that the expiratory oscillator plays in respiratory control. Additionally, it will provide a crucial foundation for future investigations into its modulation, plasticity, and potential therapeutic targets for respiratory disorders.

Methods

Experimental Subjects

Thirty-five adult male Sprague Dawley rats weighing 340.4 g ± 13.2 were used for the study. The rats were housed in a controlled environment with a 12-hour light-dark cycle and were allowed unrestricted access to food and water. All experimental procedures followed the guidelines set by the Canadian Council of Animal Care and received approval from the Animal Care and Use Committee (ACUC) of the University of Alberta (AUP#461).

Surgical Preparation

Prior to the experiment, the rats were initially anesthetized using 5% isoflurane in air for induction, followed by 1-3% isoflurane for maintaining a surgical plane of anesthesia. During this time, we cannulated the femoral vein to facilitate the gradual administration of urethane (1.5-1.7 g/kg body weight) for inducing permanent and irreversible anesthesia. We assessed the depth of anesthesia by monitoring the absence of the withdrawal reflex. We cannulated the trachea and used a flow head connected to a transducer (DP10310N154D, GM Instruments, UK) to detect respiratory flow using a carrier demodulator (CD15, Validyne Engineering, CA, US). We provided supplemental oxygen (30%) throughout the experiment and connected a gas analyzer (ML 206-1704, ADInstruments, CO, US) to the tracheal tube. The gas analyzer measured fractional concentration of O₂. Based on this and the flow rate at the level of the trachea (minute ventilation), we calculated O₂ consumption according to Depocas and Hart (1957). Paired electromyogram (EMG) wire electrodes (AS633, Cooner Wire, CA, US) were inserted into the oblique ABD and diaphragm (DIA) muscles. The wires were connected to differential amplifiers (model 1700, AM Systems, WA, US), and we sampled the activity at a rate of 1kHz using the Powerlab 16/30 system (ML880, AD Instruments, CO, US). We performed vagotomy by resecting a 2 mm portion of the vagus nerve at the mid-cervical level, and the rats’ body temperature was maintained at a constant level of 37±1°C using a servo-controlled heating pad (55-7022, Harvard Apparatus, MA, US). During the entire experimental procedure, rats breathed spontaneously and end tidal CO2 was not adjusted through the experimental protocol.

Bicuculline Injections

Following the surgical preparation, we placed the instrumented rats in a prone position on a Kopf stereotaxic frame (model 962, Kopf Instruments, CA, US). We trimmed and aseptically cleaned the incision area and achieved access to the brainstem by partially removing the occipital bone. We determined the coordinates for injection into the ventral medulla with Bregma positioned 5 mm below lambda. The bicuculline/fluorobeads solution (200 µM in HEPES, Bicuculline methochloride, D46M4623V, Sigma-Aldrich, MA, US; n=28) or HEPES/fluorobead (F8797, Invitrogen, MA, US) (CTRL, n=7) was bilaterally pressure-injected into specific coordinates using a sharp glass micropipette (30 µm tip diameter) at a rate of 100 nL/min (total volume injected was 200 nL per side). Each rat was injected at 2.8 mm lateral from the midline and at a specific rostrocaudal (RC) coordinate based on the following groups: -0.2 mm from the caudal tip of the facial nucleus (VIIc) (n=5), +0.1 mm from VIIc (n=7), +0.4 mm from VIIc (n=5), +0.6 mm from VIIc (n=6), +0.8 mm from VIIc (n=5), and CTRL (n=7). We recorded the physiological responses to the injection for 20-25 min.

Histology Procedures

Upon completion of the experiments, we transcardially perfused the rats with saline (0.9% NaCl), followed by 4% paraformaldehyde in phosphate buffer. We then collected, postfixed, cryoprotected, and sectioned the brains in a cryostat (CM1950, Leica, Germany) to obtain 30 µm thick transverse slices. To prepare the brain tissue for immunohistochemistry, we washed the brain sections multiple times with phosphate- buffered saline (PBS) to remove the cryoprotectant solution (15% sucrose, 15% ethyleneglycol, 0.5% PVP-40 in 0.5X-PBS). All immunostaining procedures were conducted at room temperature. We incubated the sections in a blocking solution (0.3% Triton X-100 and 10% normal donkey serum (NDS) in PBS) for 1 hour to enhance membrane permeability to antibodies and minimize nonspecific binding. After blocking, we incubated the tissue overnight with the primary antibody solution, which consisted of PBS, 1% NDS, and 0.3% Triton X-100. In order to identify specific brain structures relevant to this study, we used the following primary antibodies: anti-choline acetyltransferase (ChAT; goat; 1:800; RRID: AB_20797, EMD Millipore, ON, Canada), anti-PHOX2B (PHOX2B; mouse; 1:100; RRID:AB_2813765, Santa Cruz Biotechnology, TX, US), anti-cFos (cFos; rabbit; 1:500; RRID:AB_2247211, Cell Signaling Technology, MA, US), and anti-tyrosine hydroxylase (TH; chicken; 1:800; RRID:AB_570923, EMD Millipore, ON, Canada). The next day, we washed the tissue three times in PBS and incubated it for 2 hours in a secondary antibody solution containing PBS, 1% NDS, and specific secondary antibodies (1:200; Cy3-donkey anti-rabbit, RRID:AB_2307443; Cy5- donkey anti-mouse, RRID:AB_2340820; Cy2-donkey anti-chicken, RRID:AB_2340370; rhodamine-red-donkey anti-goat, RRID:AB_2340423; Jackson ImmunoResearch Laboratories, PA, US). Following the incubation with the secondary antibodies, we washed the sections three times and mounted them with Fluorsave mounting medium (345789-20ML, EMD Millipore, ON, Canada). We observed the slides under an Evos FL fluorescent microscope (AMF4300, Thermofisher Scientific, MA, US) and acquired TIFF files for cell counting and analysis of the rostrocaudal, mediolateral, and dorsoventral coordinates of the injection sites using ImageJ. We examined serial sections with a 120 µm interval, spanning from 400 µm caudal to 1000 µm rostral to the caudal tip of the facial nucleus (VIIc). To group the animals according to the location of the injections, we calculated a 3D distance from the VIIc using the rostrocaudal, mediolateral, and dorsoventral coordinates of each injection site.

Data Acquisition and Analysis

EMG signals were amplified at ×10,000 and filtered between 300 Hz and 1 kHz (model 1700, AM-systems, WA, US). The data was sampled at a rate of 1 kHz (using an automatic anti-aliasing filter at 500Hz) using the PowerLab 16/35 acquisition system and analyzed using LabChart7 Pro (ADInstruments, CO, US), Excel 2013, Origin 9 (OriginLab Corp., Northampton, MA, USA), as well as custom scripts written in MATLAB (The Mathworks Inc.). All raw EMG data were digitally rectified and integrated using a time-constant decay of 0.08s. Respiratory airflow was further used to determine respiratory rate, V_T, and Ve. The tidal volume was obtained by integrating the airflow amplitude and converting it to milliliters using a five-point calibration curve (0.5–5 ml range). The coupling between ABD and DIA signals was measured as a ratio and analyzed by quantifying the number of bursts of activity observed for the ABD and DIA EMG signals during the first 10 minutes of the response, excluding time bins at end of the response (due to fading and waning of the ABD activity in those instances). The ABD delay was calculated as the time between the end of the second injection and the beginning of the ABD response (negative values mean that the response started following the first injection before the second injection occurred). Whereas ABD duration was calculated by measuring the time between the first and the last ABD burst elicited by bicuculline.

Airflow and integrated EMG signals (Figure 1 A) were standardized across rats by first zeroing signals in the rest period phases between breaths as recorded during baseline conditions and then converting the remaining amplitude values into standard deviation units (SD) based on a time-collapsed amplitude distribution. Baseline data was averaged over a period of 5 minutes. Experimental data was tracked in time bins of 2- minute duration from the baseline period prior to injections and spanned 20 min of recording post-injection. Mean-cycle measurements for each signal were computed by averaging across all cycles within a given time bin (∼300 cycles in baseline, ∼100 cycles per response time bin). We then used the average calculations of respiratory rate (RR), tidal volume (VT), Minute Ventilation (Ve), expiratory ABD amplitude, expiratory ABD area, VO₂, VE/VO₂ to obtain values relative to the baseline period. Peak responses were identified as the time bin that produced the strongest changes relative to baseline. To assess the changes within each of the phases of the breathing cycle, each mean- cycle was subdivided into the Late-E, inspiratory, Post-I, and resting phases using the zero-crossings of the mean airflow signal (Green lines, Figure 1B). The inspiratory phase begins with an inward (positive going) intake of air as airflow increases above zero and ends when airflow crosses zero towards the negative direction. This negative going deflection begins the Post-I phase, which lasts until airflow returns to zero during the rest period. The Late-E phase is defined in two ways. In cases where a significant negative deflection in airflow was observed previous to inspiration, the late-E phase was determined to have begun when airflow decreased below zero and ended when airflow crossed back above zero to begin the subsequent inspiratory phase. When airflow did not decrease below zero prior to inspiration, as was often the case during baseline conditions, the late-E phase was determined as the last quarter of the expiratory period (Figure 1B). The Area under the curve (AUC) was measured during baseline and was subtracted from the corresponding AUC of the response for each time bin (Figure 1C). This AUC measure was computed as the sum of the signal in a given respiratory phase as all signals were sampled at the same rate. Note that areas calculated below the zero (0) line, as would be expected from a negative airflow during expiration, yields negative AUC values. Any changes to the duration of the three respiratory phases, as well as changes to the entire respiratory period were independently assessed using the triggers described above.

Figure 1.

We further contrasted the overall changes to the three respiratory signals elicited across injection locations by combining them into a multi-dimensional (3D) space (Figure 1D). In a manner similar to respiratory pressure-volume loops, each point in time is defined as a dot in this 3D space based on the 3 simultaneous physiological measurements of respiration we recorded (airflow, together with the ∫DIA EMG, and ∫ABD EMG). Each respiratory cycle could then be plotted as a trajectory of points in this space, represented as a loop. Example loops for the baseline and the first 2-minute time bin of the elicited response are shown in Figure 1D and are also projected onto separate 2D planes in Figure 1E. Our approach was to measure how the trajectories of response loops changed relative to the baseline loop during each phase of the respiratory cycle, as per the color scheme (same as in Figure 1B). Given these loops may be deformed as a result of changes to any one of the three signals which form the basis of this 3D space, this method determines the total extent to which breathing is altered in response to stimulation. We quantified changes to respiratory trajectories elicited by bicuculline injection using the sum of the “straight-line” or “Euclidean” distances between points on the response and baseline loops in each respiratory phase. An example of the Euclidean distance between the peak of inspiration in baseline and response loops is shown in Figure 1F (pink line). To compute this distance for all time points, we resampled the response to ensure that that baseline and response mean-cycles had the same number of comparison phase points in a given respiratory phase. Distances were then calculated between these baseline and response loops at each timepoint t using the Euclidean distance formula:

where AirDiff_t, DiaDiff_t and AbdDiff_t are the differences between the response and baseline mean-cycles at timepoint t. Averaging the distances measured at every timepoint in a given respiratory phase yields the mean Euclidean distance measure used in the analyses displayed in Figure 6G.

We extended our multivariate analysis of respiration to account for the cycle-by-cycle variability in the recorded signals. At each timepoint, we compared how far away the response mean-cycle is from the baseline mean-cycle relative to the covariance across all baseline cycles. Thus, this measure can be interpreted as having units of standard deviation, as it describes how likely it is to observe the response given the distribution of baseline measurements gathered across many respiratory cycles. This measure is known as the Mahalanobis distance and is exemplified as the pink distance relative to the black cloud of points in the baseline mean-cycle of Figure 1F. Similarly to above, this distance is computed for each timepoint t:

where R_t is the response mean-cycle at timepoint t, µ_t is the baseline mean-cycle at timepoint t, and ∑_t is the covariance of all baseline cycles at timepoint t. Similar to the Euclidean distance above, we the calculated Mahalanobis distances across all timepoints with each respiratory phase and then averaged them across cycles to assess the mean Mahalanobis distance.

Statistics

We tested the assumptions of homogeneity of variance and normal distribution using the Brown–Forsythe (α = 0.05) and Shapiro–Wilk (α = 0.05) tests, respectively. Comparisons between injection locations were carried out across time bins using a two- way repeated-measures analysis of variance (two-way RM ANOVA; one repeated factor; α = 0.05) where appropriate. In case of significant main effects, we conducted pairwise multiple comparisons post-hoc procedures (Bonferroni method; α = 0.05).

Analysis of the peak response across injection locations, was done using an analysis of variance (one-way ANOVA; α = 0.05) in the case of absolute values, or a one-way RM ANOVA (α = 0.05) in the case of values reported relative to baseline conditions. Upon finding significant main effects, we used Tukey or Bonferroni post hoc test (α = 0.05) respectively. Effect size was computed within the analyses of variance using η² (Eta squared). In cases where assumptions of normality were rejected (p < 0.05), we compared responses across injection location via a Kruskal-Wallis test and following with a post-hoc Dunn test, adjusting the alpha value for multiple comparisons with Sidak’s method.

Results

Histological Analysis

For each rat, we precisely identified the core of the injection sites by assessing the rostro-caudal coordinates obtained from sections containing fluorobeads. Based on the core coordinates, experiments were divided into 5 groups in which injection spanned from -0.2 mm caudal to the tip of the facial nucleus (VIIc) to +0.8 mm rostral to VIIc (as illustrated in Fig 2 B). All injections were consistently positioned ventrolateral to the facial nucleus, and there was no overlap observed between the fluorobead-marked injection sites and PHOX2B+ cells within the more ventro-medial retrotrapezoid nucleus (RTN), as depicted in Figures 2 A-B. Additionally, the injection sites were confirmed not to overlap with the soma of TH+ cells in either the caudal C1 or the more rostral A5 areas, as shown in Figure 2 A-B.

Figure 2.

To assess the cellular responses surrounding the injection sites, we quantified both cFos+ and cFos+/PHOX2B+ cells within a defined region surrounding the fluorobead- marked areas in the pFL (as outlined in Figure 2 A). In rats injected with HEPES buffer (control group), we observed an average activation of 44.7 ± 4.0 cells per hemisection along the rostrocaudal axis of the ventral medulla (n=7). Among these activated cells, an average of 2.3 ± 1.0 cells per hemisection were also positive for PHOX2B, as illustrated in the inset of Figure 2 C. In contrast, rats injected with bicuculline displayed a distinct response pattern, characterized by a peak activation with an average of 104.5 ± 5.1 cFos+ cells localized at the core of each injection site (-0.2 mm = 89.7; +0.1mm = 123.1; +0.4mm = 110.2; +0.6mm = 98.3; +0.8mm = 101.2; Figure 2C). Importantly, the number of cFos+ cells decreased to 44.8 ± 1.2 per section beyond the boundaries of the injection area for each experimental group. On the other hand, the number of cFos+/PHOX2B+ cells remained consistent between the control (CTRL) animals (2.3 ± 1.0 cFos+/PHOX2B+ cells per hemisection) and those injected with bicuculline (2.4 ± 0.4 cFos+/PHOX2B+ cells). These findings strongly suggest that bicuculline specifically activated cells within the vicinity of the injection sites which spread ∼300 μm (Figure 2C, horizontal lines) and did not activate PHOX2B+ cells in the RTN area, beyond their baseline level of activity.

Characteristics of ABD signals across the rostrocaudal axis of pFL

We successfully elicited ABD responses with the administration of bicuculline at all tested injection sites. The ABD response was most reliably triggered at locations rostral to the tip of VIIc. This observation is supported by the fact that the caudal location (-0.2 mm from VIIc) elicited responses in only 3 out of 5 rats, whereas all other groups consistently induced responses in all rats (n=5-7 rats /group). The elicited response was characterized by a late-E component in the ∫ABD EMG trace across all experimental groups (Figure 2D), as well as a downward deflection in airflow, signifying forced exhalation of the respiratory reserve volume preceding inspiration (i.e., active expiration; Figure 2D). Remarkably, in the most rostral groups (+0.6 and +0.8 mm rostral to VIIc), the ∫ABD EMG traces exhibited a tonic expiratory component, which preceded the late- E peak and was absent in the most caudal locations (-0.2 mm, +0.1 mm, and +0.4 mm). Additionally, 5 out of 6 rats in the group with injection sites at +0.6 mm rostral to VIIc, exhibited a post-I ABD peak, a feature that was absent in all other groups (Figure 2D).

Temporal dynamics of the ABD response reveal a longer-lasting and fastest response at rostral locations

The temporal characteristics of the ᶴABD EMG signal elicited by bicuculline injections at different rostrocaudal locations exhibited notable variations. Overall, responses had a shorter duration in the most caudal locations (-0.2 mm; +0.1 mm; + 0.4 mm) compared to the rostral locations (+0.6 mm; +0.8 mm; Figure 3A). In all groups, the ABD response initiated following the second bicuculline injection, except at location +0.6 mm, where it commenced following the first injection and became more pronounced with the second injection (Figure 3 A). Across all groups, the most robust ABD responses, as measured by ∫ABD EMG amplitude, were observed within the first 2 min of the response, declining to zero in the last 2 min (18 min post-injection) in the most caudal groups (-0.2 mm; +0.1 mm; + 0.4 mm). In contrast, these responses remained present, albeit weaker, in the most rostral groups (+0.6 mm; +0.8 mm) after 20 min from injection (Figure 3 B-C).

Figure 3.

The coupling of the ᶴABD EMG and ᶴDIA EMG peak signals, which reflects the robustness of the response and the coupling between inspiration and active expiration, was notably weaker at the most caudal location (-0.2 mm = 0.6 ± 0.2) in comparison to the most rostral groups (+0.4 mm = 0.96 ± 0.02; +0.6 mm = 0.89 ± 0.05; +0.8 mm = 0.97 ± 0.02; One-Way ANOVA p = 0.015; η² = 0.43; Tukey -0.2 mm vs +0.4 mm: p = 0.024; -0.2 mm vs +0.6 mm: p = 0.048; -0.2 mm vs +0.8 mm: p = 0.029; Figure 3 D). Similarly, the ABD response duration was longer at the two most rostral locations (+0.6 mm = 17.6 ± 2.7 min; +0.8 = 17.1 ± 3.3 min) compared to the most caudal group (-0.2 mm = 2.4 ± 1.1 min; One-Way ANOVA p = 0.043; η² = 0.41; Tukey -0.2 mm vs +0.6 mm: p = 0.048; -0.2 mm vs +0.8 mm: p = 0.041; Figure 3E). Lastly, the group with injection sites located at +0.6 mm initiated responses sooner than all other groups (-0.2 mm = 20.3 ± 13.4 sec; +0.1 mm = 32.5 ± 20.6 sec; +0.4 = 40.1 ± 28.7 sec; +0.8 = 23.1 ± 19.8 sec), with an average response delay of 88.7 ± 32.3 sec previous to the second injection (but after the first one), which is represented as negative values (One-Way ANOVA p= 0.041; η² = 0.57; Tukey -0.2 mm vs +0.6 mm: p = 0.039; +0.1 mm vs +0.6 mm: p = 0.040; +0.4 mm vs +0.6 mm: p = 0.045; +0.8 mm vs +0.6 mm: p = 0.049; Figure 3 F). In summary, these results underscore that injections at the most rostral locations (+0.6 mm and +0.8 mm) elicit more robust and prolonged ABD responses than the caudal location (-0.2 mm). Furthermore, bicuculline injections at +0.6 mm from VIIc initiated the fastest response, suggesting proximity to the cells responsible for ABD recruitment following bicuculline injection.

Bicuculline injection elicited stronger respiratory effects at the rostral locations

The ABD response elicited by bicuculline injection generated a late-E downward inflection in airflow that was absent in baseline conditions (Figure 4 A). We measured the late-E peak amplitude and area to assess the strength of the response across different rostro-caudal locations. On average, the late-E peak amplitude and area values reached a maximum between 2 and 4 min post-second injection for all rostro- caudal locations except the most rostral injection site (+0.8 mm), which peaked between 6 and 8 min post-second injection (Figure 4 B,D). Although the maximum late-E peak amplitude in the two most rostral groups (+0.6 mm = -0.033 ± 0.007 V; +0.8 mm = -0.027 ± 0.011 V), were not different from the remaining locations (-0.2 mm= -0.014 ±0.004 V; +0.1 mm = -0.021 ± 0.007 V; +0.4 mm = -0.023 ± 0.004 V; One-Way ANOVA p= 0.41; Figure 4 C), the large effect size (η² = 0.18) would suggest that further investigation is warranted. Similarly, the late-E peak area in the two most rostral locations (+0.6 mm = -0.0056 ± 0.0010 V.s; +0.8 mm = -0.0047 ± 0.0019 V.s), were not different from the rest of the injection sites (-0.2 mm= -0.0014 ± 0.0007 V.s; +0.1 mm = - 0.0035 ± 0.0015 V.s; +0.4 mm = -0.0038 ± 0.0011 V.s; One-Way ANOVA p= 0.30; Figure 4 E), but there was a large effect size (η² = 0.2).

Figure 4.

The injection of bicuculline along the rostrocaudal axis of the ventral medulla induced a drop in respiratory frequency in all the injection sites tested (Figure 4 F), but this drop in frequency was only significant in comparison to baseline in the two most caudal groups (One-Way repeated measures ANOVA p = 0.003; η² = 0.46; Bonferroni: -0.02 mm baseline vs injection, p = 0.03; +0.1 mm baseline vs injection, p = 0.005; Figure 4 G). The respiratory rate reached a minimum value of 37.08 ± 1.36 bpm (12% drop from baseline) at 2 min post-injection in the most caudal location (-0.2 mm), at 4 min for the middle groups (+0.1 mm = 37.2 ± 2.0 bpm – 17% drop from baseline; +0.4 mm = 40.1 ± 1.6 bpm – 9% drop from baseline) and at 6 min for the most rostral injection sites (+0.6 mm = 38.6 ± 2.6 bpm – 10% drop from baseline; +0.8 mm = 39.1 ± 2.2 bpm – 11% drop from baseline; Figure 4 F). Interestingly, the duration of inspiration during the response decreased in all groups relative to baseline respiration (T_i response = 0.279 ± 0.034s, T_i baseline = 0.318 ± 0.043s, Wilcoxon rank sum: Z = 3.24, p = 0.001). Contrary to this decrease in inspiratory duration, the total expiratory time increased in all groups and remained elevated compared to baseline (T_E response = 1.313 ± 0.188s, T_E baseline = 1.029 ± 0.161s, Wilcoxon rank sum: Z = 4.49, p = 0.001). Overall, these results suggest that the most caudal locations (-0.2 mm and +0.1 mm) experienced a faster and more drastic drop in respiratory frequency following bicuculline injection, and this reduction in frequency was driven by an increase in expiratory period accompanied by a reduction in the inspiratory time.

The amount of inspired air per each breath, measured as V_T, increased when compared to baseline in all the injection sites tested except the most caudal location (-0.2 mm from VIIc) (One-Way repeated Measures ANOVA, p = 0.01; η² = 0.44 Figure 4 H-I). In all the groups, the increase in V_T peaked at 4 min post-injection (Figure 4 H). Of particular interest is that the peak V_T observed in the +0.6 mm location (V_T = 10.8 ± 0.5 ml/kg – 29% increase from baseline) was higher than that observed in the two most caudal groups (-0.2 mm = 7.0 ± 0.4 ml/kg – 8% increase from baseline and +0.1 mm = 8.0 ± 0.3 ml/kg – 16% increase from baseline; One-Way ANOVA p = 0.01; Bonferroni p = 0.03; Figure 4 I). On the other hand, V_E decreased compared to baseline in the most caudal group (-0.2 mm = 255.4 ± 16.2 ml x min^-1 x kg^-1 – 12% drop from baseline, One- Way Repeated Measures ANOVA, p= 0.001; η² = 0.59), whereas at the +0.6 mm location, V_E increased following bicuculline injection (V_E = 414.3 ± 22.7 – 16% increase from baseline, One-Way Repeated Measures ANOVA, p = 0.001; η² = 0.59 Figure 4 J- K). Interestingly, the change produced in V_E after injection of bicuculline at the most rostral locations (+0.4 mm; +0.6 mm and +0.8 mm) was larger than that produced after injection at -0.2 mm from VIIc (One-Way ANOVA p = 0.001, η² = 0.59; Bonferroni p = 0.0004; Figure 4 K). These results suggest that the V_E was mostly driven by the drop in respiratory frequency at the two most caudal locations (-0.2 mm and +0.1 mm), whereas V_T drove the increase in V_E observed in the most rostral groups (+0.4 mm; +0.6 mm; +0.8 mm).

Injection of bicuculline decreased oxygen consumption at the rostral locations

To evaluate the metabolic effects of the respiratory changes induced by bicuculline injections, we measured the amount of oxygen consumption (V_O2). Overall, V_O2 remained similar to baseline at the 3 most caudal locations (-0.2 mm; +0.1 mm and +0.4 mm; Figure 4 L-M) and dropped from baseline only at the most rostral sites (+0.6 mm and +0.8 mm) (One-Way Repeated Measures ANOVA p = 0.0001; η² = 0.59). The minimum V_O2 achieved at the rostral groups (+0.6 mm = 11.8 ± 0.8 ml x min^-1 x kg^-1 – 33% drop from baseline and +0.8 mm = 12.2 ± 6.1 ml x min^-1 x kg^-1 – 25% drop from baseline) was also different from the values achieved at the most caudal locations (-0.2 mm = 17.3 ± 3.3 ml x min^-1 x kg^-1 – 10% drop from baseline; +0.1 mm = 17.6 ± 1.3 ml x min^-1 x kg^-1 – 17% drop from baseline; One-Way Repeated Measures ANOVA p = 0.0001, η² = 0.59, Bonferroni p = 0.005; Figure 4 M). Interestingly, the ratio of V_E/V_O2 increased only at the two most rostral locations(One-Way Repeated Measures ANOVA p = 0.001; η² = 0.55; Figure 4 N-O), peaking at 6-8 min post-second injection (+0.6 mm = 34.8 ± 3.0 – 73% increase from BL; +0.8 mm = 38.7 ± 12.8 – 82% increase from baseline; Figure 4 N). In summary, these results suggest that bicuculline injection and the resulting respiratory effects led to a reduction in oxygen consumption in the most rostral locations. Furthermore, the reduction in V_O2 paired with the increases in V_E observed at the most rostral injection sites, resulted in hyperventilation, an effect which was not observed in the two most caudal groups.

Rostral bicuculline injection induces more prominent changes to all phases of the respiratory cycle

To further differentiate the responses elicited by bicuculline injections, we subdivided our analysis by respiratory phase (late-E, inspiratory, and post-I) and measured the differences in the area of respiratory signals (airflow, ∫DIA EMG and ∫ABD EMG) as per the algorithm described in the methods section (Figure 5 A-C).

Figure 5.

Figure 6.

Injections at all locations induced negative deflections in the airflow signal (Figure 5D) and were concomitant with a potent increase of the ∫ABD EMG signal within the late-E phase (Figure 5 D,F). The expiratory airflow reached significance (α = 0.05) when comparing the second most rostral location (+0.6mm) relative to control injections at time points from 4 min (KW test: H(5) = 14.66, p = 0.012, Dunn: p < 0.001) to 12 min post-injection (KW test: H(5) = 18.35, p = 0.003, Dunn: p < 0.001; Figure 5D). Moreover, airflow responses were more negative following bicuculline injections in the +0.6mm group compared to the -0.2mm group, and remained different at several time points from 8 min (KW test: H(4) = 11.41, p = 0.022, Dunn: p < 0.001) to the last 20 min of the post-injection period (KW test: H(4) = 10.65, p = 0.031, Dunn: p = 0.002; Figure 5D). In a similar manner, increases in the ᶴABD EMG at the 2 most rostral locations (+0.6 mm: Purple, +0.8 mm: Orange) differed relative to those elicited by control injections and at the most caudal injection site from 8 min to 14 min postinjection (-0.2 mm: Red) (KW test: H(5) = 20.81, p <0.001, Dunn: +0.6mm: p < 0.001; +0.8mm: p < 0.001; Figure 5F) with the +0.8mm group remaining elevated at 16 min (KW test: H(4) = 13.56, p = 0.009, Dunn: +0.8: p = 0.002; Figure 5F). Thus, our results in the Late-E phase demonstrate that the rostral regions of the ventral medulla are the most sensitive to bicuculline- induced increases in expiratory airflow, presumably elicited via a prominent activation of ABD activity.

For the inspiratory phase, bicuculline injections evoked strong increases in the respiratory airflow, which were strongest in the most rostral locations (+0.6mm, +0.8mm) and relatively weaker in more caudal locations (-0.2mm). Inspiratory airflow in the +0.6 location increased relative to control injections for time points extending from 2 min (KW test: H(4) = 13.25, p = 0.021, Dunn: p < 0.001) to 12 min post-injection (KW test: H(4) = 14.58, p = 0.012, Dunn: p = 0.001; Figure 5D). Furthermore, we also observed increases in the airflow of the +0.6 mm location when compared to the +0.4 mm group, from 10 min (KW test: H(4) = 10.50, p = 0.033, Dunn: p = 0.004) to 12 min (KW test: H(5) = 10.93, p = 0.027, Dunn: p = 0.003) post-injection (Figure 5D). Contributions of the ABD muscle to the respiratory changes during inspiration (Figure 5F, inspiratory phase) were determined to be negligible as the raw ABD EMG is not active during this phase and the remaining signal in the ᶴABD EMG can be accounted for by exponential decay induced by the integration algorithm used to pre-process these signals (see methods section). These results illustrate that the increases in inspiratory airflow evoked by bicuculline injections are stronger when targeting more rostral injection locations. However, these changes are relatively weak and, on average, remain altered for a shorter period of time relative to the effects noted above for the Late-E phase.

Negative deflections in the airflow signal during the post-I phase were observed in several responses and were most pronounced in the 2 most rostral injection groups consistent with the results above (Purple and Orange, Figure 5D). Bicuculline injections into the most rostral group (+0.8mm: Orange) evoked substantial increases in ∫ABD EMG activity during the post-I phase which were elevated during the 6 minute (KW test: H(4) = 10.37, p = 0.035, Dunn: p = 0.002) and 16 minute time bins (KW test: H(4) = 9.52, p = 0.049, Dunn: p = 0.001) when compared to injections into the most caudal group (-0.2mm). These effects were unlikely to be related to any changes in the DIA signal, as the DIA muscle activity was quiescent during the post-I period (Figure 5A) and the ᶴDIA EMG is contaminated by the pre-processing as noted above for the ABD muscle during inspiration (see methods section).

Altogether, the area-under-the-curve analysis in Figure 5 shows that airflow is substantially altered in all 3 phases of the respiratory cycle and the ABD muscle is strongly activated in several groups in the late-E and post-I periods with injections of bicuculline along the pFL axis. Importantly, these effects are stronger when targeting the rostral areas of the pFL (+0.6mm and +0.8 mm) compared to injections that target the more caudal regions.

Bicuculline responses following injections along the rostrocaudal axis of the ventral medulla can be differentiated by the deviations of their 3D trajectories at different phases of the respiratory cycle

We next sought to determine whether responses to bicuculline injections across the rostro-caudal axis of the pFL could be more precisely differentiated from each other when combining the responses in all three recorded signals (airflow, together with ʃDIA and ʃABD EMG) evoked within the periods of late-E, inspiration, and post-I. Our approach was to represent all three signals as axes in a 3-dimensional space wherein each point corresponds to an individual time point during our recordings. When plotting the measurements across an entire respiratory cycle, this multidimensional representation resembles a looped trajectory as shown in Figure 1D for cycles recorded during the baseline and post-injection periods. Indeed, each phase of the cycle (late-E, inspiration, and post-I) could be recognized as a particular region along the 3D trajectory. Inspiration, for example, occupies the space where airflow and ∫DIA EMG are positive (Figure 1D: red). In baseline conditions, the loop is largely confined to the plane defined by the airflow-∫DIA axes, as the ABD muscles are not activated during resting states (Figure 1E, left). Following bicuculline injections, respiratory loops became more complex, involving marked perturbations occurring in this first plane, together with novel deformations extending along the ∫ABD EMG axis as shown in the right-hand plane in Figure 1E, described by airflow and the ʃABD EMG. The relative amplitudes of these complex deviations from the baseline cycle could then be used to characterize the degree of change caused by our injections across pFL sites.

Figures 6 A through E demonstrate how the trajectories computed during the baseline (shown in black) are altered as the response evolves across each time bin spanning the full 20-minute post-injection period for all injection locations. Supplementary videos of each example recording in Figure 6, along with their cycle-by-cycle variability, further clarify their 3D shape as the trajectory is rotated about the ∫DIA EMG axis (Supplementary Materials 1-5). We noted several key deformations to baseline loops including a “bulb” of activation on the outermost edge of the loop during the peak of inspiration when airflow is at its maximum. These “bulbs” are similarly confined to the ᶴDIA EMG-airflow axis as both of these signals are involved in shaping their overall prominence. Baseline loops are also perturbed during the late-E period as a result of ᶴABD EMG-airflow activation which can be seen as “tails” which extend upwards in the right-most plots of Figure 6A-E. Abdominal-airflow activation during the post-I period manifests as an upwards “foot-like” extension of the cycle in left-most portion of the Figure 6 A-E, most prominently displayed in the representative recording of Figure 6D. The example trajectories in Figure 6A-E demonstrate how the deformations in each phase can be used as features to distinguish the responses of each injection location.

Consistent with our analyses above, the late-E tails of the response loops elicited by caudal injection locations are relatively smaller compared to those observed in rostral locations (Figure 6A). Moreover, these “tails” in the -0.2mm, +0.1mm, and +0.4mm groups last for a shorter duration than the two most rostral locations as can be determined by their more blue-yellow color (Figure 6A-E: rightmost panel). Inspiratory bulbs were more similar across groups and were most prominent in the +0.6mm group (Figure 6A-E: middle panel). Lastly, post-I “feet” are largely absent in the example responses of the three most caudal groups (-0.2mm, +0.1mm, +0.4mm) and are particularly noticeable in the +0.6mm and +0.8mm groups.

Quantifying the extent of these 3 features as the Euclidean distance between bicuculline injection and baseline in each respiratory phases is displayed in Figure 6F. Compared to control injections of HEPES buffer, bicuculline evoked pronounced Late-E “tails”, as well as inspiratory “bulbs”, and post-I “feet” in all injection groups which lasted throughout the 20 min response, with the notable exception of the most caudal group (-0.2mm) which returned to a baseline trajectory by 8 min post-injection (Figure 6G). Consistent with the results in Figure 5, the most rostral injections at +0.6mm and +0.8mm elicited the largest late-E deformations. The 3D approach applied here further differentiated the responses of these two groups, as the late-E “tails” of the +0.6mm group were longer than those observed in the -0.2mm group early in the response from 8 min (KW test: H(4) = 12.22, p = 0.016, Dunn: p < 0.001) to 12 min (KW test: H(4) = 12.34, p = 0.015, Dunn: p = 0.002), while those in the +0.8mm group peaked later in the response and were longer than in the -0.2mm group from 12 min (KW test: H(4) = 12.34, p = 0.015, Dunn: p = 0.001) to 16 min (KW test: H(4) = 10.89, p = 0.028, Dunn: p = 0.001; Figure 6G). These results are consistent with the observations in Figure 3F, which indicate that responses in the +0.6mm group begin earlier relative to other groups. Moreover, the +0.6mm group produced the largest deformations to the inspiratory “bulb” which were greater than those in the -0.2mm group from 8 min (KW test: H(4) = 12.59, p = 0.013, Dunn: p < 0.001) to 14 min (KW test: H(4) = 10.27, p = 0.036, Dunn: p = 0.001; Figure 6G). On the contrary, the largest post-I “feet” were observed in the +0.8mm group and were greater than the -0.2mm group from 4 min (KW test: H(4) = 10.01, p = 0.040, Dunn: p = 0.002) to 16 min (KW test: H(4) = 9.70, p = 0.046, Dunn: p = 0.001), excluding at the 10 minute (KW test: H(4) = 8.93, p = 0.063) and 14 minute time bins (KW test: H(4) = 8.10, p = 0.088; Figure 6G).

As demonstrated in the representative graphs in Figure 6A-E, the variability of the respiratory trajectories within the baseline or response periods (shown as clouds around the trajectory “loop”), differs between experiments. To account for the impact of this variability on our quantification of trajectory deformations, we computed the mean Mahalanobis distance within the late-E, inspiratory, and post-I phases (see methods). Our results revealed a gradient of response intensities relative to respiratory variability, with the two most rostral groups producing the strongest deformations to breathing trajectories (Figure 6H). Similar to our results above, these groups differed in the time course and relative strength of their responses in each phase. Injections at +0.6mm produced more pronounced late-E “tails” compared to caudal injection responses from 6 min (KW test: H(4) = 11.80, p = 0.019, Dunn: p = 0.001) for 12 min (KW test: H(4) = 12.24, p = 0.016, Dunn: p = 0.001; Figure 6H). Late-E “tails” from the +0.8mm group were more prominent than “tails” evoked by caudal bicuculline injections for much greater proportion of the 20 min response compared to the +0.6mm group, ranging from 6 min (KW test: H(4) = 11.80, p = 0.019, Dunn: 0.002) to 20 min post-injection (KW test: H(4) = 11.01, p = 0.026 Dunn: p < 0.001; Figure 6H). Deformations to the inspiratory “bulb” relative to respiratory variability were larger in the +0.6mm group compared to caudal responses and peaked at the 12minute time bin post-injection (KW test: H(4) = 10.63, p = 0.031, Dunn: p = 0.002)(Figure 6H). In contrast, the extent of post-I “feet” was greater in the +0.8mm group than the -0.2mm group from 6 min (KW test: H(4) = 11.00, p = 0.027, Dunn: p = 0.002) to 20 min (KW test: H(4) = 10.74, p = 0.030, Dunn: p < 0.001), where as the extent of post-I “feet” in the +0.6mm group was greater than in caudal injections only at the 8 minute (KW test: H(4) = 11.65, p = 0.020, Dunn: p = 0.002) and 12 minute time bins (KW test: H(4) = 12.23, p = 0.016, Dunn: p = 0.001). Thus, while response intensity in all respiratory periods was strongest in the two most rostral injection groups, overall changes in the late-E and post-I period were strongest in the +0.8mm group, while changes in inspiration were strongest in the +0.6mm group.

Discussion

In this study, our primary aim was to investigate and quantify the extent of the rostro- caudal distribution of the source of active expiration within the pFL. To accomplish this objective, we used bicuculline, a GABA-A receptor antagonist known for its ability to disinhibit pFL cells and generate active expiration (de Britto & Moraes, 2017; Huckstepp et al., 2015; Pagliardini et al., 2011; Silva et al., 2019). Bicuculline was administered at various rostro-caudal coordinates, and we methodically assessed the resulting respiratory responses to pinpoint the location within this region that induced the most significant and relevant changes to the respiratory cycle. Our findings offer valuable insights into the neural circuitry governing active expiration in the brainstem, with a particular focus on less-explored, more rostral areas of the pFL.

Technical considerations

Previous work has already demonstrated a temporal correlation between the onset of late-E neuron activity in the pFL and ABD activity recruitment in response to bicuculline (Pagliardini et al., 2011; de Britto and Moraes, 2017; Magalhaes et al., 2021), as well as the presence of GABAergic sIPSCs in late-E neurons (Magalhaes et al., 2012). Since bicuculline has previously shown to exhibit inhibitory effects on Ca²⁺ activated K⁺ currents inducing non-specific potentiation of NMDA currents (Johnson and Seutin, 1997), caution should be warranted in attributing our findings solely to the effects on GABAergic disinhibition as activation of late-E neurons and active expiration may be also be triggered by these secondary pharmacological effects of bicuculline.

Our experiments focused on determining the area in the pFL that is most effective in generating active expiration as measured by ABD EMG activity and expiratory flow. We did not attempt in recording single cell neuronal activity at various locations as previously shown in other studies (Pagliardini et al 2011; Magalhaes et al., 2021) as this approach would most likely find some late-E neurons across the pFL and thus not effectively discriminate between areas of the pFL. Future studies involving multi-unit recordings or imaging of cell population activities will help to determine the firing pattern and population density of bicuculline-activated cells and further determine differences in distribution and function of late-E neurons across the region of the pFL.

Histological Analysis

We adopted a stereological approach to identify the central region of the injection sites by visualizing sections containing fluorobeads. Experimental subjects were categorized based on the proximity of the injection sites, spanning from -0.2 mm to +0.8 mm from VIIc. We took measures to ensure that injections were positioned ventrolateral to the facial nucleus (VIIc) and did not overlap with PHOX2B+ cells within the more ventromedial retrotrapezoid nucleus (RTN). The assessment of cellular activity, quantified through cFos staining, unveiled the existence of basal activity in control rats. This observed baseline activity is likely emanating from subthreshold physiological processes within the parafacial area which do not culminate in ABD activity. Analysis of the cFos staining confirmed focal activation of neurons in the pFL of rats injected with bicuculline and minimal cFos expression in the PHOX2B+ cells in all groups as compared to the control group. These results confirm the very limited mediolateral spread of the drug from the core site of injection and support previous findings support the hypothesis that the majority of PHOX2B+ cells are more ventrally located in the ventral parafacial area (pFV, Huckstepp et al., 2015) and PHOX2B+ cell recruitment is not necessary for active expiration (de Britto & Moraes, 2017; Magalhães et al., 2021).

Characteristics of ABD Signals

At rest, respiratory activity does not present with active expiration (i.e, expiratory flow below its functional residual capacity in conjunction with expiratory-related ABD muscle recruitment) and expiratory flow occurs due to passive recoil of chest wall with little to none contribution of abdominal muscle activity. Active expiration and abdominal recruitment can be spontaneously observed during sleep (in particular REM sleep, Andrews and Pagliardini, 2015; Pisanski et al., 2019) and can be triggered during increased respiratory drive (e.g. Hypercapnia, RTN stimulation, Abbott et al., 2011). Although never assessed in freely moving, unanesthetized rodents, bicuculline has been extensively used to generate active expiration and late-E neuron activity in both juvenile and adult anesthetized rats (Pagliardini et al., 2011; Huckstepp et al., 2015 Huckstepp et al., 2016; Huckstepp et al., 2018; De Britto and Moraes, 2017; Magalhaes et al., 2021). Our efforts successfully elicited expiratory ABD responses and active expiration at all tested injection sites, consistent with previous reports (Boutin et al., 2017; de Britto et al., 2020; de Britto & Moraes, 2017; Huckstepp et al., 2015; Korsak et al., 2018; Magalhães et al., 2021; Pagliardini et al., 2011; Silva et al., 2019; Zoccal et al., 2018). Previous work has primarily examined the core of pFL at locations ranging from -0.2 mm to +1.0 mm from the caudal tip of the facial nucleus, a range similar to the one explored in the current study (-0.2 mm to +0.8 mm from VIIc). However, it is important to note that those studies have traditionally used the presence or absence of active expiration in response to a stimulus as the only measure to assess active expiration. Our study reveals that all locations within that range in which injections were made exhibited active expiration with various intensities. Consequently, we employed a combination of diverse respiratory measures to assess the injection site that induced the most robust respiratory effect.

Remarkably, the ABD response was consistently generated at locations rostral to VIIc, aligning with previous observations that did not find a reliable triggering of the ABD response at caudal positions (Pagliardini et al., 2011). This response was characterized by a late-E component at all rostro-caudal locations and, in the most rostral groups, an additional post-I component that was absent in the injections at more caudal locations. Although the configuration of the ABD response, whether late-E only, post-I, or a combination of the two, has not been extensively discussed in the literature, these distinctions may provide valuable insights into the involvement of distinct neuronal circuits within the pFL and in connection with other structures of the respiratory network in mediating the ABD response. For example, sample traces from the injection of bicuculline/strychnine from 0 to +1.0 mm from VIIc with 4-5% CO₂ exposure exhibited a tonic ABD response throughout the entire expiratory phase (without post-I or late-E peaks) (Silva et al., 2019), whereas sample traces from injections spanning from -0.2 mm to +0.4 mm from VIIc seemed to indicate a combination of tonic and late-expiratory ABD responses (de Britto & Moraes, 2017; Pagliardini et al., 2011). Interestingly, the latter studies also observed tonic increases in ABD EMG (without a late-expiratory peak), but only during baseline conditions and not in response to drug injection. In our study, we also observed a post-I ABD activation that only appeared at the most rostral locations. This post-I peak was previously observed in neonate rats during ABD activation occurring throughout apneic periods induced by systemic fentanyl injection (Janczewski, Onimaru, Homma, & Feldman, 2002). A similar biphasic pattern (late- E/post-I) was also observed in facial nerve recordings of in vitro preparations and associated with neuron activity in the parafacial area (Onimaru & Homma, 2003; Onimaru, Kumagawa, & Homma, 2006). Therefore, although it has previously been described, it is currently unknown the exact mechanism by which this post-I activity in the ABD muscles is generated. For example, the interplay between the rostral pFL and brainstem structures generating post-inspiratory activity, such as the proposed post- inspiratory oscillator (PiCo; Anderson et al., 2016) or pontine respiratory networks could be reasonably involved in this process. The presence of these unique response featured at different rostro-caudal locations suggests a complex interplay among neural populations within the ventral medulla, necessitating further investigation beyond the scope of the present study.

Temporal dynamics of the ABD response reveal a longer-lasting and shortest latency response at rostral locations

The temporal dynamics of the ABD response exhibited significant variations along the rostro-caudal axis. Responses at rostral locations displayed prolonged durations and stronger coupling with DIA signals compared to caudal locations. Previous research employing bicuculline injections to induce active expiration have reported an increase in the coupling between ABD and DIA signals +0.1mm to +0.3 mm from VIIc in juvenile rats (de Britto & Moraes, 2017). Similarly, we observed good coupling at locations rostral to +0.1 mm. Such coupling has been documented in other studies as an indicator of the robustness of the ABD response, either during hypercapnia/sleep or chemogenetic modulation during sleep (Leirão, Silva, Gargaglioni, & da Silva, 2018; Pisanski et al., 2019). Therefore, in our study, we used it as a measure of the strength of the achieved activation of active expiration. Additionally, prior studies have indicated that the duration of effects induced by bicuculine/strychnine injections ranges from 16 to 30 min (Pagliardini et al., 2011; Silva et al., 2019). However, these studies did not specify whether there were disparities in the duration of ABD responses along the rostro-caudal axis. In our study, we have conclusively established that there are indeed variations in the robustness (as denoted by ABD/DIA coupling) and duration of the ABD responses, with the more rostral locations yielding the most potent responses in this regard. Importantly, the group with injection sites at +0.6 mm from VIIc exhibited the swiftest response onset, suggesting that this area is the most critical for the generation of active expiration, either through direct activation of the expiratory oscillator or, alternatively, for providing a strong tonic drive to late-E neurons located elsewhere.

Bicuculline injection elicited stronger respiratory effects at rostral locations in pFL

Bicuculline injections induced significant alterations in a number of respiratory parameters, reminiscent of previous studies exploring the effects of bicuculline/strychnine injections in the ventral medulla (de Britto & Moraes, 2017; Pagliardini et al., 2011; Silva et al., 2019). We observed a drop in respiratory frequency concomitant with the emergence of ABD recruitment in the caudal locations, similar to what was observed in those locations through pharmacological activation of pFL (Boutin et al., 2017; Pagliardini et al., 2011). However, in the most rostral locations in our study, the respiratory frequency remained unchanged. This dichotomy in the frequency response despite the presence of ABD recruitment in both cases may suggest that the slowing of the respiratory period in the most caudal locations may be a consequence of disinhibition of the adjacent Bötzinger complex. Additional evidence to support this view is that the nadir of respiratory frequency was achieved at increasing periods of time as the injections were located increasingly more rostrally, with the most caudal injection reaching the minimum respiratory frequency at 2 min post-injection, while at the most rostral injection site, the minimum respiratory frequency was observed at 6 min.

On the other hand, V_T increased by 16-29% in all the studied locations except -0.2 mm, similar to what was observed previously with activation of pFL at locations ranging from -0.2 mm to +0.4 mm (Boutin et al., 2017; Pagliardini et al., 2011). However, minute ventilation increased above baseline at the +0.6 mm location in our study, whereas it decreased at the most caudal location and remained unchanged in the rest of the sites. These results are in agreement with previous observations, where V_E remained unchanged at caudal locations after bicuculline injection (Pagliardini et al., 2011). However, it is important to highlight the absence of studies that have measured V_T and V_E after eliciting active expiration at more rostral locations (Silva et al., 2019). These results emphasize that the rostro-caudal organization of pFL is not as simple as previously thought and particularly highlights the importance of studying more rostral locations that have previously been neglected.

Injection of bicuculline decreased oxygen consumption at the rostral locations

In terms of metabolic consequences of bicuculline-induced respiratory changes, we observed a decrease in V_O2 at the most rostral locations, where we observed the largest changes in V_T and V_E produced by bicuculline injections. These changes in V_T and V_E were not due to physiological needs and, therefore, could be considered a form of artificial hyperventilation that would cause a drop in arterial pCO₂ (Daly Md Fau - Hazzledine & Hazzledine, 1963) and affect the gas exchange processes in the blood and peripheral tissues to reduce V_O2. Alternatively, the abdominal recruitment induced by bicuculline might not only increase V_T and V_E but might also facilitate respiratory mechanics in vivo (Bosc, Resta, Walker, & Kanagy, 2010; Giordano, 2005; Haupt, Goodman Dm Fau - Sheldon, & Sheldon, 2012; Ninane, Rypens F Fau - Yernault, Yernault Jc Fau - De Troyer, & De Troyer, 1992; Ninane, Yernault Jc Fau - de Troyer, & de Troyer, 1993), reducing the work of breathing and, subsequently, causing a reduction in O₂ consumption. This, however, remains to be determined.

Rostral pFL bicuculline injections induces more prominent changes to all phases of the respiratory cycle

In this study, we introduced a novel, area-under-the-curve (AUC) analysis algorithm, aimed at dissecting the respiratory changes evoked within each phase of the respiratory cycle (late-E, Inspiration, post-I) as a result of the locale of injections and relative to baseline conditions for each of the airflow, DIA and ABD signals. This method revealed that the two most rostral injection sites elicited the most robust activation of ABD activity and induced the most significant alterations in expiratory airflow during both the late-E and post-I phases. Similarly, during the inspiratory phase, rostral injections generated the most pronounced changes in inspiratory airflow and DIA activity - particularly when positioned +0.6mm from the VIIc. These findings align with the observed increases in tidal volume depicted in Figure 4 H-I at those specific locations. Moreover, in the post-I phase, we also observed the most prominent negative airflow deflections at the most rostral injection sites, corresponding to concurrent ABD activation during that same timeframe.

By applying this novel analytical approach, which, to the best of our knowledge, has not been previously employed in research evaluating ABD responses triggered by pFL activation (Boutin et al., 2017; de Britto & Moraes, 2017; Huckstepp et al., 2015; Korsak et al., 2018; Magalhães et al., 2021; Pagliardini et al., 2011; Silva et al., 2019; Zoccal et al., 2018), we were able to reveal disparities in the AUC of ABD and airflow signals. These disparities, which went undetected with the analyses shown in Figure 4B to E, were particularly evident during the late-E and post-I periods. Furthermore, the division of the response into three distinct respiratory phases enabled a more precise assessment of how targeted stimulation affects a particular phase of the respiratory cycle, which could be used in future research to provide valuable insights into the underlying neural mechanisms of respiratory control.

Bicuculline responses following injections along the rostro-caudal axis of the ventral medulla can be differentiated by the time course and total extent of respiratory activation during late-E, inspiration and post-I phases.

Our results certainly indicate that bicuculline injections at various pFL locations induce a complex set of interactions between DIA and ABD muscle activation, as well as changes in airflow across the inspiratory and expiratory phases of respiration. Our final goal was to develop an algorithm which could accurately represent the complete “framework” of respiratory activation across each recorded signal with respect to each these phases. To this end, we introduced another novel analysis method that tracked the multivariate changes occurring throughout the respiratory cycle. We revealed a number of features of the respiratory cycle including late-E “tails”, inspiratory “bulbs” and post-I “feet” that allowed us to characterize alterations that are specifically relevant to active expiration. Given that this analysis is designed to maintain a consistent representation across conditions within experiments, these features can be catalogued into a “framework” of the breathing cycle during both baseline and post-injection conditions. Thus, this method is both sensitive to subtle changes in the activations of each component signal, but also more importantly, it is highly robust to the variability inherent in the acquisition of respiratory physiological measurements in vivo.

By quantifying the respiratory cycle in terms of 3D measures of distance we discriminated between the rostral and caudally evoked responses, as well as between the two most rostral locations which were previously undetected in the standard analyses. Our findings were consistent with the other analyses applied here with regards to the larger extent of respiratory loop distortions in the two most rostral positions in comparison to the more caudal injections. Furthermore, it was observed that injections in the most rostral group produced the highest amplitude changes to the late- E “tails”, as well as the longest late-E “tails” and post-I “feet”. Conversely, the most prominent inspiratory “bulbs” were observed in the +0.6 mm group, aligning with the most substantial changes in V_T observed in our single dimension analysis.

Collectively, these results imply that while the combined effects of bicuculline injection on the respiratory cycle were most pronounced at the two most rostral positions, there were slight variations in the effects exhibited during each respiratory phase at this level of the ventral medulla. These nuanced effects may suggest differences in the functional characteristics, phenotypes and projections of cell populations at these specific coordinates. It is worth noting that the pFL neurons responsible for active expiration are typically quiescent in adults at rest (de Britto & Moraes, 2017; Magalhães et al., 2021; Pagliardini et al., 2011), making their phenotypic and functional characterization challenging. Recent work started investigating electrophysiological characteristics of late-expiratory neurons (Magalhães et al., 2021). However, these observations were made in neurons situated closer to the tip of the facial nucleus in juvenile rats. It will be important to investigate now whether phenotypic and electrical properties of neurons at more rostral positions share the same identity.

pFL or RTN?

An ongoing debate centers on whether the neurons responsible for the generation of ABD recruitment are an independent oscillator composed of PHOX2B-negative, late-E firing neurons (pFL) (de Britto & Moraes, 2017; Magalhães et al., 2021; Pagliardini et al., 2011; Pisanski & Pagliardini, 2019) or part of the PHOX2B+/NMB+ chemosensitive neurons of the RTN/ventral paraFacial region (Abbott, Stornetta, Coates, & Guyenet, 2011; George, Ruth, Daniel, Stephen, & Patrice, 2020; Marina et al., 2010). While there is evidence that specific stimulation of PHOX2B+/NMB+ neurons causes ABD activity and active expiration (Abbott et al., 2011; George et al., 2020), these studies have not shown recordings from late-E PHOX2B/NMB+ neurons in their investigations. However, there is concrete evidence of the existence of PHOX2B-negative late-E firing neurons that are activated in response to bicuculline injections in the ventral medulla (de Britto & Moraes, 2017; Magalhães et al., 2021; Pagliardini et al., 2011). Because our anatomical results show very little PHOX2B/cFos+ staining at injections sites rostral to +0.1 mm, is fair to conclude that the ABD activity generated at that level is likely not the direct result of RTN activation or driven by PHOX2B+ neurons.

Conclusion

The results of this study highlight the functional diversity of the lateral parafacial region and emphasize the strong role of the most rostral locations in generating active expiration. The use of our novel multi dimensional map to assess physiological responses highlighted the functional nuances in respiratory responses to bicuculline along the pFL rostrocaudal axis, and advanced our understanding of the parafacial region.

Impact statement

Multivariate analysis of the respiratory cycle shows that GABAergic disinhibition in the rostral regions of the lateral parafacial area is linked to the most significant and enduring changes in active expiration.

Authors contribution

Concept and experimental design, A.P. and S.P.; Data collection, A.P.; Data analysis, A.P. and M.P.; Manuscript preparation, A.P, M.P., C.T.D and S.P.

Supplementary File 1: Rotating 3D respiratory trajectories highlight the features of breathing in baseline and response conditions. Trajectories from the sample recordings in Figure 6A (-0.2mm) across baseline and response time bins (colored as per color bar on right. Each dot represents a time point from a given time bin. A subsample of ∼11000 timepoints from each time bin have been randomly selected and are rotated about the DIA axis to highlight the architecture of respiratory trajectories. Dashed lines indicate the origin of each axis. Axis values expressed in units of standard deviation (SD) as per methods description.

Supplementary File 2:Rotating 3D respiratory trajectories highlight the features of breathing in baseline and response conditions. Trajectories from the sample recordings in Figure 6B (+0.1mm) across baseline and response time bins (colored as per color bar on right. Each dot represents a time point from a given time bin. A subsample of ∼11000 timepoints from each time bin have been randomly selected and are rotated about the DIA axis to highlight the architecture of respiratory trajectories. Dashed lines indicate the origin of each axis. Axis values expressed in units of standard deviation (SD) as per methods description.

Supplementary File 3:Rotating 3D respiratory trajectories highlight the features of breathing in baseline and response conditions. Trajectories from the sample recordings in Figure 6AC (+0.4mm) across baseline and response time bins (colored as per color bar on right. Each dot represents a time point from a given time bin. A subsample of ∼11000 timepoints from each time bin have been randomly selected and are rotated about the DIA axis to highlight the architecture of respiratory trajectories. Dashed lines indicate the origin of each axis. Axis values expressed in units of standard deviation (SD) as per methods description.

Supplementary File 4:Rotating 3D respiratory trajectories highlight the features of breathing in baseline and response conditions. Trajectories from the sample recordings in Figure 6D (+0.6mm) across baseline and response time bins (colored as per color bar on right. Each dot represents a time point from a given time bin. A subsample of ∼11000 timepoints from each time bin have been randomly selected and are rotated about the DIA axis to highlight the architecture of respiratory trajectories. Dashed lines indicate the origin of each axis. Axis values expressed in units of standard deviation (SD) as per methods description.

Supplementary File 5:Rotating 3D respiratory trajectories highlight the features of breathing in baseline and response conditions. Trajectories from the sample recordings in Figure 6E (+0.8mm) across baseline and response time bins (colored as per color bar on right. Each dot represents a time point from a given time bin. A subsample of ∼11000 timepoints from each time bin have been randomly selected and are rotated about the DIA axis to highlight the architecture of respiratory trajectories. Dashed lines indicate the origin of each axis. Axis values expressed in units of standard deviation (SD) as per methods description.

Acknowledgements

This work was supported by CIHR Project scheme grant to SP (388717) and an NSERC Discovery operating grant to CTD (2021-02926). AP was supported by an NSERC CGS doctoral scholarship (CGSD3 – 535109 – 2019).