Introduction

Neural circuits that regulate behavior must integrate rapid sensorimotor responses with slower fluctuations in internal physiological state. Reflexive reactions unfold within milliseconds, whereas sensory gain and behavioral sensitivity can vary over seconds to minutes depending on arousal, context, and physiological state. In the nociceptive system, this coordination across timescales is particularly important because descending brainstem pathways dynamically regulate pain sensitivity and are major targets of pharmacological and neuromodulatory interventions. However, despite the central role of these circuits in controlling pain and mediating therapeutic effects, how neural activity within descending pain-control pathways is organized across multiple timescales remains poorly understood. Understanding these temporal dynamics is therefore important not only for explaining how pharmacological and neuromodulatory interventions influence descending pain control, but also for guiding the design and optimization of such therapies.

A central component of this descending pain-control system is the rostral ventromedial medulla (RVM), a key node in a modulatory pathway that regulates spinal nociceptive processing ^1,2. Through direct projections to the dorsal horn, the RVM exerts potent inhibitory and facilitatory control over no-ciceptive transmission ^3–9, dynamically adjusting pain sensitivity according to behavioral and physiological context. This descending modulation is mediated by two principal classes of RVM neurons: ON- and OFF-cells. Increased ON-cell activity facilitates nociceptive transmission and promotes hyperalgesia, whereas activation of OFF-cells suppresses no-ciceptive signaling and produces hypoalgesia. During noxious stimulation, ON-cells increase firing immediately prior to withdrawal responses, while OFF-cells pause, providing a physiological signature for identifying both populations ¹⁰ (Figure 1). Cells that do not exhibit stimulus-linked modulation are typically classified as NEUTRAL-cells. Together, these populations contribute to a descending control system that dynamically regulates nociceptive gain. Experimental activation of OFF- and ON-cells respectively results in hypoalgesia and hyperalgesia ^{6–9,11–14}. Furthermore, inactivation of ON-cells attenuates hypersensitivity in persistent pain models ^15–20. More generally, periods of ON-cell dominance have been shown to be associated with relative hypersensitivity compared to periods of OFF-cell dominance ²¹. Therefore, RVM function cannot be understood solely through reflex-linked ON- and OFF-cell responses; ongoing fluctuations in these populations are also likely to shape descending pain modulation.

Figure 1.

Consistent with this possibility, several features of RVM organization suggest that its function extends beyond reflex-linked modulation. The RVM integrates bottom-up sensory input from the spinal cord with sustained top-down influences from the periaqueductal gray, hypothalamus, and other homeostatic centers ^22–24. Studies comparing direct and indirect recruitment of spinal and trigeminal inputs through the parabrachial complex further show that RVM neurons can exhibit both transient and prolonged responses depending on the route of activation, suggesting that they integrate nociceptive signals across distinct temporal scales ²⁵. Early electrophysiological recordings demonstrated that ON- and OFF-cells show broadly anticorrelated fluctuations even in the absence of sensory stimulation ²⁶. However, these obser-vations were limited in scope, and the temporal structure of both evoked and ongoing RVM activity has never been quantitatively characterized.

We hypothesized that RVM neurons operate across multiple temporal scales, integrating fast responses associated with reflex-linked control with slower fluctuations reflecting on-going network or state-dependent modulation. To test this idea, we combined neuronal recordings from identified ON-, OFF-, and NEUTRAL-cells with probabilistic modeling to characterize both stimulus-evoked and ongoing activity in the RVM using a robust lightly anesthetized rodent model in which activation of OFF-cells has been demonstrated to produce hypoalgesia and activation of ON-cells to produce hyperalgesia ^11,12,15,27. In this well characterized regime, animals (and neurons) remain responsive to heat stimuli applied to the hindpaw, whilst preventing spontaneous movement ^10,28. Our analyses reveal multi-timescale population dynamics in which fast, multi-phase responses to noxious stimuli unfold over seconds to tens of seconds, while slower quasi-periodic fluctuations in ongoing activity occur over minutes. These oscillations, specific to ON- and OFF-cells, were predictable from past activity. Together, these results demonstrate that ongoing multi-timescale dynamics shape activity in descending pain-modulatory neurons, providing a dynamical framework for how descending pain modulation is organized across time.

Results

Canonical ON/OFF responses to noxious stimulation

To understand how descending pain-modulatory neurons operate across multiple timescales, we first examined the activity of identified RVM neurons during noxious stimulation, consistent with previous work ^{10,26,29–31}. We used a well-established setup of light anesthesia and extracellular single unit recordings (see Methods), which preserves reliable with-drawal responses to noxious heat while minimizing spontaneous movement and allowing stable recordings from identified RVM neurons across repeated trials. Recordings from ON- and OFF-cells revealed structured population responses that extended well beyond the brief epoch surrounding with-drawal behavior (Figure 1). These neurons exhibited multi-phase activity patterns comprising an initial rapid stimulus-linked change in firing followed by prolonged recovery dynamics lasting tens of seconds.

Figure 1b shows representative recordings from an ON-cell and OFF-cell recorded simultaneously in one animal, along-side two NEUTRAL-cells recorded in a second animal. During noxious heat stimulation (vertical dashed line), the ON-cell increased firing in close temporal proximity to paw with-drawal, whereas the OFF-cell exhibited a sharp pause in firing, consistent with canonical reciprocal dynamics. In contrast, NEUTRAL-cells showed no stimulus-associated changes in firing (Supplementary Figure S3), consistent with their canonical lack of nociceptive responsiveness. Notably, spontaneous switches between OFF- and ON-dominant activity were also observed during inter-trial intervals (e.g., at 180 s and 520 s), indicating that the activity of these cell classes is not solely stimulus-evoked but also reflects ongoing network dynamics.

