Abstract
Sleep is known to drive the consolidation of motor memories. During nonrapid eye movement (NREM) sleep, the close temporal proximity between slow oscillations (SOs) and spindles (“nesting” of SO-spindles) is known to be essential for consolidation, likely because it is closely associated with the reactivation of awake task activity. Interestingly, recent work has found that spindles can occur in temporal clusters or “trains.” However, it remains unclear how spindle trains are related to the nesting phenomenon. Here, we hypothesized that spindle trains are more likely when SOs co-occur in the prefrontal and motor cortex. We conducted simultaneous neural recordings in the medial prefrontal cortex (mPFC) and primary motor cortex (M1) of male rats training on the reach-to-grasp motor task. We found that intracortically recorded M1 spindles are organized into distinct temporal clusters. Notably, the occurrence of temporally precise SOs between mPFC and M1 was a strong predictor of spindle trains. Moreover, reactivation of awake task patterns is much more persistent during spindle trains in comparison with that during isolated spindles. Together, our work suggests that the precise coupling of SOs across mPFC and M1 may be a potential driver of spindle trains and persistent reactivation of motor memory during NREM sleep.
Significance Statement
We found evidence for two types of sleep spindle organization in primary motor cortex (M1)—trains and isolated. Train spindles exhibited stronger M1 slow-oscillation (SO) locking and M1 SO-spindle nesting. Precise coupling between prefrontal cortex and M1 SOs predicted M1 train spindles and was correlated with increased nesting. Spindle trains were also associated with enhanced M1 SO-locked spiking modulation and sustained task-related ensemble memory reactivation. Our results indicate that cross-area SO dynamics may be key for M1 spindle trains and sustained memory reactivations. This advances our understanding of how global activity patterns during sleep may promote consolidation of memories.
Introduction
Sleep is known to be important for memory consolidation (Buzsáki, 1989; Stickgold, 2005; Marshall et al., 2006; Marshall and Born, 2007; Tononi and Cirelli, 2014; Poe, 2017). There is also growing evidence for the importance of sleep in the consolidation of new motor memories (Walker et al., 2002; Rasch and Born, 2013; Ramanathan et al., 2015; Gulati et al., 2017; Boutin et al., 2018, 2024; J. Kim et al., 2019, 2023; Klinzing et al., 2019). Non-rapid eye movement (NREM) sleep, in particular, has been shown to be an important window for the consolidation of motor memories by driving the reactivation of awake task-related neural ensembles (Walker et al., 2002; Ramanathan et al., 2015; J. Kim et al., 2019, 2023). Such reactivation is particularly prevalent during spindles. Spindles have been found to occur across regions known to be important for learning, in general, and motor learning, in particular (Ramanathan et al., 2015; Fogel et al., 2017; Vahdat et al., 2017; Boutin et al., 2018; J. Kim et al., 2019; Silversmith et al., 2020). Importantly, there are studies showing that reactivations during “nested” spindles, which are temporally locked to slow oscillation (SO) UP states, are essential for consolidation (J. Kim et al., 2019). Despite work demonstrating the importance of SO-spindle nesting, recent findings suggest that spindles exhibit inherent rhythmicity in their spectral power that is associated with motor memory consolidation during sleep (Antony et al., 2018; Boutin et al., 2024). The clustering of spindles over time, also referred to as “spindle trains,” suggests that there may be longer timescales of functional organization, i.e., beyond the notion of nesting between a single SO and a single spindle. It remains unclear how SOs, which tend to originate in the prefrontal cortex, interact with spindle trains. Moreover, it remains unclear how they affect reactivation events.
While much attention has been devoted to studying the association between “local” SOs (i.e., in primary cortical areas such as the motor cortex), spindles, and motor learning, less is known about the role of SOs originating from the medial prefrontal cortex (mPFC) in coordinating local SOs and spindles (J. Kim et al., 2019; Silversmith et al., 2020). It is not clear whether SOs, which are known to start in mPFC and propagate toward posterior areas of the cortical mantle (Massimini, 2004), might have a role in the organization of spindle trains. For example, propagation of mPFC SOs toward posterior cortical areas may prime thalamocortical circuits for the generation of cortical spindles (Latchoumane et al., 2017; Hay et al., 2021). In this case, we would expect that mPFC SOs would exhibit an increased propensity toward close temporal coupling with spindles and with train spindles, in particular. Based on the known importance of SO-spindle nesting in driving memory consolidation, this might also suggest a role for spindle trains in memory consolidation.
Here we performed simultaneous recordings in mPFC and primary motor cortex (M1) in rats during NREM sleep after reach-to-grasp (RTG) task training. We first identified spindle trains and distinguished them from isolated spindles. Strikingly, the co-occurrence of mPFC SOs with M1 SOs was closely associated with the occurrence of spindle trains in M1. Notably, we also found evidence for coincident SOs in mPFC and M1 during a spindle train, suggesting that mPFC and M1 interactions are continuous during spindle trains. Importantly, there were significantly longer and sustained reactivation events during a spindle train, suggesting that spindle trains may be an important component of sleep-dependent processing. Overall, our results indicate that following a subset of mPFC SOs, there appears to be globally coordinated SOs in mPFC and M1 that are colocked with spindle trains in M1.
Materials and Methods
Animals
All animal training and surgical procedures were conducted in accordance with protocols approved by the Institutional Animal Care and Use Committee at the San Francisco Veterans Affairs Medical Center. Adult male Long–Evans rats sourced from Charles River Laboratories and aged between 3 and 6 months old (n = 6, 300–400 g) were housed in a controlled temperature room with a 12 h light/12 h dark cycle. All experiments were conducted during the light cycle (lights on at 6 A.M., lights off at 6 P.M.). The data presented in this work has been previously utilized in other published work (J. Kim et al., 2023).
Electrode implantation surgeries
Prior to electrode implantation, each animal's hand preference was determined by observing their arm preference during reaching behavior in response to the experimenter holding a pellet of food with tweezers while animals were allowed to roam freely in the behavior recording chamber. Electrode implantation surgeries were then performed under isoflurane anesthesia (5% for induction, 1–3% for maintenance), and body temperature was maintained at 37°C with a Kent Scientific Right Temp heating pad. Preoperative medications included atropine sulfate (0.02 mg kg−1 of body weight, i.p.), dexamethasone (0.5 mg kg−1, i.p.), and lidocaine (3 cc, infiltration around scalp incision). A single 32-channel Tucker-Davis Technologies (TDT) microwire array was inserted into Layer 5 of the M1 in the forelimb area centered at 3.0 mm medial–lateral, 0.5 mm anterior–posterior, and 1.4–1.8 mm dorsal–ventral (all measurements relative to the bregma). A second 32-channel array was placed into mPFC centered at 1.0 mm medial–lateral, 3.5–4.0 mm anterior–posterior, and 3.5–4.0 mm dorsal–ventral (all measurements relative to the bregma). No particular area of mPFC was targeted as we were interested in SO propagation dynamics across the brain rather than the content encoded by mPFC during behavior. Reference and ground wires were twisted around a skull screw inserted over the cerebellum. Final electrode depth was determined using intraoperative recordings to maximize spiking and local field potential (LFP) signal quality at the time of array implantation. The craniotomy was closed using KwikCast and a headcap was built using acrylic dental cement. Postoperative analgesia was provided through buprenorphine (0.02 mg kg−1, i.p.), meloxicam (0.2 mg kg−1, in water for 5 d), and dexamethasone (0.5 mg kg−1, in water for 5 d). Antibiotic coverage was provided through trimethoprim sulfadiazine at 15 mg kg−1 in water for 5 d. All animals were allowed to recover for at least 7 d before the start of experiments.
RTG task training and sleep recordings
After animals recovered from microwire implantation surgeries, they were handled 3–5 d before commencement of food restriction (started 2 d before RTG training start). Animals were fed 50 g kg−1 of chow at the end of the day (∼5 P.M.), and we ensured that each animal's body weight did not drop below 90% of their initial weight. Animals were also acclimated to a custom-made, robotic behavioral box during this time period. Animals were trained to reach for the pellet during the RTG by allowing the automated pellet presentation system (Wong et al., 2015) to cycle through several rounds of RTG trials (details of trial structure have been previously described elsewhere; Wong et al., 2015; Lemke et al., 2019, 2021; Kondapavulur et al., 2022; J. Kim et al., 2023). No human intervention was used during the pretraining period. Animals typically showed a successful trial 1–3 d after the start of automated training, and only behavior sessions starting on that day and going forward were included in the learning curve and neural data. It took 12–13 d (13 d in five animals and 12 d in one animal) to reach plateau performance (see J. Kim et al., 2023 for learning curve). Each recording sessions consisted of a 2 h pretraining sleep session, 100 RTG trials, and then another 2 h post-training sleep session. LFP and spiking activity were recorded continuously during all phases of the task. We used LFP from the post-training sleep periods, and all segments of detected NREM sleep were used for analyses. Average duration of NREM sleep was 56.6 ± 13.1 min per sleep session [mean ± standard deviation (SD) from N = 18 sleep sessions]. There was no differential weighting of data by time within each sleep period. The last three learning sessions were used for analyses, as they exhibited stable reach success rates, but sessions from individual rats were otherwise independent. The reach success rate was scored as the percentage of trials that the rat was able to successfully retrieve the pellet with its arm and place the reward in its mouth; this metric was calculated by manual scoring of the RTG behavior video recordings.
Electrophysiology and spike sorting
We recorded synchronous mPFC and M1 extracellular neural activity, including single units and LFPs, using a TDT RZ2 system. Spiking data were sampled at 24,414 Hz, and LFP data were sampled at 1,018 Hz (followed by online bandpass filtering from 0.1 to 500 Hz with a notch filter at 60 Hz). For spike sorting, snippets of online data that crossed a high signal-to-noise threshold (at least 3.5 SD from the signal's 300 Hz high-pass filtered mean) were deemed spiking events. These snippets were then spike sorted using MountainSort with each day's recording having its own single-unit sorting (i.e., we made no attempts to track singe units across days). Following automated sorting, minimal manual merges and rejections of putatively identified single-unit clusters was performed. Only clearly identifiable units along with sessions of a single day, with good waveforms and a high signal-to-noise ratio, were used.
LFP postprocessing and sleep stage segmentation
Broken or noisy LFP channels were noted online for every recording day and then excluded during off-line analysis. After rejection of bad channels, 23–32 M1 LFP channels per session were used for LFP analyses (29.7 ± 2.2, mean number of channels ± SD, N = 18 recording sessions). In PFC, 18–32 LFP channels per session were used for LFP analyses (28.3 ± 3.2, mean number of channels ± SD). Each LFP channel's data were Z-scored, and then a mean was taken across all channels (but separately for the mPFC and M1 electrode arrays). This mean trace was then used for all subsequent analyses. In order to segment wake versus NREM sleep, we followed the procedure outlined by J. Kim et al. (2019). In brief, each LFP channel was segmented into nonoverlapping 6 s windows. In each window, the power spectral density was computed and averaged over the delta/SO (0.1–4 Hz) and gamma (30–60 Hz) frequency bands. Then a k-means classifier was used to cluster epochs into two clusters, NREM sleep and REM/awake. Only long (0.30 s, five consecutive windows) epochs of sleep were analyzed.
SO and spindle detection
After segmentation of sleep LFP into NREM versus wake epochs, we next proceeded to segment SOs and spindles separately in mPFC and M1 LFP. The SO detection algorithm is similar to those used previously (J. Kim et al., 2019). To detect the ∼1 Hz SO, the mean LFP signal was filtered in a low-frequency band (second-order, zero phase-shifted, high-pass Butterworth filter with a cutoff at 0.1 Hz followed by a fifth-order, zero phase-shifted, low-pass Butterworth filter with a cutoff at 4 Hz). Next, all positive-to-negative zero crossings during NREM sleep were identified, along with the previous peaks, the following troughs, and the surrounding negative-to-positive zero crossings. Each identified epoch was considered an SO if the peak was in the top 85% of peaks, the trough was in the top 40% of troughs, and the time between the negative-to-positive zero crossings was 0.300 ms but did not exceed 1 s.
The spindle detection applied here is similar to the algorithm used in J. Kim et al. (2023). The mean LFP signal across mPFC and M1 was separately filtered in the spindle band (10–16 Hz) using a zero phase-shifted, third-order Butterworth filter. A smoothed envelope was calculated by computing the magnitude of the Hilbert transform of this signal and then convolving it with a Gaussian window (α = 2.5). Next, we determined two thresholds for spindle detection based on the mean and SD of the spindle band envelope during NREM sleep (lower, 1.5 SD; upper, 2.5 SD). Epochs in which the spindle envelope exceeded the upper threshold for at least one sample and the spindle power exceeded the lower threshold for at least 500 ms were considered spindles. In other words, the smoothed instantaneous envelope of the filtered spindle LFP must exceed the 1.5 * SDenvelope threshold for at least 500 ms, and at least one or more points of the envelope which are already above the 1.5 * SDenvelope threshold must also exceed the 2.5 * SDenvelope threshold. Finally, spindles that were sufficiently close in time (<300 ms) were combined. For each spindle epoch, the peak of the spindle band LFP was identified. All identified spindles were aligned to this peak for generating average spindle waveforms and measuring SO-spindle nesting. As a final step for both SO and spindle detection, only those events which occurred during NREM sleep were kept for further analyses.
Measurement of time between mPFC–mPFC, M1–M1, and mPFC–M1 SOs
To determine whether the timing relationship between mPFC-to-M1 SOs was different as compared with the timing relationship between mPFC-to-mPFC and M1-to-M1 SOs, we first calculated the time between each mPFC SO to its immediately following SO (and the same was repeated for M1 SOs). We also calculated the time between each mPFC SO and the very next M1 SO that followed it; here the search for the next M1 SO following each mPFC SO was done in a causal way given the literature suggesting that, in humans and rodents, SOs propagate from anterior to posterior across the cortical mantle (Massimini, 2004; Niethard et al., 2018), and thus we only selected for M1 SOs that followed mPFC SOs. In order to simplify visualization and prevent the counting of SOs across two contiguous NREM epochs, we set a maximum cutoff of 10 s for the interevent histograms and for further analyses that relied on SO–SO interevent times.
Mixture of Gaussian clustering of mPFC SO to M1 SO timing
In order to cluster the mPFC–M1 SO time distribution into three discrete categories—two capturing the bimodal “ultra short” peak at <100 ms plus another “short” peak at ∼1–2 s and a third cluster which we use as a data-driven control condition—we fit a mixture of Gaussian (MOG) model to the natural log-transformed distribution of mPFC–M1 SO times. We fit the model using scikit-learn (Pedregosa et al., 2011)'s GaussianMixture class with parameters n_components, 3, and covariance_type, “full.” We then took the mean and variance across all animals for each of the three Gaussian clusters and exponentiated the averages to get the mean and variance for the original, nonexponentiated mPFC–M1 SO distribution. Then we segmented the mPFC–M1 SO epochs by grouping individual data points of the mPFC–M1 SO distribution that were within the mean ± 1 SD of each cluster center. To calculate the percentage of M1 SOs locked to mPFC SOs within the time bounds of the clusters, we calculated the fraction of all M1 SOs (within a given sleep session) which were preceded by an mPFC SO within the set timeframe of each cluster.
Spindle train segmentation
In contrast to human studies of EEG spindle trains (Boutin et al., 2022, 2024), there are no established criteria for segmentation of spindle trains in intracortical rodent LFP signals. We thus developed criterion for determining train spindle status in our dataset by using expectation–maximization to fit an MOG model (Pedregosa et al., 2011; with n = 2 clusters) to the natural log-transformed distribution of time between spindle events. We then exponentiated the resulting mean and SD of the fitted Gaussians and plotted them on the regular, base-10 axis of spindle interevent times (Fig. 1D). By calculating the intersection of the two Gaussian curves, we arrived at the equivalent of a Bayesian decision boundary segmenting train versus isolated spindles into two groups (mean ± SD of all intersections, 2.78 ± 0.05 s). Our percentage of train (∼55%) versus isolated (∼45%) spindles matches prior reports based on human EEG data (Boutin et al., 2022; Solano et al., 2022). Hence, we define spindle trains as ≥2 spindles ≤2.78 s apart.
Locking of M1 spindles relative to both M1 and mPFC SOs
For analyses shown in Figures 3D and 4B which involve the three-way relationship between mPFC SOs, M1 SOs, and M1 spindles, we start by finding the nearest M1 SO to each M1 spindle. Here, consistent with prior literature, we allow M1 spindles to occur either before or after the UP state of M1 SOs (J. Kim et al., 2019); however, we still enforce the formerly mentioned constraint of ensuring that mPFC SOs precede M1 SOs. Moreover, if a single M1 spindle–M1 SO complex could lock to multiple mPFC SOs, we only locked onto the mPFC SO that was closest to the M1 SO that each respective M1 spindle was most proximal to.
M1 SO-spindle nesting quantification
For those analyses reporting the percentage of M1 spindles closely locked to M1 SOs, we calculated the percentage of these closely locked, otherwise known as “nested,” spindles by considering a spindle as nested if its peak amplitude occurred within – 0.5 to +1.0 s of a M1 SO UP state (J. Kim et al., 2019). To report the percentage of M1 spindles that were nested, we then simply calculated the number of M1 spindles with an M1 SO UP state that was within the spindle nesting interval above. In order to generate the spindle nesting time histograms in Figure 4C, we measured the time from each M1 spindle peak to that peak's nearest M1 SO UP state. To calculate the probability of a M1 spindle following mPFC SO (irrespective of whether there was an intervening M1 SO, Fig. 3D), we binned time into 50 ms bins and calculated the probability of a M1 spindle within each bin separately for train and isolated spindles. The resultant distribution was then smoothed with 10-point Gaussian window (α = 2.5) for improved visualization. To report the continuous probability histogram of M1 spindle-to-M1 SO nesting in Figure 3F, we binned time into 75 ms chunks and then calculated the percentage of train versus isolated spindles nested to M1 SOs. The resultant distribution was then lightly smoothed with a 30-point Gaussian window (α = 2.5) for improved visualization.
Quantifying mPFC–M1 SO co-occurrence within spindle trains
To explore the role of mPFC–M1 SO synchronization within a spindle train, we selected all PFC and M1 SOs that occurred between the first and last spindle of each spindle train. We then determined, in a causal manner, the time between each PFC SO and its next associated M1 SO. This distribution of timing values was subsequently visualized as a proportion histogram, binned in 62.5 ms bins (Fig. 6D).
SO-locked PETH generation
To calculate spiking modulation during the DOWN/UP states of M1 SOs, we first calculated the PETH of all M1 units relative to each M1 SO detected during a given sleep session on a −1.50 to +0.75 s window around the UP state of each M1 SO. We then binned the spiking activity of each unit and for each SO using 1 ms bins, smoothed the resultant spike trains with a 50 ms Gaussian window (α = 2.5), and binned the data into 100 ms bins for further analysis. This was followed by taking the Z-score of each smoothed PETH and then separating the PETHs based on the mPFC–M1 SO timescale categories. We then calculated an average, Z-scored M1 SO-locked PETH. In order to quantify minimum and maximum spiking modulation of M1 SO DOWN and UP states, respectively, we calculated the minimum DOWN state-locked mean Z-scored firing rate from −0.45 to −0.25 s relative to the M1 SO UP state as well as the maximum UP state-locked mean Z-scored firing rate from −0.20 to 0.00 s relative to the M1 SO UP state. The mean ± SEM of these values was then plotted in Figure 7C.
GPFA trajectory calculation
In order to track reactivation of task-related neural during sleep, we first used Gaussian-process factor analysis (GPFA; Yu et al., 2009) to reduce the dimensionality of the neural data during the awake reaching task. GPFA was performed with a time bin of 15 ms, and six factors determined leave-one-out cross-validation approach. For each recording session, we concatenated binned M1 spike trains into a neuron by time (N × T) matrix for all RTG trials and then Z-scored the resultant spike train. Spike counts used for GPFA trajectory calculation were selected specifically around the RTG period and adjusted for the individual reach onset-to-pellet-touch timing of each animal in the dataset (−200 to 400 ms from reach onset for four animals; −200 to 600 ms from reach onset for one animal; and −200 to 1,000 ms from reach onset for the last animal). After fitting GPFA trajectories to these matrices, we kept the top three factors for further analyses as they accounted for >85% of shared variance explained in the M1 spiking activity. GPFA analyses were only performed within-day as we do not hold the same single units across days.
Reactivation analysis
The neural spike trains related to reaching were dimensionally reduced using GPFA. The resulting GPFA factors were then used to measure task-related activity reactivation during train and isolated spindles in post-training NREM sleep. A previously published method was employed to measure the magnitude of reactivation values (J. Kim et al., 2023). During a spindle peak in post-training sleep, the binned spike trains (bin length, 15 ms) were projected onto the template space, which was defined by the top three GPFA factors derived from neural activity during the reaching task. The projection was a linear combination of Z-scored binned spike trains from post-training sleep with the template space calculated from the neural activity during reach training. We searched for motor memory reactivations during post-training NREM sleep by projecting binned spike trains onto a “template space” defined by GPFA factors from the reaching task. To account for temporal compression or extension of reactivations, we explored different window sizes (75–495 ms) and time lags (−95 to 420 ms) from spindle peaks. For each spindle, we selected the combination of size and lag that maximally correlated with the session-specific reach template, resulting in a neural trajectory reactivation with a “reactivation correlation (R)” value. Here reactivation R is the correlation coefficient between a neural trajectory reactivation during a spindle and the reach template. We then compared reactivation R values between train spindles (2–5 consecutive spindles) and isolated spindles, where we simulated 1–4 subsequent events following an isolated spindle to match the interevent intervals and counts of train spindles in each session. This allowed us to investigate the impact of task-related activity during post-training NREM sleep on motor memory reactivations during train versus isolated spindles.
Statistics
All plots and statistical analyses are performed across N = 18 sleep sessions pooled across six animals (three sessions per animal). Unless otherwise noted, all differences in means are compared using one-/two-way ANOVAs with Tukey multiple-comparison corrected post hoc tests. Differences in distributions are measured using Kolmogorov–Smirnoff tests. All analyses were performed using Python v3.10 except for the GPFA trajectory reactivation analyses which were calculated in MATLAB R2019a.
Results
Identification and properties of train spindles
Simultaneous recordings were conducted in mPFC and M1 (Fig. 1A) from six separate rats trained to perform a skilled RTG task. Each recording session started with 2 h of pretraining sleep, followed by 100 RTG trials, and concluded with 2 h of post-trainings sleep. Three post-training sleep recording sessions were used for all analyses from each of six individual rats (N = 18 sleep sessions total). After segmenting periods of NREM sleep from wakefulness plus detecting SOs (in mPFC and M1) and sleep spindles (in M1), we observed periods of sleep LFP that were enriched in temporally adjacent sleep spindles (Fig. 1B, spindle trains) and periods of LFP with a relative paucity of sleep spindle incidence (Fig. 1C, isolated spindles). Indeed, human studies have identified spindle trains using extracortical EEG (Boutin et al., 2022, 2024). Here, we established criteria for the segmentation of spindle trains in intracortical rodent LFP signals. We identified the membership of individual spindles to spindle trains by clustering the distribution of interspindle event times using an MOG (Fig. 1D density histogram of interspindle times). The mean Bayesian decision boundary (Fig. 1D, boxplot; N = 18 sleep sessions) between train and isolated spindle distributions across animals yields a time of 2.78 s between spindle events, and we thus define spindle trains as ≥2 spindles ≤2.78 s apart. When applied across all animals and sleep sessions, we found that 45% of all spindles are classified as “isolated” spindles, while the remainder are classified as “train” spindles (Fig. 1E). Interestingly, our approach to segmenting train spindles, when applied to intracortical LFP signals, matches the percentage of train versus isolated spindles that are typically seen in human EEG data (Boutin et al., 2022; Solano et al., 2022). Figure 1F demonstrates the time between adjacent sleep spindles within trains, while Figure 1G shows the time between the end of one spindle train and the start of the next spindle train. Spindle trains were relatively evenly distributed across NREM sleep in each post-training sleep epoch with only a slightly increased prevalence during the first 10–20% of sleep and the last 30% of sleep (Fig. 1H).
Identification and properties of spindle trains. A–C, Simultaneous sleep recordings in mPFC (blue square) and M1 (orange square) with three representative traces of train and isolated sleep spindles shown in B and C, respectively. Traces shown as Z-scored LFP and vertical scale bars represent 3 Z-values in B and C. D, Histogram of time between all spindles in M1 (maroon line, 18 sleep sessions). MOG clustering with n clusters = 2 was used to segment train from isolated spindles. Horizontal boxplot captures the point of intersection (blue scatters, 18 sleep sessions). Gaussian distribution from each animal shown in gray (for train spindles) and black (for isolated spindles). E, Histogram showing proportion of isolated spindles and train spindles across 18 sleep sessions. F, Density histogram of interspindle event time for M1 train spindles (data shown as mean ± SEM with data from individual sleep sessions in gray lines). G, Density histogram of time between the end of one spindle train and the start of another (data shown as mean ± SEM with data from individual sleep sessions in gray lines). H, Incidence of spindle trains as a function of time in NREM sleep for each sleep session (data shown as mean ± SEM with data from individual sleep sessions in gray lines). I, Left panel, Comparison of spindle amplitude for isolated (pink) and train spindles (maroon). Right panel, Comparison of spindle duration for isolated (pink) and train spindles (maroon); *p = 0.012, two-sample t-test. J, Bar graph, Proportion of train versus isolated spindles closely nested to an M1 SO UP state, i.e., within −0.5 to +1.0 s spindle peak to M1 SO UP state (data shown as mean ± SEM). Inset histogram, Histogram of time from each M1 SO UP state to the nearest M1 spindle (maroon, train spindles; pink, isolated spindles; ***p < 0.001, Kolmogorov–Smirnoff test, 18 sleep sessions).
While isolated and train spindles did not have statistically significant differences in amplitude (Fig. 1I, left violin plots with inset boxplots; isolated spindles, 0.685 ± 0.019 Z-scored amplitude; train spindles, 0.696 ± 0.022 Z-scored amplitude; data as mean ± SEM; two-sample t test, ts = −0.390; p = 0.699; df = 33.314; N = 18 sleep sessions), isolated spindles were slightly longer duration as compared with train spindles (Fig. 1I, right violin plots with inset boxplots; isolated spindles, 0.840 ± 0.021 s; train spindles, 0.775 ± 0.013 s; data as mean ± SEM; two-sample t test, ts = 2.690; p = 0.012; df = 27.764; N = 18 sleep sessions).
We next examined the timing relationship between train and isolated spindles and their closest corresponding local M1 SO. We found that train spindles were more likely to nest just before the local M1 SO UP state as compared with isolated spindles which were more likely to nest immediately after the M1 SO UP state (Fig. 1J inset; Kolmogorov–Smirnoff test; test statistic, 0.789; p < 0.001). Next, we quantified the percentage of M1 train versus isolated spindles that nest within [−0.5, +1.0] s of an M1 SO UP state (i.e., the canonical definition of SO-spindle nesting; J. Kim et al., 2019). Here we found that train spindles were statistically significantly more likely to nest with the UP state of M1 SOs (Fig. 1J, train spindles, 39.0 ± 0.76%; isolated spindles, 33.5 ± 1.02%; t test, ts = 2.944; p = 0.006).
mPFC SO proximity predicts nesting of train spindles to M1 SOs
Having observed increased nesting of train spindles to M1 SOs, we next sought to determine if (1) mPFC SOs promote the generation of M1 train spindles and (2) whether mPFC–M1 SO–SO coupling promotes nesting between M1 SO and M1 spindles. Figure 2A shows a representative snippet of mPFC and M1 LFP activity around an isolated M1 spindle. Plotted is the low-frequency (0.1–4 Hz) and spindle band (10–15 Hz) LFP in both mPFC and M1. Identified SOs are shown in purple for mPFC and maroon for M1 (Fig. 2A,B, top and middle trace). Identified M1 isolated spindles are labeled in pink on Figure 2A (bottom trace). Figure 2B follows the same plotting convention as Figure 2A though for a segment of LFP activity centered around a M1 spindle train (the same animal and sleep session plotted as in Fig. 1B, spindle trains in maroon on bottom trace). Here we noticed the presence of mPFC SOs in close temporal proximity to train spindles as compared with the absence of mPFC SOs for the example of an isolated M1 spindle.
Interactions between spindle trains and mPFC SO–M1 SO coupling. A, Representative snippet of mPFC and M1 LFP activity around an isolated M1 spindle. Low-frequency (0.1–4 Hz) and spindle band (10–15 Hz) LFP in both mPFC (blue top line) and M1 (orange middle and bottom lines). Identified SOs are shown in purple for mPFC and maroon for M1 (top and middle trace). M1 isolated spindles labeled in pink (on bottom trace). B, The same plotting convention as A though for a segment of LFP activity centered around a M1 spindle train (the same animal and sleep session plotted as in A, spindle trains in maroon on bottom trace). Traces shown as Z-scored LFP and vertical scale bars represent 3 Z-values. C, Schematic illustrating the measurement of temporal proximity of mPFC SO to M1 spindle interevent times for train and isolated spindles. D, Continuous probability distribution showing proportion of M1 train versus isolated spindles as a function of time from the nearest mPFC SO. Maroon line represents train M1 spindles, pink line for isolated M1 spindles; solid lines with shading represent mean ± SEM; faint lines in the background represent sleep sessions from individual rat sessions (18 sleep sessions, maroon for train spindles and black for isolated spindles). E, Schematic illustrating the measurement of M1 SO-spindle nesting as a function of proximity to the nearest preceding mPFC SO. F, A plot with the percentage of M1 SO-spindle nesting events relative to time from the nearest mPFC SO (same plotting convention as in D).
We next sought to quantify the relationship between the mPFC SOs relative to M1 spindles (either train or isolated, Fig. 2C). When measured continuously from the UP state of an mPFC SO, we noticed that M1 train spindles were much more likely to occur in approximately the first 1 s after the mPFC SO UP state in comparison with the M1 isolated spindles (Fig. 2D). We further examined the relationship between the temporal proximity of mPFC SOs to the nesting percentage of local M1 SOs and M1 spindles (either train or isolated, Fig. 2E). A similar analysis examining the fraction of M1 spindles nested to M1 SOs, as a function of both temporal proximity to mPFC SOs and spindle train membership, revealed a similar pattern—the closer an mPFC SO was to a M1 SO-spindle complex, the greater the likelihood of M1 train spindle nesting to the UP state of a local M1 SO (Fig. 2F).
We next sought to quantify the impact of mPFC–M1 SO temporal proximity on local M1 SO-spindle nesting. We began by first calculating the distribution of mPFC-to-M1 SO time differences (Fig. 3A; histogram of mPFC–M1 SO time from 0 to 10 s on the main panel; the inset shows histogram of values from 0 to 0.5 s). We calculated this distribution in a temporally “causal” direction—i.e., an M1 SO was not allowed to occur after a given mPFC SO in alignment with previously published work showing that SOs propagate from frontal cortical areas posteriorly across the cortical mantle (Massimini, 2004; Oyanedel et al., 2020). Here we observed a striking bimodality in mPFC–M1 SO timing with a peak from 0 to 0.1 s and then a second peak from ∼1 to 2 s. In order to segment these two classes of mPFC–M1 SO interactions, we performed an MOG clustering on the log-transformed distribution of mPFC–M1 SO timing (Fig. 4A–E). We fit the MOG model with n = 3 clusters; two clusters to capture the ultra short (<100 ms) and short (1–2 s) timing interactions plus a third long cluster fit by the model that we use as a data-driven control. Figure 3A demonstrates the resultant discretization of the mPFC–M1 SO timing distribution (average cluster center ± 1 SD: ultra short [0.006, 0.073], short [0.754, 2.585], and long [3.094, 7.44]).
mPFC SO proximity predicts nesting of train spindles to M1 SOs. A, Histogram of time between mPFC SOs and M1 SOs. An MOG clustering was used to separate mPFC–M1 SO timing interactions into three clusters (ultra short, short, and long latency clusters are colored dark green, salad green, and light green, respectively). B, Comparison of the percentage of M1 train and isolated spindles nested for each of the three clusters (data shown as mean ± SEM across 18 sleep sessions). C, Histogram of the time between M1 SOs and M1 spindles for the three classes.
Characteristics of mPFC SO–M1 SO timing. A, Histogram of natural log-transformed values of causal times between mPFC SOs and M1 SOs (dark gray bars). Mixture of Gaussian fitted distributions plotted in dark green for ultra short, salad green for short, and light green for long mPFC–M1 SO timing interactions. Red line plots the mean of all three Gaussian distributions summed together. B, Boxplots showing centers of each Gaussian cluster across animals and session following MOG clustering (fit with n clusters = 3, plotted on a base-10 scale) of the mPFC–M1 SO time distribution. C, Proportion of M1 SOs locked to mPFC SOs at each of the three timescales identified using MOG clustering. D, Continuous probability distribution showing chance of M1 SO UP state occurring at a given time point relative to mPFC SO UP states. E, Visual examples of mPFC SO–M1 SO coupling and its relationship to M1 SO-spindle nesting. Top set of three traces show ultra short mPFC SO–M1 SO coupling, middle set of three traces shows short mPFC SO–M1 SO coupling, and bottom set of three traces show long mPFC SO–M1 SO coupling. mPFC SOs in purple, M1 SOs in maroon, and nested spindles maroon. All traces shown as Z-scored LFP and vertical scale bar represents 3 Z-values.
We next analyzed whether mPFC SO–M1 SO latency impacted the nesting of train versus isolated spindles to M1 SOs. In Figure 3B, we noticed that train spindles were especially strongly nested to M1 SOs when the M1 SOs were preceded by an ultra short or short (but not long) latency mPFC SO (two-way ANOVA, interaction of latency and train status, F, 24.198; p = 3.759 × 10−11; N = 18 sleep sessions). Figure 4E provides a visual example, illustrated with traces of raw data, of the nesting phenomenon summarized in Figure 3B. Tukey post hocs on the interaction between train status and mPFC–M1 SO latency demonstrated statistically significant difference for ultra short (p < 0.001) and short (p < 0.001) but not long (p = 0.638) mPFC–M1 SO timescales. These results suggest that the proximity of an mPFC SO to an M1 SO-spindle complex tends to be a strong predictor of train spindle nesting to the UP state of M1 SOs which is further reinforced in Figure 3C where we notice stronger clustering of M1 SO to M1 spindle times right at the UP state of M1 SOs for ultra short and short mPFC SO–M1 SO pairs.
Co-occurrence of mPFC and M1 SOs within spindle trains
What might be a mechanistic consequence of spindle trains in cross-area dialogue during NREM sleep? We explored synchronization between mPFC and M1 SOs when bracketed by the timing of the spindle train itself. We began by isolating spindle trains that were two, three, or four spindles in length. Then we measured the interval from each mPFC SO to the very next M1 SO but only for those mPFC–M1 SO pairs that occurred within the peak-to-peak time bounds of each adjacent spindle pair (Fig. 5A shows a schematic demonstrating the calculation of mPFC–M1 SO timing within adjacent train spindles). Figure 5B presents a raster of mPFC SO UP states, M1 SO UP states, and M1 spindles for 90 randomly selected spindle trains. The raster plots are sorted by both the length of the spindle train (in terms of the number of spindles within a given spindle train) and the time from the first to the last spindle peak. Figure 5C specifically plots rasters for spindle trains that contain mPFC–M1 SO pairs within the boundaries of the train (the same legend as Fig. 5B). Quantification of the mPFC SO–M1 SO time (in seconds) is shown as a proportion histogram in Figure 5D (data shown as mean ± SEM for 18 sleep sessions); note the presence of almost exclusively ultra short mPFC–M1 SO locking times thus further suggesting that spindles trains are especially enriched in tight temporal coupling between mPFC and M1 SOs.
Temporally precise interactions between mPFC and M1 SOs during spindle trains. A, Schematic of approach for quantifying mPFC SO to M1 SO timing bracketed by M1 spindle trains. Timing between SOs is measured from the trough of an mPFC SO UP state to the trough of a M1 SO UP state (gray arrow). B, Ninety randomly selected rasters showing timing of mPFC SO UP states (orange ticks), M1 SO UP states (green ticks), and M1 spindle trains (maroon ticks). Rasters are sorted in ascending order by number of spindles inside each spindle train (i.e., “two spindle trains” rasters represent spindle trains containing two spindles and accordingly for three and four spindle trains) and by the duration of the whole spindle train. C, Example randomly selected rasters for spindle trains containing temporally synchronous mPFC–M1 SOs. Same plotting convention as B, each subplot shows trains of different lengths. D, Proportion histogram for time between the UP states of each mPFC SO to M1 SO pair during a spindle train. Data shown as mean ± SEM.
mPFC–M1 SO coupling modulates M1 SO-locked spiking
Prior work has shown a strong burst of local cortical spiking activity during UP states of SOs (Steriade et al., 1993), during reactivations of task-related neural activity (Ramanathan et al., 2015), and during spindles nested within SOs (Silversmith et al., 2020; Lemke et al., 2021). However, no work has previously quantified the impact of cross-area SO–SO coupling on the modulation of cortical spiking depth. Here we analyzed spiking of M1 neurons in response to our three classes of mPFC–M1 SO locking times. Figure 6A shows a representative PETH of M1 single-unit firing activity in response to a local M1 SO (Fig. 6A, SO-triggered average M1 LFP with superimposed single-unit firing activity locked to the UP state of M1 SO; note the transient cortical silencing during the SO DOWN state followed by rebound during the SO UP state).
mPFC–M1 SO proximity modulates SO-locked spike modulation of M1 neurons. A, Example PETH for an M1 unit locked to SOs from the same sleep session (average LFP-triggered SO trace in a gray line; individual spikes in black rasters). B, Average smoothed and Z-scored firing rates locked to the UP state of M1 SOs; each panel presents data for M1 SOs locked to one of our three mPFC SO latency categories (data shown as mean ± SEM for 18 sleep sessions). Black bars highlight regions of interest used to calculate the minimum and maximum Z-scored firing during M1 SO DOWN and UP state, respectively. C, Quantification of DOWN and UP state firing across M1 SOs locked to mPFC SOs of different latencies (data shown as mean ± SEM for 18 sleep sessions).
When quantifying the mean M1 SO-locked firing rate of all M1 neurons across sleep sessions, we found that M1 SOs which follow an ultra short or short latency mPFC SO manifest with the greatest degree of cortical neuron silencing during the SO DOWN state (Fig. 6B, left and middle panels). Conversely, M1 SOs locked to long latency mPFC SOs exhibited the least amount of DOWN state spiking modulation (Fig. 6B, right panel). We next sought to quantify the modulation depth of both UP and DOWN states. Thus, we calculated the minimum DOWN state-locked mean firing rate (from −0.45 to −0.25 s, Fig. 6B black bars below SO DOWN states) as well the maximum UP state-locked mean firing rate (from −0.20 to 0.00 s, Fig. 6B black bars above SO UP states). Quantifying the minimum DOWN state firing rates, we find that the Z-scored firing rate (mean ± SEM) is −0.379 ± 0.009, −0.364 ± 0.007, and −0.329 ± 0.008 for ultra short, short, and long latency M1 SOs (Fig. 6C, data shown as mean ± SEM Z-scored firing rate). Quantifying the maximum UP state firing rate (mean ± SEM), we find that the Z-scored firing is 0.073 ± 0.021, 0.032 ± 0.018, and 0.06 ± 0.02 for ultra short, short, and long latency M1 SOs. A one-way ANOVA revealed a statistically significant main effect of mPFC SO latency on the minimum M1 DOWN state firing rate (F = 9.576; p < 0.001; N = 18 sleep sessions). Tukey post hoc analyses further revealed a statistically significant difference between ultra short and long (p < 0.001) as well as short and long (p = 0.007) but no differences between ultra short and short (p = 0.410) minimum M1 DOWN state firing rate. There were no statistically significant differences in maximum UP state-locked Z–scored firing rates (one-way ANOVA, F = 1.187; p = 0.306; N = 18 sleep sessions). These results thus suggest that mPFC SO proximity to an M1 SO modulates the spiking of M1 neurons.
Spindle trains modulate persistence of reactivations during NREM sleep
SO-spindle nesting has been implicated causally and correlationally in the preservation of memory reactivations linked with task performance (Ramanathan et al., 2015; Gulati et al., 2017; Latchoumane et al., 2017; J. Kim et al., 2019; Schreiner et al., 2021). Conceptually, the goal of reactivation analysis is to quantify the reactivation of task-related neural ensemble activity that occurred during the immediate RTG time period (Fig. 7A, left panel represents a schematic showing single-unit firing activity during the RTG period; Fig. 7C, left heatmap shows representative trial-averaged PETH across single units) with spiking activity that occurs in post-training sleep during train spindles, isolated spindles, and simulated events. The simulated events were spaced after an isolated spindle and used as a time-matched comparison with spindles in a train of the same length. We used GPFA to illustrate the perireach neural activity as a low-dimensional neural trajectory (Yu et al., 2009; J. Kim et al., 2023). We estimated reactivations during spindles by projecting binned spike activity onto the top three reach-related GPFA factors (Fig. 7B, schematic of a neural trajectory reactivation trace).
Spindle trains are associated with persistent memory reactivations. A, Schematic representation of reactivation analysis. Left, Online single-unit firing activity locked to the reach trajectory is used to create a template that is then convolved with off-line firing during spindles. Right, Schematic representation of firing during spindles. We hypothesize train spindles (maroon) as having the strongest reactivation of task-related activity versus isolated spindles (pink) and temporal epochs spaced after isolated spindles (gray). B, Schematic illustrating the calculation of GPFA reactivations. The top three GPFA factors calculated on peri-reach single-unit activity is used to project a concatenated matrix of train and isolated spindle single-unit firing counts onto the three factors. C, Left panel, Trial-averaged heatmap of single-unit Z–scored firing rates. Middle, 3D plot of trial-averaged GPFA trajectory calculated on reach activity (reach template, gray line) and for reactivation activity during a single spindle (spindle reactivation, green line). Note the “compressed” GPFA trajectory seen during spindle reactivations which is consistent with prior literature. Right, Spindle-averaged heatmap of single-unit Z–scored firing rates during spindles. D, Left, Spindle trains preserve the reactivation of neural ensembles during post-training sleep. Dashed line captures linear regression fit across train spindles of length two through five. Right, Isolated spindles followed by equidistant time points (i.e., “randomly delayed events” set 1.5 s apart to match the average time between train spindles) exhibit a decay in neural reactivation. Dashed line captures linear regression line fit across train spindles one through five. The solid line and shaded area represent mean ± SEM for normalized reactivation values (data for 18 sleep sessions). Green lines above plots are representative GPFA trajectories for train, isolated, and randomly delayed spindle events. R, Pearson's correlation coefficient for linear regression fit; p, p value for linear regression fit. Panel B adapted with permission from J. Kim et al. (2023).
We then measured a reactivation strength (i.e., correlation coefficient, Rs) by comparing a neural trajectory during sleep spindles to a template neural trajectory during reach. Figure 7C demonstrates the overlap of neural trajectory calculated using the perireach PETH (Fig. 7C, middle panel) and a neural reactivation trajectory calculated from spindle-associated PETH (Fig. 7C, middle panel). The green lines at the top of Figure 7D show representative neural trajectories quantifying spindle-driven reactivations for train, isolated, and simulated events. We found that train spindles sustained the reactivation of reach-related neural activity across time (Fig. 7D, left panel; linear regression for normalized reactivation strength across spindle trains of up to five spindles in length; Pearson's r = 0.177; p = 0.094; N = 18 sleep sessions). In contrast, if we examine reactivation of neural activity during simulated events following isolated spindles, we found that reactivation of reach-related neural activity exhibits a statistically significant decrease across time (Fig. 7D, right panel; linear regression for normalized reactivation strength across isolated spindles and equidistant points selected to match spindle trains of up to five spindles in length; Pearson's r = −0.249; p = 1.8 × 10−2; N = 18 sleep sessions). Together, these results suggest that spindle trains represent unique epochs of neural activity for sustained task-related memory reactivation.
Discussion
Global synchronization of SOs during learning
SOs have broadly been implicated in the consolidation of new memories. This has been demonstrated both observationally (Stickgold, 2005; Neske, 2016; J. Kim et al., 2023) as well as causally in experiments involving either targeted memory reactivation (TMR) in humans (Marshall et al., 2006; Rasch et al., 2007; Ngo et al., 2013; Ngo and Staresina, 2022) and optogenetic perturbations in animals (Latchoumane et al., 2017; J. Kim et al., 2019). Additional work suggesting SOs benefit memory formation has been shown in studies of aging which suggest that older adults have fewer SOs and with electrophysiological properties that differ from that of younger controls (Helfrich et al., 2019; Bouchard et al., 2021).
Previous work has also shown that SOs originate in frontal cortical areas of both humans (Massimini, 2004) and rodents (Ruiz-Mejias et al., 2011; D. Kim et al., 2015; Niethard et al., 2018; J. Kim et al., 2023), and their origination in frontal areas is followed by (1) propagation of the SO across the cortical mantle posteriorly (Massimini, 2004; Niethard et al., 2018; Oyanedel et al., 2020) and (2) by increased SO–SO coupling in posterior brain regions that are task-relevant (J. Kim et al., 2023). Here we found that SO coincidence across mPFC and M1 exhibited a bimodal timing distribution with an ultra short peak at <100 ms followed by a short peak around ∼1–2 s. It is possible that transient cross-area SO synchronization is driven by subcortical targets with the thalamus being especially well positioned given its capacity to coordinate SO DOWN state transitions through distributed projections to Layer 1 of the cortex (Hay et al., 2021). Another possibility is that cortical somatostatin-positive interneurons in Layers 2/3 provide the necessary bulk population signaling onto a localized neighborhood of apical dendrites of Layer 5 pyramidal neurons, thus synchronizing the phases of SOs across mPFC and M1 in our intracortical LFP recordings (Niethard et al., 2018; Peyrache and Seibt, 2020).
How then can we explain the short latency mPFC–M1 SO interactions? One possibility is that these interactions are a simple co-occurrence; however, our finding of increased train spindle nesting to M1 SOs, including for those M1 SOs antecedent to short latency mPFC SOs, suggests that these mPFC–M1 SO interactions have functional implications. It is possible that short latency mPFC–M1 SO locking occurs through either synchronous subcortical input (e.g., via hippocampal sharp-wave ripples triggering cortical SOs; Peyrache et al., 2009; Staresina et al., 2015; Oyanedel et al., 2020; J. Kim et al., 2023) or a previously reported inherent rhythmicity seen in the SO band (Sanchez-Vives and McCormick, 2000; Ngo et al., 2015). For the latter, it could be that an SO propagating from mPFC leads to localized changes in cortical physiology, such as alterations in dendritic or somatic excitability or in the overall excitation/inhibition balance in the cortex (Niethard et al., 2018), which encourages the generation of a follow-up, localized SO in M1.
Spindle trains preserve memory reactivations in a global context
The unitary nesting of a single SO with a single spindle (i.e., SO-spindle nesting) has previously been linked to the consolidation of motor skill and the sustained preservation of task-related neural ensembles across a night of sleep (Ramanathan et al., 2015; Latchoumane et al., 2017; J. Kim et al., 2019, 2023). Here we extend our understanding of unitary SO-spindle nesting to the concept of spindle trains by showing that (1) train spindles are much more likely to nest with M1 SOs locked to ultra short and short latency mPFC SOs and (2) memory reactivations are preserved across a train of spindles (vs isolated spindles and the time intervals after them). In general, spindles have an inherent rhythmicity to their spectral power envelope, and TMR techniques fail to both solicit the generation of a new spindle and promote off-line gains in task performance when TMR is delivered without consideration for the train phenomenon (Antony et al., 2018). Here we found that train spindles were strongly associated with temporally close mPFC SOs as it relates to both train spindle incidence across time and their nesting to M1 SOs.
It has also been hypothesized that spindle trains prevent interference during memory consolidation (Boutin et al., 2024); others argue that SOs serve as a global coordinator of systems memory consolidation (Peyrache and Seibt, 2020). Our work provides evidence for both perspectives simultaneously: the strong entrainment of train spindles to M1 SOs in close proximity to mPFC SOs is a potential mechanism for the motor system to bind top–down action context with precise kinematic control in a way that minimizes interference from previously learned motor engrams. Moreover, our finding of the preservation of memory reactivations during train spindles, when combined with the fact that spindle peaks are especially enriched in single-unit firing activity (Silversmith et al., 2020; which can lead to spike-timing dependent plasticity; Rosanova, 2005), suggests that spindle trains are well suited to the local consolidation of a recently practiced motor engram. However, future causal experiments are required to fully test this hypothesis. Lastly, it is well established that the learning of a motor task leads to the emergence of neural sequences (Lebedev et al., 2020; Zhou et al., 2020). It is intriguing to speculate that distinct spindles within a spindle train help consolidate the submovements of a given behavior. While our single-unit count is insufficient to dissociate reactivations of submovements, recent advances in neural probes capable of recording hundreds of single units simultaneously (Jun et al., 2017) pave the way for studying if individual spindles within a spindle train subserve the neural substrate for behavioral chunking.
Footnotes
We thank Nikhilesh Natraj and Aviv Mizrahi-Kliger for their critical feedback on the manuscript and input on computational methodology. This work is supported by National Institutes of Health (NIH) 1 F31 NS127514-01 to D.D.; NIH K99NS119737 to J.K.; NIH RF1NS132913, 3R01NS112424, and VHA 5I01RX001640 to K.G.
↵*D.D. and J.K. contributed equally to this work.
The authors declare no competing financial interests.
- Correspondence should be addressed to Karunesh Ganguly at karunesh.ganguly{at}ucsf.edu.
