A spiral attractor network drives locomotion in Aplysia

The neural control of motor behaviour arises from the joint activity of large neuron populations. Unknown is what underlying dynamical system generates this joint activity. Here we show that the network-wide activity driving locomotion of the sea-slug Aplysia is a low-dimensional spiral attractor. We imaged large populations at singlespike resolution from the Aplysia’s pedal ganglion during fictive locomotion. Evoking locomotion rapidly moved the population activity from irregular spontaneous activity into a low-dimensional, periodic, decaying, orbit. This orbit was a true attractor: the activity returned to the same orbit after transient perturbation; and repeatedly evoking locomotion caused the activity to converge on the same low-dimensional orbit. To show this attractor was the locomotion program, we accurately decoded simultaneous recordings of neck motorneuron activity directly from the low-dimensional orbit. Our results provide direct evidence that neural circuits are periodic attractors. Consequently, they support the long-held hypothesis that population activity is an emergent property of a simpler underlying system. Neural control of motor behaviour arises from the mass action of neuron populations1–6. While individual neuron activity can correlate with specific aspects of movement7,8, the intermittent participation of individual neurons across repeated movements9–12 suggests that only the collective population activity, and not specifics of single neuron firing, are key to motor control. Thus, a parsimonious theory of motor control requires unravelling the underlying system embodied by a population13–18. The idea that high-dimensional population activity is an emergent property of an underlying low-dimensional dynamical system has a long history in neural theory19,20. A key hypothesis is that stable repetitions of population activity arise from an underlying, low-dimensional attractor. Stable repetitions of low-dimensional population activity have been described in arm movements15,16, navigation systems21,22, and sensory encoding23,24, but the underlying dynamical system has not been determined. It is thus unknown whether high-dimensional population activity can emerge from a true low-dimensional attractor.

The idea that high-dimensional population activity is an emergent property of an underlying low-dimensional dynamical system has a long history in neural theory 19,20 . A key hypothesis is that stable repetitions of population activity arise from an underlying, low-dimensional attractor. Stable repetitions of low-dimensional population activity have been described in arm movements 15,16 , navigation systems 21,22 , and sensory encoding 23,24 , but the underlying dynamical system has not been determined. It is thus unknown whether high-dimensional population activity can emerge from a true low-dimensional attractor.
To directly test for the presence of an underlying attractor in motor control, we imaged population activity at single-neuron, single-spike resolution in the pedal ganglion of Aplysia during fictive locomotion (Fig 1a). The pedal ganglion presents an ideal target for this problem as it contains the pattern generator 25,26 , motorneurons 27,28 and modulatory neurons 29,30 underlying locomotion. Using this model system, we show here that evoking fictive locomotion caused heterogenous population activity, but under which always lay a low-dimensional, slowly decaying periodic orbit. This periodic trajectory met the convergence and perturbation criteria for an attractor. Crucially, we identify the attractor as a stable, decaying spiral, and show there is this same, single attractor within the pedal ganglion of all tested animals. While individual neurons varied their participation in the attractor between bouts of locomotion, we could nonetheless decode motorneuron activity directly from the low-dimensional orbit.
Our results thus provide experimental support for the long-standing theory that neural population activity is a high-dimensional emergent property of a low-dimensional system.

Results
We sequentially evoked 3 bouts of fictive locomotion in each of 10 isolated central nervous system preparations (Fig 1b). Each bout of locomotion was evoked by short stimulation of the tail nerve P9, mimicking a sensory stimulus to the tail that elicits the escape locomotion response 27 ; in intact animals, a strong tail stimulus typically elicits a twopart escape behavior consisting of several cycles of a vigorous arching gallop, followed by several minutes of a more sedate rhythmic crawl 25,31 . Recorded populations from the pedal ganglion comprised 120-180 neurons each, representing ≈10% of the network in each recording. The population recordings captured rich, varied single neuron dynamics within the ganglion's network following the stimulus (Fig 1c). A dominant, slow (≤ 0.1 Hz) oscillation in neural firing (Fig. 1d) is consistent with the periodic activity necessary to generate rhythmic locomotion. But the variety of single neuron dynamics 32 (Fig 1c) and the slowly decaying population firing rate (Fig. 1e) post-stimulus hint at a more complex underlying dynamical system driving locomotion than a simple, consistent oscillator.
We show here that the joint activity of the population meets the necessary conditions for a periodic attractor (Fig. 1f). We identified these as: (1) applying a driving force causes the system's dynamics to fall onto a stable, periodic orbit; (2) repeatedly driving the system causes convergence of the dynamics to the same orbit; and (3) the system should return to the periodic orbit after the end of transient perturbation. Supplementary  Fig. 1 demonstrates these conditions in a dynamical model.

Joint population activity forms a low-dimensional periodic orbit
We first show that under the heterogenous population activity evoked by the tail-nerve stimulation there is a low dimensional periodic trajectory, consistent with there being a periodic attractor in the pedal ganglion. Projections of a population's joint activity into three dimensions typically showed that stimulation caused a strong deviation from the spontaneous state, which then settled into repeated loops (Fig. 2a). Capturing a significant proportion of the population variance generally required 4-8 embedding dimensions (Fig. 2b), representing a dimension reduction by more than a factor of 10 compared to the number of neurons. However, we cannot directly visualise this space; therefore we cannot tell by visual inspection if the low-dimensional trajectory repeatedly returns to the same position, and so is truly periodic.  To determine whether population activity in higher dimensions reached a stable periodic orbit, we made use of the idea of recurrence 33,34 . For each time-point in the lowdimensional trajectory of the population's activity, we check if the trajectory passes close to the same point in the future (Fig. 2c). If so, then the current time-point recurs, indicating that the joint activity of the population revisits the same state at least once. The time between the current time-point and when it recurs gives us the period of recurrence. A strongly periodic system would thus be revealed by its population's trajectory having many recurrent points with similar recurrence periods; random or chaotic dynamics, by contrast, would not show a single clustered recurrence period.
Plotting recurrent time-points showed that the evoked low-dimensional population activity typically recurred with a regular period (example in Fig. 2d). We found strongly periodic recurrence on the scale of 10-15 s in many but not all of the 30 evoked population responses (Fig. 2e,f). This reflected the range of stimulation responses from strongly periodic activity across the population to noisy, stuttering, irregular activity ( Supplementary  Fig. 2). Nonetheless, despite this heterogeneity across stimulus responses, the activity of almost all populations was dominated by a single periodic orbit (Fig. 2e), robust to the choice of threshold for defining recurrence ( Supplementary Fig. 3).