Population responses are summarized in Figure 1c, which shows peri-stimulus time histograms (PSTHs) for 70 ON-cells and 60 OFF-cells, aligned to the time of paw withdrawal (Figure 1c, upper) or to the onset of heat stimulation (Figure 1c, lower). In total, 185 ON-cell and 155 OFF-cell trials were aligned to withdrawal, while 182 ON-cell and 151 OFF-cell trials were aligned to heat onset the models used more, because they only used intervals up to 100 s, here we are plotting longer intervals to 150 s(2 to 3 trials/cell, trial breakdown per animal is provided in the Supplementary Information). PSTHs showed a clear lag and shallower slopes when aligned to heat onset compared to withdrawal, confirming that our dataset captures the canonical, reciprocal ON/OFF responses associated with nocifensive reflexes.

Notably, ON- and OFF-cell activity also evolved on longer timescales that extended beyond the time of withdrawal, suggesting both fast and slow dynamics. We therefore modeled ON- and OFF-cell population dynamics using Bayesian regression to characterize these response and recovery timescales.

Evoked ON/OFF activity unfolds over multiple timescales

To characterize the evoked response and recovery dynamics of ON- and OFF-cell populations, we applied a piecewise non-linear Bayesian regression model with a Poisson likelihood to the summed spike counts per time bin across the entire ON- and OFF-cell populations, aligned initially by with-drawal time (Figure 2a; see Methods). We chose a window from 10 seconds before, to 100 seconds post withdrawal (or heat onset, discussed later). This resulted in 195 ON-cell and 170 OFF-cell trials aligned to withdrawal, and 197 ON-cell and 172 OFF-cell trials aligned to heat onset. Cells were not distinguished by their firing rate pre-trial and all trials with a withdrawal were included. We focused on population-level structure rather than animal- or cell-specific effects, given the small number of trials per animal and the absence of a principled basis for differential weighting. Under the assumption of exchangeability, pooling across cells and animals reduces uncorrelated noise and preserves shared population dynamics. We used Bayesian highest density intervals (HDIs), defined as the minimum width Bayesian credible intervals ³². A 97% HDI means there is a 97% probability that the unobserved parameter lies within the given interval, given the model and data.

Figure 2.

For all analyses of evoked activity, we decomposed each trial into two segments: a response phase, defined as the monotonic departure of firing rate from the pre-withdrawal baseline (increase for ON-cells, decrease for OFF-cells), and a recovery phase, defined as the monotonic return toward a post-withdrawal baseline. These were separated by a single transition timepoint t_s, corresponding to the extremum (peak for ON-cells, nadir for OFF-cells) before recovery began. A null model assuming a constant mean spike count performed significantly worse than all piecewise models, confirming that both an evoked departure and a distinct recovery process were required to explain the data (see Supplementary Figure S8). Since all trials were aligned to withdrawal, variability in the withdrawal latency (time between heat onset and withdrawal) cannot shift post-withdrawal activity into the pre-withdrawal window, making an averaging artifact an unlikely source of the observed recovery dynamics.

The fitted switching timepoint between response and recov-ery differed for ONand OFF-cells (Figure 2g): 0.18 s before withdrawal (97% HDI: 0.27–0.095 s before), and s after withdrawal (97% HDI: 0.85–1.9 s after). Thus ON-cells peak slightly before withdrawal, while OFF-cells enter a brief post-withdrawal depression before beginning their recovery phase. This relative timing is visible in previous literature but has not been explicitly modeled ³⁰. The response phase was modeled as an exponential rise for ON-cells and a sigmoidal decrease for OFF-cells. We also considered a Gaussian cumulative density function (er f) decrease for OFF-cells, which has lighter tails than the sigmoid, and a different biological interpretation (see Methods, Discussion); this fitted significantly worse than the sigmoid across all recovery models (see Supplemental Gaussian CDF Fits, Tables S9 and S10). The recovery phase was parameterized using alternative functional forms: single vs. double exponential recovery for ON-cells, and double exponential vs. linear + exponential recovery for OFF-cells (Figure 2b).

As shown in Figure 2c,g, ON-cell recovery was best captured by a double exponential model, consisting of a superposition of two distinct decay rates (timescales). ON-cell firing rose on a sub-second timescale before paw withdrawal, with time constant τ_pre = 0.653 s (97% HDI: 0.570–0.739 s). After paw withdrawal, the recovery exhibited both a fast component τ_fast = 4.157 s (97% HDI: 3.582–4.722 s) and a slow component τ_slow = 164.223 s (97% HDI: 129.571–199.809 s). A single exponential model (Figure 2e) could not simultaneously capture the sharp evoked peak and the prolonged decay, instead over- and undershooting different portions of the trajectory.

For OFF-cells, recovery was best described by the combined linear and exponential model (Figure 2d,h). OFF-cell firing declined sharply before withdrawal, with a sigmoidal response slope of k = 2.852 s⁻¹ (97% HDI: 2.333–3.392 s⁻¹). The post-withdrawal recovery contained a fast exponential component τ = 4.760 s (97% HDI: 3.711–5.861 s), together with a slower linear drift of k_linear = 0.186 counts per second (97% HDI: 0.174–0.196 counts per second). Although the double exponential alternative fit the early part of the recovery (Figure 2f), it failed to asymptote to baseline over the recorded interval, making the linear+exponential model the better descriptor of the OFF-cell dynamics we observed.

We next quantified the timing of the response and recovery phases of the evoked changes in firing rate in terms of percentage deviations from baseline (see Supplemental Figure S1 for full posteriors). OFF-cell population firing had already decreased by 10% from baseline at 1.5 s before the paw withdrawal (97% HDI: 1.746–1.306 s), and ON-cell firing increased by 10% on a similar timescale, at 1.7 s before withdrawal (97% HDI: 1.9–1.5 s), indicating that both populations depart from baseline well before the behavioral response. The two classes also crossed 50% of their total change within a narrow pre-withdrawal window (OFF: 726 ms before; 97% HDI: 810–645 ms; ON: 630 ms before; 97% HDI: 690–570 ms), showing that the dominant portion of the evoked response is temporally locked to the upcoming reflex. Notably, ON-cell firing reached 90% of its peak well before paw withdrawal (240 ms before, 97% HDI: 320–180 ms), while OFF-cell suppression did not reach 90% until shortly after the behavioral response (41 ms post, 97% HDI: 95 ms pre –180 ms post), reflecting a temporal asymmetry in how the two populations complete their evoked transitions.

Comparison of these values to the transition time for OFF-cells described above shows that OFF-cell firing continues to decrease for several seconds post-withdrawal before beginning to recover. At the population level, ON-cell activation both precedes withdrawal and completes its evoked component rapidly, whereas OFF-cell suppression is slower to reach its minimum and persists into the post-withdrawal period. Consistent with this ordering, ON-cell activation reaches its peak before OFF-cell firing reaches full suppression, indicating that the shift toward an ON-dominant state begins before the OFF-dominant state has fully terminated. Because this timing was observed in a pseudo-population rather than within simultaneously recorded ON/OFF pairs, it should be interpreted as a population-level temporal organization rather than a claim about pairwise interactions.

We used separate baselines for the pre- and post-withdrawal activity, which allowed us to quantify recovery on the same percentage scale, accounting for any slow modulation of tonic firing relative to the trial length. ON-cells returned to 50% of their post-withdrawal baseline by 15.39 s after the withdrawal (97% HDI: 9.373-22.43 s), but required more than an order of magnitude longer to approach full recovery, reaching 90% of their baseline by 222 s (97% HDI: 183.6-248.8 s). OFF-cells increased back to 50% of their baseline by 18 s post-withdrawal (97% HDI: 15-21 s) but to 90% of their base-line by 83 s (97% HDI: 82-84 s). These timescales are considerably longer than the typical duration of EMG activity ³⁰, indicating that the prolonged temporal structure cannot be explained as a trivial consequence of the motor output itself.

An important consequence of modeling ON- and OFF-cell activity parametrically is that it allows the sharpness of the evoked transitions to be quantified directly, using the slope or time constant of the transitions, rather than inferred qualitatively from peri-stimulus averages (as in Figure 1c) ³³. For OFF-cells, the slope of the sigmoidal fits to firing-rate decreases was significantly steeper for the withdrawal-aligned data (2.852 total spikes per bin, 97% HDI: 2.333–3.392 spikes per bin), compared to the heat-aligned data (1.747 total spikes per bin, 97% HDI: 1.476–2.041 spikes per bin) (Supplementary Tables S6 and S7). Similarly, ON-cells exhibited faster increases in firing when aligned to withdrawal, with exponential time constants of 0.653 s (97% HDI: 0.570–0.739 s), compared to 1.278 s (97% HDI: 1.129–1.426 s) when aligned to heat onset (Supplementary Tables S3 and S5). In contrast, NEUTRAL-cells (n = 32; 69 withdrawal-aligned and 70 heat-aligned trials) exhibited no discernible modulation associated with either heat onset or withdrawal behavior (Supplementary Figure S3). Importantly, these neurons were recorded under the same anesthetic and physiological conditions as the ON and OFF populations, supporting the cell-type specificity of stimulus-linked modulation.

Taken together, these results shed new light on the traditional reflex-locked view and demonstrate that ON- and OFF-cell dynamics unfold on multiple, cell-type–specific timescales that extend well beyond paw withdrawal. This led us to examine whether ON- and OFF-cells exhibit ongoing rhythmic activity in the absence of stimulation, and how such slow oscillations might interact with stimulus-evoked responses.

RVM ON- and OFF-cells exhibit slow quasi-periodic structure during ongoing activity

To investigate the ongoing activity of RVM neurons, we analyzed ON-, OFF-, and NEUTRAL-cells (N_ON = 25, N_OFF = 25, N_NEUTRAL = 32) recorded in the absence of external stimulation (Figure 3a).

Figure 3.

Even without sensory input, ON- and OFF-cells exhibited slow, quasi-periodic fluctuations in firing rates (examples in Figure 3a). To quantify these apparent rhythms, we computed power spectral density (PSD) estimates for each cell class (Figure 3b,c,d). Both ON- and OFF-cells exhibited clear low-frequency spectral peaks, often accompanied by harmonics, indicating quasi-periodic dynamics with a characteristic period of approximately five minutes s, full width at half maximum (FWHM) = 16.92 s;. NEUTRAL-cells, in contrast, lacked consistent spectral peaks, instead showing broad low-frequency components and prominent high-frequency power (> 10 Hz), reflecting regular, tonic firing. These results reveal a slow oscillating modulation seemingly specific to ON- and OFF-cells.

We next summarized firing variability over the full recording period (Figure 3e). NEUTRAL-cells fired regularly and rapidly (Coefficient of Variation = 0.35 ±0.26; mean Inter-Spike Interval = 0.17± 0.25 s, ~10 Hz). In contrast, ON- and OFF-cells exhibited highly variable, burst-like firing patterns, with coefficients of variation well above 1 (CV_ON = 4.91± 4.32, CV_OFF = 3.70 ±3.09) and broad interspike interval (ISI) distributions (ISI_ON = 1.25± 4.80 s, ISI_OFF = 0.84± 1.36 s; Supplementary Table S1). Because relatively few neurons were recorded per animal, we did not attempt further sub-clustering.

Together, these analyses indicate that ON- and OFF-cells exhibit structured, quasi-periodic fluctuations in ongoing activity, whereas NEUTRAL-cells maintain high-rate, non-rhythmic firing. However, frequency-domain measures such as PSD and summary statistics capture only average power and variability; they discard temporal phase information and cannot resolve how firing unfolds over time.

To address this, we next used Gaussian process (GP) modeling to characterize ongoing dynamics directly in the time domain. We implemented two complementary GP approaches: a purely periodic GP with a Gaussian likelihood to identify dominant rhythmic timescales, and a generalized Poisson GP with a SoftPlus link function to test whether these slow fluctuations were structured and predictable on a per-neuron basis.

A periodic Gaussian process reveals dominant slow oscillations in ON- and OFF-cell activity