Joint population activity meets the conditions for a periodic attractor
The trajectory of the population activity describes a circumscribed region of low-dimensional space -the manifold -within which it remains. To truly be an attractor ( Fig. 1f; Supplementary Fig. 1), the system's dynamics should rapidly reach and stay on the manifold when evoked; reach the same manifold every time it is evoked; and return to the manifold when perturbed. We found that almost all evoked population responses quickly reached the attractor, within one oscillation period (Fig. 3a), and were then dominated by recurrence, indicating they stayed on the manifold.
The key problem in analysing any putative attractor from experimental data is identifying when the experimentally-measured dynamics are or are not on the attractor, whether due to perturbations of the system or noise in the measurements. The set of recurrent time-points allowed us to identify when the joint population activity was most likely on the attractor.
To determine if sequentially-evoked responses from the same preparation reached the same manifold, we projected all 3 population responses into the same set of embedding dimensions, using only the recurrent points ( Fig. 3b; Supplementary Fig. 4 shows these results are robust to other projections). Falling on the same manifold would mean that every recurrent point in one population would also appear in both others, if noiseless. Consequently, the maximum distance between any randomly chosen recurrent point in population response A and the closest recurrent point in population response B should be small. We defined small here as being shorter than the expected distance between a recurrent point in A and the closest point on a random walk of the activity in the same embedding dimensions. Despite the inherent noise and limited duration of the recordings, this is exactly what we found: pairs of evoked population responses from the same preparation fell close to each other throughout (Fig. 3c), well in excess of the expected agreement between random projections of the data onto the same embedding dimensions.
We also checked that this convergence to the same manifold came from different initial conditions. The initiating stimulation is a rough kick to the system -indeed a fictive locomotion bout can be initiated with a variety of stimulation parameters 32 -applied to ongoing spontaneous activity. Together, the stimulation and the state of spontaneous activity when it is applied should give different initial conditions from which the locomotion  manifold is reached. We found that the stimulation caused population responses within the same preparation to diverge far more than in either the spontaneous activity or after coalescing to the manifold (Fig. 3d). Thus, a wide range of initial driven dynamics in the pedal ganglion population converged onto the same manifold. Previous studies have used the consistency of pairwise correlations between neurons across conditions as indirect evidence for the convergence of population activity to an underlying attractor 22,35 . The intuition here is that neurons whose activity contributes to the same portion of the manifold will have simultaneous spiking, and so their activity will correlate across repeated visits to the manifold. To check this, we computed the pairwise similarity between all neurons within an evoked population response (Fig. 3e), then correlated these similarity matrices between responses from the same preparation. We found that pair-wise similarity is indeed well-preserved across population resposnes in the same preparation (Fig. 3f). This also shows that the apparent convergence to the same manifold is not an artifact of our choice of low-dimensional projection.
In many population responses, we noticed spontaneous perturbations of the lowdimensional dynamics away from the trajectory (examples in Supplementary Fig. 5), indicated by sudden falls in the density of recurrent points (Fig. 3g). That is, perturbations could be detected by runs of contiguous points on the population trajectory that were not recurrent. In most cases (90%), the population dynamics returned to a recurrent state after the spontaneous perturbation (Fig. 3h); the two perturbations that did not were consistent with the end of the evoked fictive locomotion and a return to spontaneous activity ( Supplementary Fig. 5). Of those that returned, all but three clearly returned to the same manifold (Fig. 3i); for those three, the spontaneous perturbation appeared sufficient to move the population dynamics into a different periodic attractor ( Supplementary  Fig. 5). Potentially, these are the known transitions from the escape gallop to normal crawling 31 . The low dimensional dynamics of the pedal ganglion thus meet the stability, manifold convergence, and perturbation criteria of a periodic attractor network.