To identify the characteristic timescales underlying the slow fluctuations observed during ongoing activity, we applied a periodic GP model with a Gaussian likelihood to the smoothed firing rates of each neuron, scanning across 200 candidate periods (0.1–1500 s). This model quantifies how strongly each frequency explains the observed firing pattern, assuming only that the waveform is smooth and periodic (Figure 4). ON- and OFF-cells showed clear minima in the negative log marginal likelihood (NLML) near 300 s, with additional subharmonic peaks at 600 and 900 s, indicating quasi-periodic dynamics with a dominant 5-minute cycle. In contrast, NEUTRAL-cells (with one exception) displayed flat NLML profiles with no consistent minima, without evidence of stable rhythmicity. These results establish that slow oscillations are a robust and reproducible feature of ON- and OFF-cell activity.

Figure 4.

A Poisson Gaussian process predicts slow quasi-periodic fluctuations in ongoing RVM firing

We next asked whether these slow oscillations were predictable—that is, whether past firing reliably forecasted future activity. To test this, we used a Poisson GP model with a SoftPlus link function, allowing the variance of firing rates to scale with their mean through a dispersion parameter (α). Each model was trained on all but the final 300 s of data and tasked with predicting the held-out test segment (Figure 5a). The fitted models (solid lines) accurately captured slow fluctuations during training and generalized into the test window (dotted lines) robustly for ON-cells, variably for OFF-cells, and not for NEUTRAL-cells.

Figure 5.

Training and test accuracies, quantified by pseudo-R², were computed for each cell (Figure 5b). We performed two-sided T-tests to assess significance (R² ≠0) using sign flipping permutations of the R² value of all cells per animal, and the median R² at each permutation (see Supplemental Figure S4 for distributions). This accounted for the potential correlation of ON- and OFF-cells recorded from the same animal. All three cell classes achieved significant training accuracy (p^ON = 0.0047, p^OFF = 0.0004, p^NEUTRAL < 0.0001). Comparing R² values using a similarly permuted and Bonferroni corrected Mann Whitney U test revealed no significant differences between training R² for the three groups (see Supplemental Figure S6). This was expected, since even apparently aperiodic firing rates can contain some variance which can be captured by the GP. Therefore, we turned to the predictive accuracy to confirm the periodicity of each cell class, based on the predictions given by the periodic kernel.

By contrast, test predictability, quantified by pseudo-R², dissociated the three populations (Figure 5b). ON-cells showed significant generalization to unseen data (p = 0.0067), with the large majority of cells exceeding baseline performance (i.e., >0). OFF-cells, as a population, did not show reliable predictive structure (p = 0.3127), reflecting a mixture of well-predicted and poorly predicted neurons. NEUTRAL-cells performed significantly below baseline (p = 0.0033), demonstrating a lack of periodicity. Group comparisons revealed significant differences between NEUTRAL-cells and both ON- and OFF-cells (p < 0.001 and p = 0.006, respectively), while ON- and OFF-cells did not differ significantly at test (p = 0.939). Kruskal–Wallis tests confirmed significant group differences at test time (p < 0.001) but not during training (p = 0.117). Thus, reproducible long-timescale predictability was most evident in ON-cells, while OFF-cells exhibit oscillatory structure that is present but less consistent across neurons.

Posterior hyperparameters (Figure 5c; Table 1) revealed further distinctions among cell classes. Estimated oscillation periods for both ON- and OFF-cells clustered near 300 s, corroborating the frequency-domain analyses. Although NEUTRAL-cells also displayed nominal period clustering near 300 s, their low latent variance and poor predictive performance indicate that this reflects a weak parameter constraint rather than true oscillatory structure.

Table 1.

ON- and OFF-cells exhibited positive dispersion parameters (α > 0; p_ON = 0.0011, p_OFF = 0.0105), consistent with bursty, over-dispersed firing. In contrast, NEUTRAL-cells showed negative dispersion (p_NEUTRAL < 0.001), indicating highly regular firing with variance below Poisson expectation. Latent GP variance further separated the classes: ON-cells showed the largest modulation amplitude, OFF-cells intermediate variance, and NEUTRAL-cells near-flat latent structure.

Finally, single-cycle phase-normalized reconstructions (Figure 5d) revealed cell-specific oscillatory signatures in ON- and OFF-cells, whereas NEUTRAL-cells showed near-constant phase profiles, consistent with an absence of structured modulation.

Together, these results demonstrate that slow, quasi-periodic fluctuations in RVM activity are most consistently expressed and predictably organized in ON-cells, present but heterogeneous in OFF-cells, and absent in NEUTRAL-cells. This structure indicates that ongoing RVM dynamics are not fully captured by stochastic firing variability, but include structured components operating on multi-minute timescales.

Low-frequency coherence between RVM activity and heart rate

The RVM is implicated not only in nociceptive modulation but also in autonomic control, influencing cardiovascular and respiratory outputs ^34–39. Given this dual role, we next asked whether the slow RVM oscillations identified above might couple to an autonomic variable—specifically, heart rate, a key homeostatic parameter that itself exhibits rhythmic fluctuations.

Previous studies using small samples of ON- and OFF-cells reported correlations between ON/OFF-cell firing and cardiovascular measures such as heart rate and blood pressure ^40,41. Here, we extended these findings by testing whether slow RVM oscillations and heart rate exhibit a consistent phase relationship using coherence analysis, a frequency-domain measure that detects shared periodic structure even in the presence of noise or constant phase differences.

The neuronal firing rates and simultaneously recorded heart rates (Figure 6a) were segmented and demeaned, and coherence estimates were computed for each cell (Figure 6b). A subset of cells (5/32 NEUTRAL, 4/25 OFF, 8/25 ON) showed significant low-frequency coherence with heart rate (Figure 6c). Although the peak coherence frequencies did not exactly match the dominant ON- and OFF-cell oscillation period (~ 300 s), spectral peaks in heart rate were observed at multiples of this period (Figure 6d), suggesting that RVM activity and cardiac rhythms share a common underlying regulatory timescale.

Figure 6.

These results provide preliminary evidence for partial coupling between descending pain-modulating activity and autonomic rhythms, consistent with the integrative role of the RVM in coordinating sensory and homeostatic processes.

Discussion

Pain-regulatory circuits must coordinate rapid defensive responses with slower fluctuations in internal physiological state. Our results show that neurons in the rostral ventro-medial medulla (RVM), a central hub of descending pain control, operate across distinct temporal scales. By combining neuronal recordings with probabilistic modeling, we find that identified ON- and OFF-cells exhibit both fast reflex-linked responses and slower quasi-periodic fluctuations in ongoing activity. Importantly, these slow oscillatory features were not observed in NEUTRAL-cells recorded under the same anesthetic and physiological conditions, arguing against a purely global effect of anesthesia, recording conditions, or nonspecific physiological drift. The coexistence of these dynamics within the ON- and OFF-neurons indicates that descending pain modulatory activity spans both rapid stimulus-linked responses and slower quasi-periodic fluctuations, potentially contributing to context-dependent regulation of nociceptive processing. Our findings should therefore be interpreted primarily as evidence that descending pain-modulatory activity exhibits structured multi-timescale organization, although it is not definitive evidence for a specific mechanistic origin or functional role of these dynamics.

Multi-timescale neural dynamics are well documented in cortical and neuromodulatory systems, where they are often attributed to large recurrent networks and global brain-state fluctuations ^42,43. Our findings extend this principle to a brainstem circuit that directly regulates nociceptive transmission. In this framework, ON- and OFF-cells are not simple reciprocal antagonists. Rather, they contribute distinct and temporally staggered components of descending modulation—consistent with their different pharmacology and the fact that each population can be manipulated independently ¹. Given the limited local connectivity among ON- and OFF-cells ⁴⁴, the observed dynamics are unlikely to arise solely from local microcircuits and may instead reflect broader feedback loops linking the RVM with spinal, midbrain, and hypothalamic control systems.

Bayesian modeling of withdrawal-aligned activity revealed that ON- and OFF-cells exhibit structured, multi-phase responses rather than simple bursts or pauses. Population activity was more tightly aligned to the withdrawal reflex than to the onset of noxious stimulation, confirming classical observations ¹⁰ while providing a quantitative description of the temporal structure of these responses. Importantly, the extrema of the two cell classes were temporally offset: ON-cell activity peaked 180 ms before the withdrawal, whereas OFF-cell suppression reached its nadir >2 seconds after. Thus, although the two populations are broadly anticorrelated, their dynamics are not simply mirrored but instead unfold in a temporally staggered manner.

We investigated both a sigmoidal and a Gaussian CDF function for the initial response of OFF-cells. Since we fitted our models to the summed population activity over animals and trials, we expect the slope of response function to be indicative of the general mechanism underlying the response. With this in mind, a sigmoidal decrease satisfies the logistic differential equation, in which there is recurrent inhibition (see Supplemental Gaussian CDF Fits). In contrast, a Gaussian CDF response suggests that OFF-cells cease firing independently around some response time t₀. These two distributions are very difficult to distinguish between given the small timescale of the initial response, and the differences between the distributions are most pronounced in their tails, making the visual fits to data appear very similar. We observed that the sigmoidal function fit better to the response in terms of the ELPD; this is likely due to the heavier tails of the sigmoidal function fitting the noisy data better, and may change with a larger sample size or hierarchical model.

Analysis of the recovery phase revealed multiple timescales in both populations. ON-cell recovery was best captured by a double-exponential function, whereas OFF-cell recovery was best described by a combined linear and exponential trajectory. Neither population returned to its pre-stimulus firing rate within 100 s. This separation of fast and slow components suggests distinct biological processes underlying the initial response and the prolonged return toward baseline. The abrupt OFF-cell pause and sharp ON-cell peak are consistent with positive-feedback mechanisms that support rapid, precisely timed actions, whereas the slower recovery component aligns with the prolonged adjustments in nociceptive responsiveness observed behaviorally. Consistent with this, blocking the OFF-cell pause suppresses withdrawal ⁴⁵, and attenuating the ON-cell burst reduces reflex magnitude ^13,46. Moreover, the slow post-withdrawal drift toward baseline correlates with increased sensitivity to subsequent stimuli ⁴⁷ and with altered ON/OFF-cell balance in persistent pain models (e.g., spinal nerve ligation, CFA inflammation ^48,49). Together, these findings suggest that the initial response and its hysteresis-like continuation serve complementary roles: coordinating the reflex itself while shaping a longer-lasting nociceptive state in the face of potential threat.

Beyond event-related responses, we identified slow, quasi-periodic oscillations in ongoing firing, with a dominant period near five minutes. These oscillations, evident in both the power spectrum and Gaussian process models, were consistent across animals and cell classes but not evident in NEUTRAL-cells, suggesting that they reflect cell-classspecific structure within the pain-modulatory network rather than a global effect (e.g., anesthesia). The presence of statistically predictable low-frequency structure implies that the RVM alternates between facilitative (ON-dominant) and inhibitory (OFF-dominant) modes even in the absence of noxious input. One interpretation of these dynamics is that ON- and OFF-cell populations track slowly varying latent physiological states that influence nociceptive gain over time. Although we did not observe phase-dependent effects on withdrawal latency (Supplemental Figure S2), the experiment was not designed to test this directly. Prior work demonstrating increased nociceptive sensitivity during ON-cell dominance ^3,21 suggests that such internal dynamics would bias behavioral responsiveness under natural conditions.

An important unresolved question is whether the slow fluctuations observed here primarily reflect intrinsic dynamics within descending pain-modulatory circuits or coupling to broader physiological state variables such as autonomic, thermoregulatory, arousal, or homeostatic processes. These possibilities are not mutually exclusive. The absence of comparable low-frequency structure in NEUTRAL-cells argues against a purely global anesthetic or physiological origin, but the observed coherence between subsets of RVM neurons and heart-rate fluctuations suggests that at least part of the slow structure may reflect coordinated systems-level state changes.