Heterogenous population activity arises from a common attractor
While these results show the existence of a periodic orbit on an attractor in the evoked population responses, they cannot address whether these arise from the same putative attractor within and, crucially, between animals. To determine if there is a common underlying attractor despite the heterogeneity in spiking patterns across the population responses ( Supplementary Fig. 2), we introduced a statistical approach to quantifying the low-dimensional trajectory. We fitted a linear model of the local dynamics around each time point (see Methods), which allowed us to then estimate the local dynamics based on the eigenvalues of that model: whether the trajectory is expanding away or contracting towards this point, and whether or not it is rotating.
The time-series of eigenvalues typically showed long periods of similar magnitude eigenvalues, corresponding to the recurrent points (Fig 4a). Consequently, by averaging over the local dynamics only for recurrent points, we could potentially capture the dynamics of the underlying attractor. Doing so, we found that the evoked population responses had highly clustered eigenvalues (Fig. 4b,c), and thus highly similar underlying dynamics despite the apparent heterogeneity of spike-train patterns between them. The dominance of negative complex eigenvalues implies the pedal ganglion network implements a low-dimensional system containing a stable spiral -a decaying periodic orbit.
The identification of a stable spiral makes a clear prediction for what should and should not change over time in the dynamics of the population. The negative complex eigenvalues mean that the magnitude of the orbit decays over time, corresponding to the  Fig. 4). d Distances between pairs of population responses from the same preparation in three states: the end of spontaneous activity (at stimulus onset); between stimulation onset and coalescence (the maximum distance between the pair); and after both had coalesced (both reaching the putative attractor manifold; data from panel c). e Example neuron activity similarity matrices for consecutively evoked population responses. Neurons are ordered according to their total similarity in stimulation 2. f Correlation between pairs of neuron similarity matrices (Data) compared to the expected correlation between pairs of matrices with the same total similarity per neuron (Control). Values below the diagonal indicate conserved pairwise correlations between pairs of population responses within the same preparation. The two pairs on the diagonal are response pairs (1,3) and (2,3) in preparation 7; this correctly identified the unique presence of a random walk in response 3 ( Supplementary Fig. 4). g Spontaneous divergence from the trajectory. For one population response, here we plot the density of recurrence points (top) and the mean recurrence delay in 5s sliding windows. Coalescence time: grey line. The sustained "divergent" period of low recurrence (grey shading) shows the population spontaneously diverged from its ongoing trajectory, before returning. Black dot: pre-divergence window (panel i). h Breakdown of spontaneous perturbations across all population responses. Returned: population activity became stably recurrent after the perturbation. i Returning to the same manifold. Proportion of the recurrent points in the pre-divergence window that recurred after the divergent period, indicating a return to the same manifold or to a different manifold.  decreasing population spike rate in most evoked responses (Fig. 1e). However, a stable spiral indicates only a decrease in magnitude; it does not mean the orbital period is also slowing. Consequently, the presence of a stable spiral attractor predicts that the magnitude and period of the orbit are dissociable properties in the pedal ganglion network.
We checked this prediction using the linear model. The linear model estimated a mean orbital period of around 10 s (Fig. 4c), consistent with the directly-derived estimate from the recurrent points (Fig. 2f). This indicated the linear model was correctly capturing the local dynamics of each program. But our linear model also gave us a time-series of estimates of the local orbital period (Fig 4d), which we could use to check whether the orbital period was changing during each evoked response. We found that the population responses included all possible changes in periodic orbit: slowing, speeding up, and not changing (Fig 4e). As predicted there was no relationship between the contraction of the periodic orbit and its change in period (Fig 4f).
The locomotion motor program can be decoded from the low-dimensional orbit.
Collectively, these periodic, decaying dynamics are ethologically consistent with locomotion that comprises a repeated sequence of movements that decays in intensity over time 31,36,37 . If this putative low-dimensional periodic attractor is the "motor program" for locomotion, then we should be able to decode the locomotion muscle commands from its trajectory. In 3 of the 10 preparations we were able to simultaneously record activity from the P10 nerve that projects to the neck muscles 38 for all three evoked population responses. The spiking of axons in this nerve should thus correspond to the specific neck contraction portion of the cyclical escape locomotion response. We thus sought to decode the spiking of P10 directly from the low-dimensional population trajectory (Fig. 5a).
We first confirmed that each recorded neural population did not appear to contain any motorneurons with axons in P10, which could make the decoding potentially trivial ( Supplementary Fig. 6). To then decode P10 activity, we used a statistical model that predicts the firing rate of nerve P10 at each time point, by weighting and summing the recent history (up to 100 ms) of the trajectory in the low dimensional space, and using a non-linearity to convert this weighted sum into a firing rate. We controlled for over-fitting using cross-validation forecasting: we fit the model to a 40 s window of trajectory data, and predicted the next 10 s of P10 activity (Fig. 5b). By sliding the window over the data, we could assess the quality of the forecast over the entire recording (Fig. 5c).
The model could accurately fit ( Supplementary Fig. 7) and forecast (Fig. 5b) P10 activity from the low-dimensional trajectory in all 9 population responses. Emphasising the quality of the model, in Fig. 5d we plot example forecasts of the entire P10 recording based on fitting only to the first 40s window, each example taken from the extremes we obtained for the fit-quality metrics. Notably, in one recording the population response shutdown half-way through; yet despite the model being fit only to the 40s window containing strong oscillations, it correctly forecasts the collapse of P10 activity, and its slight rise in firing rate thereafter. Within the limits of this sample of data, the low dimensional trajectory of the periodic attractor appears to encode the motor program.
To confirm this, we asked whether the motor program -as represented by the P10 activity, -was truly low-dimensional. The successful decoding of future P10 activity was achieved despite needing only 3-5 embedding dimensions to account for 80% variance in the population activity for these nine recordings ( Supplementary Fig. 8

Variable neuron participation in stable motor programs
If the low-dimensional trajectory described by the joint activity of the population just is the locomotion motor program, then how crucial to this program are the firing of individual neurons [9][10][11][12]39 ? Having quantified the motor program as the low-dimensional activity trajectory, we could uniquely ask how much each neuron participated in each evoked program. We quantified each neuron's participation as the absolute sum of its weights on the principal axes (eigenvectors): large total weights indicate a dominant contribution to the low-dimensional trajectory, and small weights indicate little contribution. So quantified, participation includes both firing rate and synchrony effects on each neuron's contribution ( Supplementary Fig. 9). Every population response had a long-tailed distribution of participation (Fig. 6a), indicating that a minority of neurons dominated the dynamics of any given response. Nonetheless, these neurons were not fixed: many with high participation in one population response showed low participation in another (Fig. 6b,c). To rule out noise effects on the variability of participation (for example, due to the finite duration of recording), we fitted a noise model to the change in participation separately for each preparation (Fig. 6d,e). Every preparation's pedal ganglion contained neurons whose change in participation between responses well-exceeded that predicted by the noise model (Fig. 6f). Consequently, the contribution of single neurons was consistently and strongly variable between population responses in the same preparation.
These data show that a neuron's role within the locomotion motor program is not fixed, but leave open the question of whether single neuron variability causes variation in the program itself. In our analysis, variation between sequentially-evoked population responses is quantified by the distance between their low-dimensional projections (as in Fig. 3c). We found that the distance between a pair of population responses did not correlate with either the total change in neuron participation between the two responses ( Fig. 6g) or the distance between their participation distributions (Fig. 6h). The execution of the motor program is thus robust to the participation of individual neurons.

Participation maps identify potential locations of the pattern generator network
To get some insight into the physical substrate of the attractor, we plotted maps of the participation of each neuron in each preparation. We found that neurons with strong participation across the three evoked population responses were robustly located in the caudo-lateral quadrant of the ganglion (Fig. 7a,b). Maps of the right ganglion also indicated strong participation in the rostro-medial quadrant; due to the low numbers of maps for each side, it is unclear whether this is a true asymmetry of the ganglia or simply reflects sampling variation. Neurons with highly variable participation between population responses (Fig. 7c,d) were similarly found in the caudo-lateral quadrants of both ganglia. Strongly participating neurons were thus confined to specific distributed regions of the pedal ganglion's network.
These data are consistent with a network-level distribution of the attractor, with a particularly strong contribution from the caudo-lateral quadrant. Encouragingly, from a different data-set we previously described this region as containing neural ensembles that generated a cyclical packet of neural activity, which moved in phase with activity from the Neurons beyond thresholds (grey lines) of mean ±3SD of the fitted model were identified as strongly variable. f Proportion of identified strongly variable neurons per preparation. g Distance between pairs of population responses as a function of the total change in neuron participation between them. Each dot is a pair of responses from one preparation; the distance between them is given relative to the distance between the random projections, allowing comparison between preparations (Fig. 3c). Black dots are excluded outliers, corresponding to the pairs containing response 1 in preparation 4 with apparent chaotic activity ( Supplementary Fig. 4). h Distance between pairs of population responses as a function of the distance between the distributions of participation (panel a). Conventions as for panel g. neck-projecting P10 nerve 32 . Consequently, both those data and our new data support our hypothesis that the pattern generator for locomotion is predominantly located in the caudo-lateral network.

Discussion
Locomotion networks provide a tractable basis for testing theories of neural dynamics 6,13,14,32,[40][41][42] , as they couple complex dynamics with clearly defined outputs. We took advantage of this to test the hypothesis that neural circuits form attractor networks. Across all 30 evoked population responses examined here, there was a remarkable heterogeneity of spike-train patterns, from visually evident widespread oscillations to noisy, stuttering oscillations in a minority of neurons. Yet our analysis shows that underpinning this heterogeneity is the same dynamical system: a low-dimensional, decaying, periodic orbit. We found a remarkably consistent periodicity and rate of orbital decay across evoked responses within a preparation and between preparations. The stability of these dynamics, and the convergence of population activity to the same manifold, suggests our recordings have captured for the first time a stable-spiral attractor neuronal network.
Treating a neural circuit as a realisation of a dynamical system takes the emphasis away from the details of individual neurons -their neurotransmitters, their ion channel repertoire -and places it instead on their collective action. This allows us to take a Marr-ian perspective 43 , which neatly separates the computational, algorithmic, and implementation levels of motor control. The computational problem here is of how to generate rhythmic locomotion for a finite duration; the algorithmic solution is a decaying periodic attractor; and the implementation of that attractor is the particular configuration of neurons in the pedal ganglion -one of many possible implementations [44][45][46][47] .
We saw the separation of these levels most clearly in the variable participation of the individual neurons between evoked bouts of fictive locomotion. The projection of the pedal ganglion network's joint activity into a low dimensional space captured the locomotion motor program independently of any single neuron's activity. Even the most strongly participating neurons in a given population response could more than halve their participation in other evoked responses. These results suggest that the pedal ganglion's pattern generator is not driven by neurons that are endogenous oscillators, as they would be expected to participate equally in every response. Rather, this variation supports the hypothesis that the periodic activity is an emergent property of the network.
The wide variation of single neuron participation between evoked bouts of fictive locomotion also raises new challenges for theories of neural network attractors. While a variety of models present solutions for self-sustaining periodic activity in a network of neurons 44,46,47 , it is unclear if they can account for the variable participation of single neurons. A further challenge is that while the variable participation of individual neurons does not affect the underlying program, clearly it takes a collective change in single neuron activity to transition between behaviours -as, for example, in the transition from galloping to crawling in Aplysia. What controls these transitions, and how they are realised by the population dynamics, is yet to be explored either experimentally or theoretically.
Our results nonetheless argue against a number of hypotheses for the implementation of rhythmic locomotion by the pedal ganglion. As noted above, such single neuron variability between sequential locomotion bouts argues against the generation of rhythmic activity by one or more independent neurons that are endogenous oscillators. Our results also argue against the existence of many stable periodic states in this network 45 . Such meta-stability would manifest as changes in periodicity following perturbation. Our results show that spontaneous divergences from the attractor overwhelmingly returned to the same attractor. We observed only three (of 17) that were consistent with moving to a different periodic state; potentially these were transitions from the escape gallop to normal crawling 31 .
How then might the pedal ganglion network implement a periodic attractor? Our data were collected from an isolated central nervous system preparation, in which the modulatory influence of neurons outside the pedal ganglion cannot be discounted 48 . Nonetheless, as the pedal ganglion contains the central pattern generator for locomotion 26 , we can suggest how that generator is realised. Our results here support the hypothesis that the periodic activity is an emergent property of the ganglion's network. We know the pedal ganglion contains a mix of interneurons and motorneurons 28 , and that the motorneurons are not synaptically coupled 27 , suggesting they read-out (and potentially feedback to) the dynamics of an interneuron network. An hypothesis consistent with our results here is that the ganglion contains a recurrent network of excitatory interneurons, centred on the caudo-lateral quadrant, which feed-forward to groups of motorneurons 32 . This recurrent network embodies the attractor, in that stimulation of the network causes a self-sustained packet of activity to sweep around it 32 . We see this as the periodic trajectory of joint population activity (cf Fig. 2a, Fig. 3b).
Our data further suggest that the pedal ganglion network supports at least two stable states, the spontaneous activity and the stable-spiral attractor. Reaching the stable-spiral attractor from the spontaneous activity required long-duration, high-voltage pedal nerve stimulation (Fig 1; ref. 32 ). In dynamical systems terms, this suggests that the spontaneous state's basin of attraction is large: most perturbations return to that state, and it takes a large perturbation to move into a different basin of attraction. The co-existence of stable spontaneous and periodic states in the same network suggests that something must reconfigure the network to sustain periodic activity; otherwise, irrespective of the stimulation, the network would always return to the spontaneous state. One possibility in the pedal ganglion is that serotonin alters the effective connections between neurons: escape galloping is both dramatically extended by systemic injection of serotonin alongside tail stimulation 37 , and evoked by stimulating serotonergic command neurons CC9/CC10 in the cerebral ganglion 48 . Future experimental work should thus test the stability of the spontaneous state, and test how manipulating serotonin affects reaching and sustaining the stable-spiral attractor.
Finding and quantifying the attractor required new analytical approaches. We introduce here the idea of using recurrence analysis to solve two problems: how to identify rotations of activity in a high-dimensional space (not necessarily periodic); and how to identify when the recorded system is and is not on the putative underlying attractor. By extracting the times when the population activity is on that attractor, we could then quantify and characterise it, including identifying transient perturbations, and estimating changes in orbital period. Crucially, these extracted attractor times let us further introduce the idea of using linear models as a statistical estimator, to identify the type of attractor, and compare the detected attractor's parameters within and between preparations. Our analysis approach thus offers a road-map for further studies of population dynamics.
The grand challenge of systems neuroscience is to find simplifying principles in the face of overwhelming complexity 5,49 . We have provided evidence here that neural population activity is a high-dimensional view of the true underlying low-dimensional dynamical system. It follows that dimension reduction approaches are not just useful descriptions of population activity 17 , but can reveal the underlying system realised by a neuronal network 18,32 .

Online Methods
Imaging Full details of the preparation are given in ref. 32 . Briefly, the cerebral, pleural and pedal ganglia were dissected out, pinned to the bottom of a chamber, and maintained at 15 − 17 • C. Imaging of neural activity used the fast voltage sensitive absorbance dye RH-155 (Anaspec), and a 464-element photodiode array (NeuroPDA-III, RedShirtImaging) sampled at 1600 Hz. Optical data from the 464 elements were bandpass filtered in Neuroplex (5 Hz high pass and 100 Hz low pass Butterworth filters), and then spike-sorted with independent component analysis in MATLAB to yield single neuron action potential traces (the independent components), as detailed in 50 . Rhythmic locomotion motor programs were elicited using 8V 5ms monophasic pulses delivered at 20Hz for 2.5s via suction electrode to pedal nerve 9. A separate suction electrode was attached to pedal nerve 10 to continuously monitor the locomotion rhythm 38 . The stimulation protocol (Fig. 1b) used short (15 mins) and long (60 mins) intervals between stimulations, as the original design also sought effects of sensitisation.
Spike-train analysis Power spectra were computed using multi-taper spectra routines from the Chronux toolbox 51 . We computed the power spectrum of each neuron's spiketrain post-stimulation, and plot means over all spectra within a recorded population, and the mean over all mean spectra. We computed the spike-density function f (t) for each neuron by convolving each spike with a Gaussian, with σ = {median ISI in population}/ √ 12, using a time-step of 10 ms. To visualise the entire population's spiking activity (Fig. 1c), we cluster neurons by the similarity of their firing patterns using our modular deconstruction toolbox 32 . Different dynamical types of ensembles were identified by the properties of their autocorrelograms -see 32 for details.
Model network We used a three-neuron network to demonstrate the dynamical properties of a periodic attractor as realised by neurons (Supplementary Fig. 1). Each neuron's membrane dynamics were given by τ aȧi = −a i (t) + c i (t) + 3 j=1 w ji r j (t) − γy i (t), with adaptation dynamics τ yẏi = −y i (t) + r i (t), and output firing rate r i (t) = max {0, a i (t)}. Weights w ji ≤ 0 give the strength of inhibitory connections between the neurons, each of which receives a driving input c i . This model, due to Matsuoka 52,53 , generates selfsustained oscillation of network firing rates given constant scalar inputs c i (t) = c, despite each neuron not being an endogenous oscillator. The time constants of membrane τ a and adaptation τ y dynamics, together with the strength of adaptation γ, determine the periodicity of the oscillations 52,53 . Here we use τ a = 0.025, τ y = 0.2, and γ = 2; input was c i = 3 throughout except where noted.
Recurrence analysis Low dimensional projections of the joint population activity were obtained for each program using standard principal components analysis, applied to the covariance matrix of the spike-density functions. The d leading eigenvectors W i of the covariance matrix define the d principal dimensions, and the d corresponding eigenvalues are proportional to the variance accounted for by each dimension. Projections (the "principal components") onto each of the chosen dimensions is given by where the sum is taken over all n neurons in the analyzed population.
We used recurrence analysis 33,34 to determine if the low-dimensional projection contained a stable periodic orbit. To do so, we checked if the low-dimensional projection P (t) = (p 1 (t), p 2 (t), . . . , p d (t)) at time t recurred at some time t + δ in the future. Recurrence was defined as the first point P (t + δ) = (p 1 (t + δ), p 2 (t + δ), . . . , p d (t + δ)) that was less than some Euclidean distance θ from P (t). The recurrence time of point P (t) is thus δs. Contiguous regions of the projection's trajectory from P (t) that remained within distance θ were excluded. Threshold θ was chosen based on the distribution of all distances between time-points, so that it was scaled to the activity levels in that particular program. Throughout we use the 10% value of that distribution as θ for robustness to noise; similar periodicity of recurrence was maintained at all tested thresholds from 2% upwards (Supplementary Figure 3).
We checked every time-point t between 5s after stimulation until 10s before the end of the recording (around 7770 points per program), determining whether it was or was not recurrent. We then constructed a histogram of the recurrence times using 1s bins to detect periodic orbits (Fig. 2e): a large peak in the histogram indicates a high frequency of the same delay between recurrent points, and thus a periodic orbit in the system. All delays less than 5s were excluded to eliminate quasi-periodic activity due to noise in otherwise contiguous trajectories. Peaks were then defined as contiguous parts of the histogram between empty bins, and which contained more than 100 recurrent points. Programs had between one and four such periodic orbits. The peak containing the greatest number of recurrent points was considered the dominant periodic orbit of the program; the majority of programs had more than 50% of their recurrent points in this peak (blue-scale vectors in Fig. 2e). The mean orbit period of the program was then estimated from the mean value of all recurrence times in that peak.
We measured the attractor's stability as the percentage of all points that were in periodic orbits. Evolving dynamics of each program were analysed using 5 s sliding windows, advanced in steps of 1 s. We defined the "coalescence" time of the attractor as the mid-point of the first window in which at least 90% of the points on the trajectory were recurrent.
Testing convergence to the same manifold To determine if sequentially-evoked programs had the same manifold, we determined how closely the trajectories of each pair of programs overlapped in the low-dimensional space. We first projected all three programs from one preparation onto the principal axes of first program, to define a common lowdimensional space. The shuffling was repeated 100 times. Mean ± 2SEM of the shuffled distances are plotted in (Fig. 3c); the error bars are too small to see.
To check the robustness of the convergence to the same manifold, we repeated this analysis starting from a common set of principal axes for the three programs, obtained using principal component analysis of their concatenated spike-density functions. We plot the results of this analysis in Supplementary Fig. 4a.
As a further robustness control, we sought evidence of the manifold convergence independent of any low-dimensional projection. We made use of the idea that if neurons are part of sequential programs on a single manifold, then the firing of pairs of neurons should have a similar time-dependence between programs 22,35 . For each pair of programs (A, B) from the same preparation, we computed the similarity matrix 54 S(A) between the spike-density functions of all neuron pairs in A, and similarly for B, giving S(B). We then computed the correlation coefficient between S(A) and S(B): if A and B are on the same manifold, so their pairwise correlations should themselves be strongly correlated. As a control we computed a null model where each neuron has same total amount of similarity as in the data, but its pairwise similarity with each neuron is randomly distributed 54 . The expected value of pairwise correlation between neurons i and j under this model is then E ij = s i s j /T , where (s i , s j ) are the total similarities for neurons i and j, and T is the total similarity in the data matrix. For comparison, we correlated S(A) with E, and plot these as the control correlations in Fig. 3e.
Testing return to the same manifold after perturbation We detected divergences of the trajectory away from the putative manifold, indicating spontaneous perturbations of population dynamics. We first defined potential perturbations after coalescence as a contiguous set of 5s windows when the density of recurrent points was below 90% and fell below 50% at least once. The window with the lowest recurrence density in this divergent period was labelled the divergent point. We removed all such divergent periods whose divergent point fell within 2 oscillation cycles of the end of the recording, to rule out a fall in recurrence due solely to the finite time horizon of the recording. For the remaining 19 divergent periods, we then determined if the population activity returned to a recurrent state after the divergent point; that is, whether the density of recurrence returned above 90% or not. The majority (17/19) did, indicating the perturbation returned to a manifold.
For those 17 that did, we then determined if the recurrent state post-divergence was the same manifold, or a different manifold. For it to be the same manifold after the spontaneous perturbation, then the trajectory before the perturbation should recur after the maximum divergence. To check this, we took the final window before the divergent period, and counted the proportion of its recurrent delays that were beyond the end of the divergent period, so indicating that the dynamics were in the same trajectory before and after the divergence. We plot this in Fig. 3h.
Statistical estimation of the attractor's parameters We introduce here a statistical approach to analysing the dynamics of low-dimensional projections of neural activity timeseries obtained from experiments. We first fitted a linear model around each point on the low-dimensional trajectory to capture the local dynamics. For each point P (t), we took the time-series of points before and after P (t) that were contiguous in time and within 2.5 × θ as its local neighbourhood; if less than 100 points met these criteria P (t) was discarded. We then fitted the dynamical modelṖ * = AP * that described the local evolution of the low-dimensional projection P * by using linear regression to find the Jacobian matrix A; to do so, we use the selected local neighbourhood time-series as P * , and their first-order difference asṖ * . The maximum eigenvalue λ = a + ib of A indicates the dominant local dynamics 55 , whether contracting or expanding (sign of the real part a of the eigenvalue), and whether oscillating or not (existence of the complex part of the eigenvalue i.e. b = 0). We fitted A to every point P (t) after the stimulation off-set, typically giving ≈ 9000 local estimates of dynamics. The dominant dynamics for the whole program were estimated by averaging over the real a and the complex b parts of the maximum eigenvalues of the models fitted to all recurrent points in the dominant periodic orbit. The linear model's estimate of the orbit rotation period was estimated from the complex part as ω = 2πb∆t, with the sampling time-step ∆t = 0.01s here. The linear model's estimate of the contraction rate is exp(a/∆t), which we express as a percentage.
Tracking changes in periodicity over a program We tracked changes in the oscillation period by first averaging the recurrence time of all recurrent points in a 5s sliding window. We then correlated the mean time with the time-point of the window to look for sustained changes in the mean period over time, considering only windows between coalescence and the final window with 90% recurrent points. We used a weighted version of Spearman's rank to weight the correlation in favour of time windows in which the trajectory was most clearly on the periodic orbit, namely those with a high proportion of recurrent points and low variation in recurrence time. The weighted rank correlation is: given vectors x and y of data rankings, and a vector of weights w, compute the weighted mean m = i w i x i / i w i and standard deviation σ xy = i w i (x i − m x )(y i − m y )/ i w i , and then the correlation ρ = σ xy / √ σ xx σ yy . We used the weight vector: w i = s −1 i Q i , where s i is the standard deviation of recurrence times in window i, and Q i is the proportion of recurrent points in window i. P-values were obtained using a permutation test with 10000 permutations.
Decoding motor output We decoded P10 activity from the low-dimensional trajectory of population activity using a generalised linear model. We first ruled out that any simultaneously recorded neuron was a motorneuron with an axon in nerve P10, by checking if any neurons had a high ratio of locking between their emitted spikes and spikes occurring at short latency in the P10 recording. Supplementary Fig. 6 shows that no neuron had a consistent, high ratio locking of its spikes with the P10 activity.
We convolved the spikes of the P10 recording with a Gaussian of the same width as the spike-density functions of the simultaneously recorded program, to estimate its continuous firing rate f 10 . We fitted the model f 10 to determine the P10 firing rate as a function of the past history of the population activity trajectory. Using a generalised linear model here allows us to transform the arbitrary coordinates of the d-dimensional projection P (t) into a strictly positive firing rate. Fitting used glmfit in MATLAB R2014. To cross-validate the model, we found the coefficients β using a 40s window of data, then forecast the P10 firing rate f * 10 using the next 10 seconds of population recording data as input to the model. Forecast error was measured as both the median absolute error and the correlation coefficient R between the actual and forecast P10 activity in the 10s window. The fitting and forecasting were repeated using a 1s step of the windows, until the final 40s+10s pair of windows available in the recording.
We tested activity histories between 50 and 200ms duration, with time-steps of 10ms, so that the largest model for a given program had d × 20 coefficients. These short windows were chosen to rule out the contributions of other potential motorneurons in the population recording that would be phase offset from neck contraction (as 200 ms is 2% of the typical period). All results were robust to the choice of history duration, so we plot results using history durations that had the smallest median absolute error in forecasting for that program.
Single neuron participation We quantified each neuron's participation in the lowdimensional projection as the L1-norm: the absolute sum of its weights on the principal axes (eigenvectors) for program m: ρ m i = d j=1 |λ m j W m j (i)|, where the sum is over the d principal axes, W m j (i) is the neuron's weight on the jth axis, and λ m j is the axis' corresponding eigenvalue. Within a program, participation for each neuron was normalised to the maximum participation in that program. To fit a noise model for the variability in participation between programs, we first computed the change in participation for each neuron between all pairs of programs in the same preparation. We then fit a Gaussian model for the noise, using an iterative maximum likelihood approach to identify the likely outliers; here the outliers are the participation changes that are inconsistent with stochastic noise. In this approach, we compute the mean and variance of the Gaussian from the data, eliminate the data-point furthest from the estimate of the mean, re-estimate the mean and variance, and compute the new log likelihood of the Gaussian model without that data-point. We iterate elimination, re-estimation, and likelihood computation until the likelihood decreases. The final model (mean and variance) found before the decrease is then the best-fit Gaussian model to the bulk of the data. Neurons whose maximum change in participation exceeded a threshold of the mean ±3SD of that best-fit model were then considered "strongly variable" neurons.
We asked whether the variation in low-dimensional dynamics of sequentially-evoked programs was a consequence of the degree of variation in single neuron participation. Between a pair of consecutively evoked programs, we quantified the variation in their low dimensional dynamics as the Hausdorff distance between them, normalised by the mean distance between their random projections. This normalisation allowed us to put all programs on a single scale measuring the closeness relative to random projections, such that 1 indicates equivalence to a random projection, < 1 indicates closer than random projections, and > 1 indicates further apart then random projections. For a given pair of programs, we quantified the variability of individual neurons' participation in two ways: by summing the change in participation of each neuron between the programs; and by computing the Hellinger distance between the two distributions of participation (one distribution per program).
Participation maps Each neuron's (x,y) location on the plane of the photodiode array could be estimated from the weight matrix from the independent component analysis of the original 464 photodiode time-series; see 32 for full details. We were able to reconstruct locations for all neurons in 8 of the 10 recorded preparations; for the other two preparations, partial corruption of the original spike-sorting analysis data prevented reconstructions of their neuron locations. We merged all left or right ganglion recordings on to a common template of the photodiode array. The marker sizes and colour codes for each neuron were proportional to the normalised maximum participation of that neuron (Fig. 7a,c) and to the range of normalised maximum participation across the three programs (Fig. 7b,d). Here we demonstrate the essential properties of a periodic attractor network. a A simple neuronal implementation of a periodic attractor network using three mutually inhibitory neurons. Model neurons are not endogenous bursters: oscillatory activity arises from the network. b Underlying periodic attractor is revealed by the recurrence of the neurons' joint activity. Top: given sustained input, the three neurons rapidly settle into an alternating oscillatory rhythm. Bottom: projection of the three neurons' firing rates onto the first two principal components (PCs). Starting from time zero (black dot), the joint activity of the network rapidly converges to a recurrent trajectory, where each corner of the triangle -extending into the short diagonal line -corresponds to the maximum rate of one neuron. The space containing the repeated trajectory is the attractor's manifold. c Convergence to the manifold from different initial conditions. From three different initial conditions for the neurons (coloured dots), each corresponding joint activity converges to the attractor manifold. d Transient perturbations cause a divergence and return to the attractor's manifold. Top: A brief increase in input to one neuron (black bar) perturbs the ongoing oscillation, which is rapidly recovered. Bottom: in the PC projection, the initial recurrent trajectory (in black) is moved away from the manifold by the perturbation (red); removing the perturbing input rapidly moves the joint activity back to the manifold (blue). e Sustained perturbations move the network to a new attractor. Top: A sustained change in input (black bar) causes one neuron (purple) to become silent and others to change their oscillation patterns. Bottom: in the PC projection, the perturbation causes the joint activity of the network to move from the initial stable manifold (black) to a new stable manifold (blue: note this is a repeating line). f Decaying input causes a spiral. Top: linearly decaying input to all neurons changes the amplitude, but not the frequency of the oscillatory activity. Bottom: in the PC projection, this decay is shown by the spiral of the joint activity's trajectory to the stable point at (0,0) -which corresponds to no activity.  Figure S2: Range of dynamics in the evoked locomotion programs. Here we quantify the complexity of the dynamics in each evoked program using the structure of the recurrence plots (5 example plots are shown). Each white dot is a pair of time-points at which the low-dimensional trajectory recurred. A diagonal line indicates a contiguous period of recurrence of the population dynamics. We thus quantified the complexity of each program by the duration of the longest diagonal line (measuring the predictability of the system) and by the proportion of recurrent points falling in long (> 5 s) diagonal lines (measuring the cleanness of the system's dynamics). We plot the joint distribution of these two quantities here, one point per program; colours indicate programs from the same preparation. These show how the range of dynamics captured by our dataset stretches from quasi-noise (bottom left), though stuttering/pausing oscillations (top left), to continuous oscillation (right). Despite this heterogeneity, almost all evoked programs revealed the same underlying dynamical system (Figs. 2-4, main text) .  Figure S4: Convergence and non-convergence to the same manifold. a Robustness of convergence to the same manifold. Here we repeat the analysis of Figure 3c, but using different low-dimensional projections. In the main text Figures 2-4, we analyse low-dimensional projections of each individual program; the comparison of distances between sequentially evoked programs was then done by projecting all three programs onto the axes of the first program. Here we derive a common set of axes, by first applying PCA to the concatenated timeseries of the three programs. We then project each program onto those common axes, before computing distances between data and control (random projection) pairs of programs (dots). We obtain the same results: sequentiallyevoked programs converge to the same manifold. b Apparent chaotic motion in the non-convergent program.
Here we show the recurrence plot for program 1 of preparation 4, the only program that was further from the other programs in the same preparation than random projections. We see that this program seemed to reach an apparently chaotic attractor about 45s after stimulation: quasi-periodic activity, with no regular structure. c Apparent random-walk in the non-correlated program. Here we show the recurrence plot for program 3 of preparation 7, the only program that did not recapitulate the pair-wise correlations of neural activity of the other programs in the same preparation. We see that this program had no clear periodic activity, but seemingly random periods of locally-similar activity, consistent with a dynamical system undergoing a random walk. a b c Figure S5: Spontaneous perturbations of ongoing programs. We plot here examples of each type of spontaneous perturbation we observed, visualised using the recurrence plot (top), and the raster plot of population activity (bottom) sorted into detected ensembles. a Spontaneous perturbation away from the periodic attractor. Blue line indicates the final 5s window with more than 50% recurrence, long before the end of the recording. b Spontaneous perturbation returning to the same manifold of the periodic attractor. The orange lines indicates the point of maximum divergence from the manifold. The recurrence plot shows that same trajectory was recovered thereafter (white patches in upper right corner), visible as the return of the phasically firing ensembles (raster plot). c Spontaneous perturbation moving to a new periodic attractor. The recurrent trajectory up to 30s poststimulus did not recur thereafter (black region in upper right of the recurrence plot), but was replaced with a new recurrent trajectory. The raster plot shows this was an apparent transition in the oscillation period, and/or firing pattern of the ensembles.  Figure S6: Ruling out P10 motorneurons in the recorded population. We aimed to decode the spiking activity of the motorneuron axons in the P10 nerve directly from the low-dimensional trajectory of population activity in the pedal ganglion (Fig. 5). If any P10 motorneurons had been captured in our population recording, these could have made a disproportionate contribution to the decoding. We thus first sought to detect them, and remove them if necessary. We reasoned that any putative P10 motorneuron should have close to a 1:1 locking between its spikes in the population activity recordings and corresponding spikes in the P10 nerve recording; in other words P (spike10|spike) ≈ 1. Moreover, we expect a putative P10 motorneuron to show ∼1:1 lockng consistently across all three programs from the same preparation. For each neuron in each recorded program with simultaneous P10 activity, we computed P (spike10|spike) separately for a range of delays (0 to 25ms) between neuron and P10 spikes, allowing for ±2ms jitter at each delay. We plot the maximum probability for each neuron in programs 1 and 2 (a), and programs 1 and 3 (b); symbol size is proportional to the mean number of spikes for that neuron in the two programs (as few spikes would indicate an unlikely motorneuron). Colours indicate neurons from the three different preparations. We see that no neuron in our dataset comes close to P (spike10|spike) ≈ 1 at any delay, and there are no neurons with consistently high locking between programs. We concluded that there are no P10 motorneurons in our population recordings. MAE: median absolute error of fit between P10 activity and best-fit model. R 2 : correlation between P10 activity and best-fit model. b Relationship between mean model fit to and forecast of the scale of activity during the data-window (fit: 40s; forecast: next 10s). Success in fitting the scale of activity was predictive of forecasting.
c Relationship between mean model fit to and forecast of the change of activity (measured by R 2 ). Success in fitting the change in activity was not always strongly predictive of forecasting the change.  Figure S8: Increasing dimensionality of state-space did not improve the P10 decoding model. Here we show the effects of increasing the embedding dimensions for the population activity before using the decoding model for P10 activity. We used enough dimensions to capture at least 90% of the variance; we used 80% in the main text. a Comparison of the number of embedding dimensions needed to explain 80% and 90% variance for the nine recordings with simultaneous P10 activity. A 10% increase in variance needed at least double the number of dimensions. b Example fit and forecast by the statistical decoding model for P10 firing rate when using more embedding dimensions. Grey bar indicates stimulation time. c Summary of model forecasts for all 9 population responses with P10 activity, when using more embedding dimensions. Dots and lines show means ± 2 SEM over all forecast windows (N = 173). MAE: median absolute error (error in scale) of forecast; R 2 : variance explained of forecast. d Comparison of the mean MAE in forecasting P10 activity when using embedding dimensions explaining 80% and 90% variance. Increasing the number of dimensions did not systematically improve the forecasting of the scale of the P10 activity. e Comparison of the mean R 2 in forecasting P10 activity when using embedding dimensions explaining 80% and 90% variance. Increasing the number of dimensions did not systematically improve the forecasting of the variance of the P10 activity.  Figure S9: Participation measures rate and synchrony. a Example correlation between the evoked firing rates of all neurons in one program and the participation of those neurons in all three programs evoked from that preparation. b Histogram of firing rate and participation correlations across all 30 programs. Firing rate explains on average (R 2 , red line) 50% of the variation in participation within a preparation; the remaining variation in participation can be accounted for by differences in synchrony between individual neurons and the low-dimensional projection.