The lightly anesthetized preparation used here has been extensively employed to characterize canonical ON- and OFF-cell physiology because it preserves reliable nocifensive reflexes while permitting stable long-duration recordings from identified neurons. Nevertheless, anesthesia likely alters global brain-state dynamics and may influence the amplitude, frequency, or coherence of the slow fluctuations described here. Although the absence of comparable low-frequency structure in NEUTRAL-cells argues against a purely nonspecific anesthetic effect, determining which aspects of these dynamics generalize to awake conditions remains an important direction for future work.

Given the limited local connectivity among ONand OFF- cells ⁵⁰, and the lack of evidence for direct mono-synaptic connectivity between the populations ⁴⁴, the low-frequency structure is unlikely to arise from intrinsic microcircuit loops and more plausibly reflects intrinsic membrane properties of these cells or broader feedback interactions involving the spinal cord, periaqueductal gray or other brainstem regions, hypothalamus, or forebrain ^{3,6,18,22,51,52}. The strong synchrony within each functional population, combined with pronounced asymmetries between them, suggests coordination via shared modulatory inputs rather than reciprocal microcircuitry.

Because recordings were obtained from small pseudo-populations rather than large simultaneously recorded ensembles, we cannot determine the extent to which these slow fluctuations are synchronized across the broader ON- and OFF-cell populations. Nevertheless, the consistency of the low-frequency structure across neurons recorded across multiple experiments suggests that these dynamics are not confined to isolated cells. Future approaches using large-scale electrophysiology or calcium imaging could help resolve the population-level organization of these slow dynamics.

Furthermore, there is substantial molecular heterogeneity within ON- and OFF-cell populations ^{6–9,53–55}. Because neurons were classified physiologically rather than molecularly, we cannot determine how the temporal dynamics described here relate to any identified molecular subclasses of ON- and OFF-cells. However, despite preventing a straightforward molecular classification of ON- and OFF-cells, such diversity may support the flexible functional organization observed here, enabling the integration of upstream control signals based on both internal (slow periodicity) and external (stimulus-evoked) variables. The RVM receives convergent non-nociceptive inputs ⁵⁶ and participates in thermoregula-tory and autonomic processes ^41,57, further supporting a role as a multimodal hub integrating information across modalities and timescales.

Pain modulation is closely intertwined with autonomic and behavioral adjustments that support homeostasis ^23,58,59. Prior work reports associations between ON/OFF-cell firing and heart rate or body temperature ^40,41. Motivated by this, we examined heart rate as a potential covarying signal. Because both heart rate and RVM activity are shaped by many independent factors which tend to decorrelate signals, simple correlation is difficult to interpret; coherence provides a more targeted measure of shared structure at the timescale of interest (see Supplemental Coherence vs correlation). A subset of neurons across all three cell classes showed significant low-frequency coherence with heart rate, and heart-rate spectra exhibited peaks at the RVM oscillation frequency and its harmonics. We interpret these findings cautiously: they indicate that the oscillation is detectable outside the RVM and therefore physiologically meaningful but do not imply direct causation. Notably, coherence in some NEUTRAL-cells—despite their lack of intrinsic rhythmicity—suggests that modulatory inputs influencing the slow oscillation distribute across RVM cell types and underscore the still-unclear functional role of NEUTRAL-cells.

In conclusion, our results demonstrate that descending pain-modulatory neurons operate across multiple timescales, from seconds-long reflex-related responses to minute-scale endogenous rhythms. This demonstrates that descending pain-modulatory activity contains structured dynamics even in the absence of noxious stimulation and places the RVM within a broader feedback loop balancing sensory processing, defensive behaviors, and internal physiological states. Understanding this temporal organization provides a quantitative framework for investigating how brainstem circuits coordinate defensive behavior and offers new insight into how pain-regulatory systems can be modulated.

Methods

Experiments

All experiments followed the guidelines of the National Institutes of Health and the Committee for Research and Ethical Issues of the International Association for the Study of Pain, and were approved by the Institutional Animal Care and Use Committee at Oregon Health & Science University (OHSU).

Male and female (n = 77 and 6 respectively) Sprague-Dawley rats (Charles River; 191-348 g, mean ± SD: 264± 27 g) were deeply anesthetized using isoflurane (4-5%) and a catheter placed in the external jugular vein for subsequent infusion of the short-acting barbiturate methohexital. They were then transferred to a stereotactic frame and a small craniotomy made to gain access to the RVM. Body temperature was monitored and maintained at 36-37 °C with a heating pad, and heart rate was monitored using EKG. Mean ± SD heart rate was 391.1 ± 34.4 bpm, with mean range over the protocol 0.62 ± 12.0 bpm across all animals.

After preparatory surgery was complete, the anesthetic plane was lowered such that a heat stimulus applied at 5-min intervals to the plantar surface of the hindpaw (radiant or contact heat delivered using a Peltier device) elicited a paw withdrawal. The stimulus consisted of a ramp from 36 to 53 °C over a period of 11 s. The stimulus was terminated when the animal withdrew the paw (mean ± SD paw withdrawal temperature: 49.0± 2.0 °C, latency: 8.7 ± 1.1 s). The paw surface was maintained at 35 °C between trials. Recordings were performed between August 2020 and August 2024.

Once the animal responded to the heat stimulus, methohexital was infused at a rate (37.8–77.5 mg/kg/hr, mean ± SD: 60.0 ± 7.6 mg/kg/hr) that maintained this lightly anesthetized state (as indicated by stable paw withdrawal response) while preventing spontaneous movement. The lightly anesthetized preparation was necessary to maintain recording stability over prolonged periods while preserving stimulus-evoked withdrawal responses ^28,60. The experiment was performed in low ambient light conditions (< 10 lux). Extracellular single-unit recordings were acquired using a stainless-steel microelectrode (Frederick Haer & Co). Signals were amplified (10k) and band-pass filtered (400 Hz to 15 kHz, Neurolog, Digitimer) before analog-to-digital conversion at 32k samples/s for real-time spike detection and monitoring using Spike2 software (CED, Cambridge, UK). Correct waveform identification was verified on an individual spike basis at the conclusion of the experiment using Spike2 template matching and cluster analysis. EMG activity (to monitor withdrawal reflexes) was also recorded using Spike2. An RVM neuron or neurons was isolated and characterized as an ON-, OFF-, or NEUTRAL-cell based on changes in firing rate associated with withdrawal of the paw from a noxious mechanical stimulus applied to the hind-paw ^26,61. ON-cells are defined by a period of activity beginning just prior to withdrawal from a noxious stimulus (or maintained firing if already active at the time of stimulus delivery). OFF-cells stop firing just prior to withdrawal (or remain silent if inactive). NEUTRAL-cells do not respond. We performed/analyzed data from two complementary protocols, each designed to capture different aspects of RVM activity during noxious stimulation or extended unstimulated periods. In the “evoked” protocol, noxious heat stimuli were delivered at 5 min intervals over a period of 16–18 min. A total of 104 cells were recorded in 53 animals: 35 OFF-cells, 45 ON-cells, and 24 NEUTRAL-cells. The “ongoing” protocol was designed specifically to examine ongoing activity in the absence of noxious stimulation. Cell activity was monitored for a period of 30 min following an initial noxious heat trial to classify the cell under study. Evoked responses were confirmed in at least two trials at the end of the experiment. Using this protocol, we recorded 62 cells in 30 animals: 25 OFF-cells, 25 ON-cells, and 12 NEUTRAL-cells. Together these protocols provided both noxious-evoked and prolonged recordings of ongoing, or “spontaneous” activity of ON-, OFF-, and NEUTRAL-cells, enabling systematic analysis of stimulus-driven dynamics and slow ongoing fluctuations. A detailed breakdown of cell numbers and animals used in each analysis is provided in Supplemental Table S11.

Recording sites in RVM were marked with an electrolytic lesion at the conclusion of the experiment. Animals were over-dosed with methohexital, and perfused transcardially with saline followed by 10% formalin. Brains were removed and post-fixed for 24 h in 10% formalin, then equilibrated for 24–72 h in 30% sucrose in PBS at 4 °C. Brains were sectioned at 40 to 60 µm, and recording sites were plotted (Supplementary Figure S11). The RVM was defined as the nucleus raphe magnus and adjacent reticular formation medial to the lateral boundary of the pyramids at the level of the facial nucleus.

Data processing

After spike sorting and manual classification of cells, a window around each trial (−10 s to 100 s) was extracted relative to the alignment marker (heat onset or withdrawal behavior). Spikes within each window were aggregated across animals and trials to form the PSTHs for ON- and OFF-cells, providing a population-level view of stimulus-evoked activity. Trials were excluded from PSTH construction if the recording ended within the window or if the next event (heat onset or withdrawal) truncated the usable interval, resulting in a variable number of trials per PSTH.

For cells with available ongoing-activity recordings, we extracted a 1500 s segment (or 960 s for the shorter NEUTRAL-cell datasets) relative to the start of the first post–ongoing-activity trial. This ensured a minimum 230 s buffer after the initial trial, allowing transient evoked activity to dissipate before analyses of ongoing dynamics. Because noxious stimulation had no effect on NEUTRAL-cell firing rates, the full recording period from all neutral cells was used for GP modeling.

All analyses were conducted in Python v3.9. Instantaneous firing rates and basic statistics were computed using the Elephant package v1.1.1⁶². Gaussian processes were implemented using GPyTorch v1.13 and GPy v1.13.2^63,64. Code interfaced with the Neo framework (v0.13) and the SonPy package (v1.9.5) to convert Spike2 (.smrx) files into Neo objects ^65,66. Data was then standardized to the Neurodata Without Borders (NWB) standard ⁶⁷. PyMC v5.12 was used for Bayesian regression ⁶⁸.

We used Bayesian regression to model stimulation-evoked responses because the models contain relatively few parameters and the Bayesian framework provides full posterior uncertainty estimates—particularly useful for the switching timepoint. For ongoing activity, we employed maximum-likelihood type-II GP regression with a generalized Poisson likelihood. This approach was required because the long recording durations (> 1500 s), the nonstandard likelihood, and the multiple local minima induced by the periodic kernel made full Bayesian sampling computationally prohibitive.

Piecewise Bayesian modeling of ON- and OFF-cells

We modeled the activity of a pseudo-population composed of cells and trials from multiple animals, rather than modeling each cell or trial individually in a multi-level model. The reasons for this were twofold: firstly, with only 3-4 trials per neuron, and roughly 2 neurons per animal, there was not enough data to estimate a full multi-level model. Secondly, as mentioned above, our primary aim was to capture population-level structure. Single-electrode recordings are essentially random samples from the population of ON- and OFF-cells, and we therefore assumed the emergence of population level effects to be due to exchangeability and uncorrelated noise averaging out across neurons.

We used piecewise Bayesian regression to fit the dynamics of ON- and OFF-cells during trials with noxious stimulation. We binned spikes using a 50 ms time bin, assuming that the total spike count for ON- or OFF-cells per bin was Poisson distributed with a variable rate λ_t, dependent on parameters θ, such that:

For ON-cells, the response phase was modeled as an exponential rise ³³. We compared single- and double-timescale recovery functions. The double exponential model consisted of:

with continuity enforced by a corrected constraint:

Priors were:

The single exponential model for ON-cells was identical except for removal of the k_fast term.

For OFF-cells we compared a sigmoidal and a gaussian CDF function as the response function ³³ and compared a double timescale recovery with an exponential + linear recovery function. We constrained the sigmoid to pass through the point (t_s, u₀), joining the left and right side of the model. The general sigmoid has four parameters, slope k, midpoint t^*, lower asymptote a and upper asymptote b, such that:

The Gaussian CDF had form of:

where erf refers to the Gauss error function. For both the Gaussian CDF and the sigmoid, we chose t^* such that f (t_s) = u₀.

We fitted a, b and k directly; however, using a reparameterization of the sigmoid, we chose t^* so that f (t_s) = u₀. Starting from the relation:

The full model consisted of:

For the double exponential case, the right hand side of the model (t > t_s) was replaced with:

Priors were set as:

and additional priors for the double exponential OFF-cell recovery d and k₃ were defined identically to c and k₂ for the linear + exponential case, respectively. Markov chain Monte Carlo sampling using PyMC was used to sample from the posterior, using 1000 tuning steps, 4 chains, and 1000 samples per chain.

This Bayesian framework provided several advantages over prior methods ³³: (i) it does not require predefining “burst” and “pause” times; (ii) it dissociates peak and nadir timing from withdrawal; (iii) it yields full posterior distributions for extrema and firing rates at any time; and (iv) it enables principled comparison of alternative timescale parameterizations for recovery dynamics.

Gaussian process modeling

Firstly, as a preliminary period scan across potential timescales of ongoing firing, a GP with a periodic covariance function (see below) and Gaussian likelihood was applied to the smoothed firing rates of each neuron. We used 200 candidate periods (0.1-500 s) and a fixed lengthscale and variance of 1.0 and 5.0, respectively. A smoothing kernel and sampling period, both of width 2 seconds, were used to transform the spikes to a rate. The first and last 2 seconds of data were removed to avoid artifacts due to the smoothing process.

Secondly, slow structure in ongoing firing was quantified using latent GP models applied to binned spike counts (5 s bins). For each neuron, we fitted a sparse variational GP (200 inducing points) with a constant mean and a periodic covariance function:

where p denotes the oscillation period and ℓ_per the periodic length scale. A narrow Gaussian prior ℓ_per~ 𝒩(0.1, 0.1) was used to avoid over-smoothing in low-rate neurons. Kernel periods were initialized at 300 s, based on the period-scan analysis.

To map latent GP values to positive firing rates, a scaled Soft-Plus transform,

with A = 10 spikes/s was used instead of the exponential link to prevent numerical instabilities and speed up convergence across both low and high-rate regimes.

Spike count variability exhibited substantial over- and under-dispersion across cell classes. We evaluated a range of like-lihoods (Poisson, Conway–Maxwell–Poisson), kernels (periodic, spectral mixture), and link functions. A generalized Poisson model ^69,70 consistently provided stable fits and accommodated the broad dispersion range of ON-, OFF-, and NEUTRAL-cells. The generalized Poisson likelihood was parameterized by a parameter θ_t and dispersion parameter α:

θ_t was related to the mean µ_t by:

We placed a prior α~ 𝒩 (0, 0.5) and computed θ_t = µ_t (1 −α) at each time point. A small contamination term (p_noise = 0.001) was added to account for rare outlier bins:

Coherence estimation

We first binned each spike train x and heart rate y using a 2 s bin size. We split the recording into 500 s segments with a 50% overlap, which gave each segment frequency resolution 1/T = 0.002 Hz. We demeaned each segment then used a multitaper spectral estimator with discrete prolate spheroidal sequences (DPSS) (standardized half bandwidth NW = 3.5, number of tapers K = 6) to calculate the coherence ⁷¹, using the formula:

Segments were then individually phase-randomized to build a surrogate coherence distribution, and used a max-statistic approach to control the Family-Wise Error Rate (FWER) across all frequencies. The significance threshold was set at α = 0.001.

This combination of parameters gave a bandwidth = B = NW/T = 0.007 Hz, meaning that the frequency f = 0.0033 Hz (T = 5 minutes) found by the Gaussian Process was smeared across the low-frequency end of our coherence estimator. It was therefore not possible to determine whether the significant coherence found was due to oscillations at the same frequency as the GP period, due to the constraints of recording length and robust spectral estimation. Between individual cells, we observed some visual and spectral differences across all cell types between 0.01 Hz “fast” coherence and 0.0033 Hz “slow” coherence but these could not be quantified.

Data availability

All data related to this publication can be found on the DANDI data repository ⁷².

Acknowledgements

FM is funded by a MRC Career Development Award (MR/T010614/1), a UKRI Advanced Pain Discovery Platform grant (MR/ W027593/1), and a EPSRC/MRC Programme Grant (UKRI1970). CA is funded by an EPSRC DTP (EP/W524633/1). The Heinricher laboratory is supported by grants from the National Institutes of Health (NS098660, NS120486) to MMH. CCdP is supported by a fellowship from the National Institute of Dental and Craniofacial Research (F31 DE030677).

Additional information

Consent for publication

All authors gave their consent for publication.

Code availability

Code is available on the Zenodo code repository and on GitHub ⁷³.

Author Contributions

This study was performed via a collaboration between the Heinricher Lab at Oregon Health and Science University nd the University of Cambridge. Carl Ashworth, Mary Heinricher and Flavia Mancini conceived the study. Melissa Martenson, Caitlynn De Preter, and Zhigang Shi performed he experiments. Carl Ashworth analyzed the data. Carl Ashworth, Mary Heinricher and Flavia Mancini wrote the manuscript. All authors reviewed the manuscript.

Funding

UK Research and Innovation (UKRI) (EP/W524633/1)

• Carl Ashworth

HHS | National Institutes of Health (NIH) (NS120486)

• Mary Heinricher

HHS | NIH | National Institute of Dental and Craniofacial Research (NIDR) (F31 DE030677)

• Caitlynn De Preter

HHS | National Institutes of Health (NIH) (NS098660)

• Mary Heinricher

UK Research and Innovation (UKRI) (MR/T010614/1)

• Flavia Mancini

UK Research and Innovation (UKRI) (MR/W027593/1)

• Flavia Mancini

UK Research and Innovation (UKRI) (UKRI1970)

• Flavia Mancini

Additional files

Supplementary Information