A freely-moving monkey treadmill model

Objective. Motor neuroscience and brain–machine interface (BMI) design is based on examining how the brain controls voluntary movement, typically by recording neural activity and behavior from animal models. Recording technologies used with these animal models have traditionally limited the range of behaviors that can be studied, and thus the generality of science and engineering research. We aim to design a freely-moving animal model using neural and behavioral recording technologies that do not constrain movement. Approach. We have established a freely-moving rhesus monkey model employing technology that transmits neural activity from an intracortical array using a head-mounted device and records behavior through computer vision using markerless motion capture. We demonstrate the flexibility and utility of this new monkey model, including the first recordings from motor cortex while rhesus monkeys walk quadrupedally on a treadmill. Main results. Using this monkey model, we show that multi-unit threshold-crossing neural activity encodes the phase of walking and that the average firing rate of the threshold crossings covaries with the speed of individual steps. On a population level, we find that neural state-space trajectories of walking at different speeds have similar rotational dynamics in some dimensions that evolve at the step rate of walking, yet robustly separate by speed in other state-space dimensions. Significance. Freely-moving animal models may allow neuroscientists to examine a wider range of behaviors and can provide a flexible experimental paradigm for examining the neural mechanisms that underlie movement generation across behaviors and environments. For BMIs, freely-moving animal models have the potential to aid prosthetic design by examining how neural encoding changes with posture, environment and other real-world context changes. Understanding this new realm of behavior in more naturalistic settings is essential for overall progress of basic motor neuroscience and for the successful translation of BMIs to people with paralysis.


Introduction
Animal models are central to investigating how the brain senses the ever-changing world and responds with a large repertoire of possible behaviors. Animal models that focus on naturally occurring sensations and behavior are increasingly desirable, as the limitations of studying only simplified sensory-motor functions and their neural underpinnings is increasingly recognized as stifling progress. For example, in sensory neuroscience, cats watching movies which capture natural scene statistics expanded our knowledge of how the visual system processes naturalistic images and time series [1]. In systems neuroscience, rhesus monkeys are often chosen as the model species because of their ability to perform a wide variety of sensory, cognitive and motor tasks including, for example, dexterous arm movements [2]. However, despite the rhesus monkeyʼs capacity for a large repertoire of naturalistic behaviors (e.g., spontaneous, non-instructed peeling and eating of a banana while swinging in a hammock), most motor studies have historically focused on a relatively small number of instructed arm movements.
Motor studies of this sort have focused on so few behaviors because of the technologies needed for performing the necessary measurements. Monkey models for motor neuroscience require two simultaneous measurements, neural and behavioral. Traditional experimental frameworks record neural activity through tethered (wired) connections (e.g., [3,4]) requiring a reduced range of motion, such as a fixed head posture and restricted arm movements, to protect both the monkey and the sensitive recording electronics, as illustrated in figure 1(a). We refer to these traditional monkey models as 'head-fixed models' because the head is often braced in place to prevent the monkey from pulling on the tethered connection or vibrating the head or chair which could introduce micro-phonic electrical noise into the neural recordings. Methods for gathering behavioral measurements can also be cumbersome for less restricted or completely free movement. A highly productive head-fixed model, and therefore posture-fixed model as well, is the reaching monkey model where monkeys perform 2D or 3D arm movements and reaching behavior is captured by tracking the endpoint of the arm by target location (e.g., [5]), with a manipulandum (e.g., [6]), or with an infra-red camera tracking system (e.g., [7,8]), by tracking markers on the whole arm [9], or by direct measurement of the arm position while the arm is in a two degree-of-freedom exoskeleton using joint-angle encoders [10,11]. These methods limit movements to the region of space where the arm can move but without changing posture and to the ranges of the behavioral sensors. These methods are sufficient for studying 2D and 3D movements such as point-to-point reaching within the aforementioned restricted workspace.
However, fixed posture 2D and 3D arm movements are a limited subset of all natural behaviors, and it is unclear how well scientific and BMI results from limited workspace experiments generalize to broader, more naturalistic behaviors in unconstrained workspaces. For example, though fairly high levels of BMI cursor control and prosthetic arm control is now possible (e.g., [12][13][14][15]), it is unclear how well these results based on head and shoulder immobilized rhesus monkey models and initial clinical trials in people with paralysis (tetraplegia) will generalize to cases where more movement remains intact (e.g., single-arm amputation, so other arm and legs can still move) and thus much larger workspaces are relevant and potential interaction among postural, movement, cognitive, and environmental context signals are possible. Moreover, recent studies demonstrate that the relationship between neural activity and movement is complex [16,17]. Thus, when motor neuroscience experiments focus on a small number of reaches, they may limit our ability to decipher this complexity in neural responses.
One potential solution is to develop additional head-fixed monkey models that examine other behaviors, such as the bipedal walking model [18]. However, examining more behaviors in isolation might not untangle the complexity of responses without knowing how neural ensembles respond across behaviors. To examine how neural ensembles in the motor cortex relate to a wide variety of voluntary movement, we propose that a more flexible monkey model is needed. We refer to these flexible models as 'freely-moving monkey models'. Compared to head-fixed monkey models, freely-moving monkey models are more difficult to design and make operational, but they may well involve less animal training. These monkey models require recording neural activity and behavior using technologies that do not encumber movements, as illustrated in figure 1(b). Therefore, traditional neural recording systems with a tethered connection must be replaced with head-mounted recording systems [19]. A number of systems have been designed for such telemetered neural recordings (i.e., wireless transmission of neural data from the electrodes/head to a nearby receiver in the room) as interest in gathering these types of measurements has increased [20][21][22][23][24][25][26][27][28][29][30][31]. Behavior must also be recorded in a larger environment than in traditional head-fixed models, for example in a home cage (in-cage model) or in the wild (inthe-wild model). One possibility for gathering behavioral measurements is by using computer vision techniques (e.g., [32,33]). Though the technological requirements for freelymoving monkey models are more complex than for headfixed monkey models, studying natural behavior may require less training time than reaching monkey models in which months can be spent teaching monkeys to learn tasks. A comparison of the technologies used between different monkey models is presented in table 1.
In this report, we present a novel freely-moving monkey model, termed the treadmill model. The treadmill monkey model uses the two key freely-moving monkey model technologies, wireless neural recordings and computer vision based behavioral measurements, to study the neural correlates of a monkey walking quadrupedally on a treadmill. This experimental framework is important in its own right potentially as a useful animal model for gait neurorehabilitation (e.g. [34]), as well as an important step towards designing even more complex in-cage and in-the-wild monkey models, which require even more complex and robust neural recording and behavioral tracking systems (see Discussion). We have previously reported, preliminarily, that neural activity encodes walking and reaching [35,36]. Here, we substantially extend these previous preliminary reports by demonstrating results in two rhesus monkeys (which is critical for knowing if the scientific results generalize across monkeys), analyze how neural activity covaries with changes in walking speeds, and analyze population level activity so as to report the first comprehensive study of the treadmill monkey model and associated scientific results. The kinematics of quadrupedal walking have been characterized in monkeys (e.g. [37]), and the neural correlates of quadrupedal walking have been studied in cat animal models [38][39][40], but, to our knowledge, this is the first time that the neural correlates of quadrupedal walking has been studied in monkeys.

Methods
Walking behavior was recorded by eight video cameras while neural activity was telemetered from chronically implanted electrode arrays. The recording environment consisted of a commercially-available exercise treadmill surrounded on four sides by transparent plastic walls held together in a custombuild extruded aluminum enclosure. The back wall was removed when the monkey was walking into the treadmill environment and clipped back in place during the experiments. The front wall had a hole to supply water and treats to the monkey. The side walls and front wall were flush against the treadmill and there was space in the back of the environment for the animal to step off the treadmill. The top of the environment was open. A diagram of the recording environment is shown in figure 2(a), pictures of the environment are shown in supplementary figure 1 (available from stacks.iop. org/jne/11/046020/mmedia), and the recordings are described in detail below.

Behavioral task
All protocols were approved by the Stanford University Institutional Animal Care and Use Committee. Two rhesus macaques (monkeys I and N) were trained, over a period of several weeks, to walk naturally and quadrupedally on a treadmill moving at various normal walking speeds. After initial training sessions in which each animal was introduced to the treadmill environment, we began recording sessions in which the animal walked while video and neural data were captured. Each recording session (one per day per monkey) lasted about 1 min. A session began by bringing a monkey to the treadmill environment using a pole attached to the collar, a standard technique for safely assisting monkeys while moving in a laboratory setting. For precautionary safety reasons, the pole remained loosely attached during walking without disturbing behavior as no force was applied. Each session was divided into 6 to 10 walking epochs in which the treadmill moved at a fixed speed, and the fixed speeds were repeated and interleaved. During an epoch, the monkey walked unencumbered on the treadmill and was free to exit off the rear end of the treadmill if desired. To encourage continuous walking, the animal was rewarded continuously with water or fruit at the front of the treadmill. Times within an epoch in which the animal exited the treadmill were not analyzed. Between each walking epoch, the animal reached from the stationary treadmill to food rewards at the front of the treadmill. We do not believe, and we had no observations suggesting otherwise, that the water and food rewards had a significant effect on walking behavior. A total of 40 (46) recording sessions were captured from monkey I (N) over the course of 12.5 (8.5) months. During these sessions, the behavior of the monkeys was carefully monitored to ensure animal safety, and there were no adverse events. These recording sessions constitute a large dataset, on par with data sets collected from monkeys with chronically-implanted arrays performing reaching tasks in 'rig'.
For analysis, walking behavior was divided into individual steps. Each step began and ended when the right arm (contralateral to implant) returned to the treadmill after swinging forward through the air. To increase the consistency of behavior in analyses and to ensure a high count of steps for individual speeds, we focused on a range of walking speeds in which the animals took infrequent breaks. This range of walking speeds was divided into three speeds categorized as slow, medium, and fast walking corresponding to 2.0, 3.0, and 4.0 kph for monkey I (2.0, 2.5, and 3.0 kph for monkey N). Monkey I is larger than monkey N, so it is natural for them to have a different set of comfortable walking speeds, and the key point in this investigation is to have a slow, medium, and fast speed for each monkey. These speeds need not be the same across monkeys as we are not making statements about neural correlates of absolute speed.

Video recordings
Walking behavior was recorded on video for offline analysis from eight camera views. Point Grey Grasshopper (GRAS-20S4M/C) cameras captured video at 24.7 frames per second at a resolution of 1624 × 1224 pixels. Video acquisition and export were performed using a 4DViews 2DX Multi-Camera system. The 4DViews camera system ensured that images were captured synchronously by all eight cameras via a synchronization line. The synchronization line was also used to synchronize the neural recordings with the video recordings as shown in black in figure 2(a) and described in the following section. Additionally, the intrinsic and extrinsic parameters of the cameras were computed using a custom multi-camera calibration method. Behavioral measurements were gathered from video frames by using computer vision techniques and by handtagging frames of interest. Computer vision techniques were used to automatically find the 3D location of the right shoulder, elbow, and wrist as the monkey walked. Unlike many neuroscience studies where markers attached to or painted on the animal aid joint tracking (e.g., [18] and [9]), we did not include markers because they can distract the animal. Instead, we developed a computer vision algorithm that is based on statistical learning, epipolar geometry (the geometry of two camera views [41]), anatomical constraints, and temporal coherence. Specifically, we trained viewdependent detectors of these arm joints in each of the (2D) images (similar to how object detectors were trained in [42]). Within each view, each arm joint (e.g., the elbow) was associated with multiple such detectors-each corresponding to the arm joint in a different 3D pose. The per-pixel detection response for an arm joint was defined to be the strongest perpixel response among all the detectors associated with the arm joint. The detection response for an arm joint for all pixels defined a likelihood map for the 2D location of the arm joint in a given camera view.
To estimate joint positions in a new unseen video sequence, we began by computing the responses of the arm joint detectors. Note that in a given video frame, an arm joint may be occluded in one view. However, our use of several cameras ensured it was always visible in at least some of the views. Next, in order to suggest candidates for the 3D locations of the arm joints, we looked for 3D locations whose projections onto the 2D images correspond to strong responses of the 2D detectors across the different views. Importantly, by using epipolar geometry, we were able to reduce the search space: we started from the 2D candidates in the image taken from the most informative camera (sideview), and then, by restricting the search along epipolar lines (given calibrated cameras), we only needed to search along 1D lines in other views. This step results in a set of 3D candidates (weighted by the responses of the detectors in all views) for each one of the arm joints (shoulder, elbow, and wrist). Finally, using temporal coherence, the structure of the kinematic tree (e.g., the shoulder location restricts the possible locations of the elbow), and anatomical constraints (obtained via biometric measurements of the animal), we optimized, using dynamic programming, for the full trajectory of the 3D locations of the arm joint across the entire video sequence. Performing inference in a batch process (over the full trajectory as opposed to frame by frame) enabled us to reduce noise and to compensate for joint-detector errors.
In addition to automated behavioral tracking, a custom MATLAB (Mathworks, Natick, MA) interface was written to hand-tag frames of interest in the video. For example, we labeled the frames that marked the transitions between each step.

Neural recordings
In each monkey, a 96-channel multielectrode array (Blackrock Microsystems, Salt Lake City, UT) was implanted in the dorsal aspect of the premotor cortex (PMd) as determined by visual anatomical landmarks. Neural signals were telemetered using one of two recording platforms for each recording session: the HermesD system or the HermesE system [22,23]. HermesD, used in 30 (40) recording sessions for monkey I (N), records 32 channels of broadband data at 12 bits per sample at 30 kSamples per second, and HermesE, used in 10 (6) recording sessions for monkey I (N), records 96 channels of broadband data at 10 bits per sample at 31.25 kSamples per second. As shown in figure 2, the respective Hermes receiver demodulated the transmitted neural data and output the result to a field programmable gate array (Spartan-3A FPGA; Xilinx, San Jose, CA) located on a ZestET1 (Orange Tree Technologies Ltd, Didcot, UK). The FPGA was programmed to check for wireless transmission errors, to remove non-informative stuffing bits, and to package the useful data onto a User Datagram Protocol (UDP) Ethernet stream to be stored for offline analysis. The FPGA also received the video camera synchronization line as an input, included the image timing information in the output UDP datastream, and output a specified pattern to synchronization LEDs that could be seen in multiple camera views. By checking the patterns on the LEDs in multiple camera views, we certified that the synchronization as measured by the hardware matched the images saved by the video recordings.
Broadband neural data were then processed to identify and extract action potential (spike) timing information for analyses. Raw broadband waveforms are plotted for 100 ms in figures 3(a)-(b) for two channels in monkey N recorded with HermesD. LFP, not analyzed in this study, was filtered out using a highpass filter with cutoff frequency of 250 Hz. Spike times were designated anytime the filtered waveform crossed a threshold of −4.0 times the RMS voltage value of the channel, similar to prior BMI research in our group and others (e.g., [43]). Filtered waveforms for 3.0 s are shown in figures 3(c)-(d) for the same two channels with the threshold for spikes marked in black. A raster plot of 10 s of walking is plotted in figure 3(e).

Dataset selection
The environment used in the treadmill monkey model (i.e., a part of a traditional rhesus monkey 'housing room' with many other metal cages also in the room) contains more sources of noise than the controlled rig environment of traditional monkey models (e.g., RF, light, and acoustically shielded rooms). In this treadmill environment, we contended with artifacts in the neural data that arose from wireless transmission errors, grounding issues, thermal sensitivities of the receiver, and interfacing issues between the Hermes systems and the neural implant. Many of these issues were mitigated over the course of our year of recordings. We chose to focus on the datasets with the fewest artifacts. Because of the time consuming nature of hand-tagging video data, we selected a limited number of representative datasets to analyze. For monkey I, we focused on 9 walking epochs from a single day (I20121024); 3 epochs were from each speed slow, medium, and fast, and neural recordings were made using HermesE. For monkey N, we focused on 6 walking epochs from two consecutive days (N20120821 and N20120822); each day had one slow, one medium, and one fast epoch, and neural recordings were made using HermesD.
The results presented in this paper are taken from representative datasets from both monkeys. While we focus on a subset of the recordings for the practical reasons mentioned above, the results presented here are similar to recording sessions from many days. In general, there is no guarantee that threshold crossings across many days will be stationary, and therefore we chose to analyze a single day for monkey I and two consecutive days for monkey N. A full analysis of the stability of waveforms is beyond the scope of this paper and has been performed previously with the same type of multi-electrode arrays [43]. Figure 4(a)-(c) plots the right shoulder, elbow, and wrist in the vertical plane for 102 frames as monkey I walks at a slow, medium, and fast pace. Because there was no 'ground-truth' measurement against which to compare to the computer vision results, we compared the computer vision 3D positions to hand-labeled 3D locations for 100 points. We found that the average Euclidean distances between the computer vision 3D data points and the hand-labeled points were 26.1 mm, 20.6 mm, and 17.1 mm for the wrist, elbow, and shoulder, respectively for monkey I. For reference, monkey Iʼs arm (shoulder to wrist) is approximately 42 cm. Errors were anisotropic in 3D coordinates. For each joint, the error in the horizontal direction perpendicular to the direction of movement was greater than the errors in the vertical direction and in the walking direction. For instance, the root-mean-square error between the hand-labeled points and computer vision points for the wrist were 9.6 mm in the direction of walking, 10.5 mm in the vertical direction, and 21.9 mm in the horizontal direction perpendicular to the direction of movement. The poorer performance of this direction is likely due to occlusion by the body. We measured the root-mean-square error of the elbow joint angle in the vertical plane to be 4.7 degrees. For monkey N, the automated tracking was less successful due to the animalʼs smaller limbs, more motion blur, and irregular walking patterns, and some points were hand-corrected to obtain more accurate tracking. For this monkey, manual corrections were required for about one frame out of every 50 frames. This represents a substantial improvement over manually determining all points. No manual corrections were required for monkey I.

Behavioral tracking using markerless computer vision
The frames that represent the beginning and end of each step are highlighted in figures 4(a)-(c). These plots represent 4.1 seconds of walking, during which monkey I takes about three steps when walking at a slow speed, about four steps when walking at a medium speed, and about five steps when walking a fast speed. To summarize the behavior over the analyzed datasets, the cumulative distribution function (CDF) of step rates, the inverse of each step duration, is plotted for monkey I in figure 4(d) and monkey N in figure 4(e). Monkey I has a more consistent step rate as shown by a steeper slope in the CDF than monkey N.

Neural spiking in PMd encodes walking phase
A central aim of motor neuroscience is understanding which behaviors are encoded by neural activity. In head-fixed monkey models, the standard approach has involved repeating the identical behavior for many trials while recording neural activity. Consistent features in neural recordings across trials are considered encodings of the behavior, and inconsistent features are considered noise. In freely-moving monkey models, there is no trial structure to exploit. However, we can classify behavioral epochs into categories (e.g. reaching or walking) and ask if the neural activity during these epochs encodes the corresponding behaviors. Neural activity encodes behavior if there are statistical differences in the neural activity between two behavioral conditions.
In other words, we narrowly define 'encoding' to indicate that the mutual information between the neural activity and the categorical behaviors is greater than zero. The manner in which information is encoded remains ambiguous. For example, motor cortex is heavily connected to both somatosensory cortex and output motor areas (e.g., the spinal cord). When we find that motor cortex encodes behavior with this method, we have not determined whether we are measuring the encoding of the proprioceptive signals or the output motor command (if a distinction can be made). Instead, we have determined that motor cortex contains information that can differentiate between the behaviors, and thus, encodes the behavior though we may not know how.
In the treadmill monkey model, we asked whether neural spiking in PMd encodes the phase of quadrupedal walking. We divided continuous walking of each monkeyʼs walking epoch (see Methods) into individual steps, which begin and end with the right arm (contralateral to implant) coming into contact with the treadmill after it is swung forward. A full step cycle is further divided into two epochs: the stance phase when the right arm is in contact with the treadmill followed by the swing phase when the right arm swings forward through the air. In figure 5, raster plots for threshold crossings are shown for 96 (32) channels for monkey I (N) while walking for 5 seconds at a medium (fast) speed. The channels are arranged by the timing of the peak firing rate within the phase of walking. The phase of walking is also marked below each raster plot with the stance phase in black and the swing phase in white. In the raster plots it is clear that a periodicity exists in the threshold crossings that closely matches the periodicity in the phases of walking and that different channels have different phase relationships with the phases of walking.
We further demonstrate the relationship between the periodicity of spiking and walking by aligning neural spiking on an individual channel to individual steps. We have chosen multi-electrode array technology to emphasize simultaneous recordings across many channels in order to investigate population dynamics. One drawback of this technological choice is that signal quality degrades over time, and we cannot always isolate single neurons. However, we were able to isolate single neurons on a few channels, and they appear to have a simple relationship with the phase of walking (see supplementary figure 2). This relationship tends to be stable across days (see supplementary figure 3). In figure 6, raster plots and firing rates are shown for a hand-sorted single-unit on a single electrode (channel) for each step as monkey I and N walk at slow, medium, and fast speeds. Individual steps are aligned at t = 0, which is when the right arm comes in contact with the treadmill after the swing phase of walking. These single units appear to have the same relationship with the phase of walking independent of speed.
Finally, we investigated to what extent PMd encodes the phase of walking. Because many channels did not have single units that we could isolate, we chose to compare threshold crossings so that we can study higher channel counts. Blurring information of multiple units on a single electrode due to threshold crossing can only hurt our measurement of the extent to which motor cortex is involved in walking. As a simple quantitative measure to determine if a channel is tuned, we compare the square root (a standard transformation for making Poisson processes more normal [44]) of the average firing rate of threshold crossings during the stance phase and the swing phase (see supplementary figure 4). We exclude neural channels with an average firing rate below 5 spikes per second in both phases of walking, with artifacts Frames where the right wrist comes in contact with the treadmill are highlighted; they represent the end of the swing phase and the beginning of the stance phase. The cumulative distribution functions (CDF) of step rate for slow (red), medium (blue) and fast (green) walking (measured from point of contact to the next point of contact) are plotted in panels (d) for monkey I (I20121024) and (e) for monkey N (N20120821 and N20120822). Monkey I walks at a more consistent rate for a given speed than monkey N. The phase of the gait are plotted below the raster plots with the stance phase in black and the swing phase is grey. To sort the channels, first we calculated the firing rate of each channel for each step (consisting of a stance phase followed by a swing phase). Because each step has a different duration, we normalized the step timing to its phase from 0 to π 2 . Then we averaged the firing rate for each channel across time-normalized steps. Finally, we found the phase where the maximum average firing rate occurred and sorted channels by this value. (b) Threshold crossings as monkey N walks at 3.0 kph (fast speed) for 5 s. due to an amplifier saturating, or that were disconnected due to interfacing issues. These channel selection criteria reduced the number of channels from 96 (32) to 76 (20) for monkey I (N). The fewer active channels in monkey N data is due to poorer signal quality of the electrode array, which degraded more quickly over time than the signal quality of the array implanted in monkey I. We found that 70/76 (11/20) channels are significantly tuned ( < p 0.01, a two-sample T-test) for monkey I (N). Thus, most channels in both monkeys encode the two phases of the gait cycle during quadrupedal walking.

Average firing rate covaries with wrist velocity during step
In addition to showing that neural activity encodes certain behaviors, motor neuroscience seeks to describe how neural activity covaries with behavioral kinematics. In both headfixed and freely moving monkey models, differences between similar movements, such as speed and trajectory, are regressed against neural activity to investigate how neural activity changes with changes in behavior. Using head-fixed reaching monkey models, it has been shown that differences in reaching speeds for similar reach trajectories are associated with differences in neural firing rates [45,46], and we hypothesize that differences in walking speeds would likewise be associated with differences in firing rates.
One might ask whether neural activity covaries with the speed or position of the wrist across all time points. This is equivalent to determining how well neural activity can predict the position of wrist over time, a measurement made by Fitzsimmons and colleagues during bipedal walking [18]. We, too, find that neural activity covaries well with wrist speed (data not shown). However, this is equivalent to saying that the neural activity is tuned to walking as we have shown in the previous section. If neural activity is tuned to the phase of walking and if the wrist position and speed co-varies systematically over each movement, then neural activity will covary with speed as well.
We tested to see if the average firing rate during individual steps is correlated significantly with the peak forward wrist speed of the step. This requires a single-trial analysis examining how well neural activity can account for trial-bytrial changes in wrist speed. For each step, we calculated the peak forward speed of the wrist from the computer vision measurements of walking behavior. Firing rates are calculated by convolving each channelʼs spike train with a Gaussian window with a standard deviation of 50 ms, and the average firing rate over the step duration was measured. Channels excluded in the previous section were again excluded. A total of 323 (173) steps are included for monkey I (N). A linear regression was performed to determine whether there is a significant relationship between the average firing rate and the peak forward speed of the wrist during walking. Example channels from each monkey are plotted in figures 7(a), (b). Of the 76 (20) channels analyzed, 52 (7) were found to have significant slopes ( < p 0.01) in the regression for monkey I (N); channels were sorted by R 2 value, and the R 2 values are plotted in figures 7(c)-(d) with slope values inset. Despite recording from PMd in both monkeys, in monkey I, most channels exhibit an increase in average firing rate for steps that have faster swings, and in monkey N, most channels exhibit a decrease in average firing rate for steps that have a faster swing. This could be due to, among other reasons, slightly different placements of the electrode array with respect to anatomical landmarks (i.e., intra-operative variation) or due to different function localization in each monkey with respect to anatomical landmarks. Of the 52 (7) channels that significantly covary with wrist speed, 48/52 (4/7) are also found to be encoded to the different phases in walking in the previous section for monkey I (N). Thus, channels that were found to significantly covary with peak wrist speed were likely to encode the phases of walking.

Neural trajectories are rotational, evolve at step rate, and separate by speed
The previous two sections examine the relationship between multi-unit spiking on individual channels and quadrupedal walking, but it is equally (if not potentially more) important to understand how the entire population relates to movement (e.g., [47]). Recent work in motor neuroscience has used dimensionality reduction techniques to view how neural activity of the population evolves. The top two or three dimensions calculated from dimensionality reduction can be plotted together, and the evolution of the neural population firing rates over time can be viewed as neural trajectories in this space (e.g., [44,48]). We performed factor analysis on multi-unit firing rates as monkey I walks at different speeds. We could not perform this analysis on data from monkey N because there are not enough active channels to build neural trajectories. A single step is plotted in figure 8(a) with the stance phase of walking in bold followed by the swing phase of walking. Neural trajectories for 10 steps at slow, medium and fast speeds are plotted in red, blue, and green in figure 8(b), respectively. The neural trajectories for steps at different speeds share similar rotational trajectories though there are small systematic differences in trajectories for different speeds. These small differences are likely caused by the number of channels that covary with step velocity as shown in the previous section. A major difference between the neural trajectories for different walking speeds is their angular frequency. The average angular frequency of the neural trajectories is measured for the 10 steps for each walking speed and plotted in figure 8(c). The average step rate (measured in steps per second) for the 10 steps is plotted in figure 8(d). Converting units from radians per second to cycles per second, it is apparent that the angular velocity of the neural trajectories occur at the same rate as the step speed.
Though prominent features of the neural trajectories seem to be consistent across walking speed, systematic differences may be embedded in the population neural trajectories. Recently it has been suggested that a separation of rotational dynamics may exist as a mechanism for neural populations to generate oscillatory output signals from inputspecified target frequencies [49]. We investigated whether there was a systematic offset in any direction by searching for a direction, z, in state-space that best separates the population between walking speeds.
There are a number of methods to separate two classes, including linear discriminant analysis, logistic regression, and support vector machines. We found these methods work well to separate the neural trajectories for slow and fast walking (data not shown). However, these methods do not incorporate the neural trajectories for walking at medium speed. Therefore, to find the separation direction z, we used least-squares to predict the walking speed from each point in the trajectory for a single walking epoch (continuous walking for approximately 40 steps) at each speed. Indeed, the neural trajectories are well separated when projected onto z, as shown in figure 8(e). This method outperformed two other previously described class separation methods on cross validated data. The overall variance in this separation direction is quite low; for example, 0.18% of the total neural variance for slow walking is in the direction of z compared to 30.87% and 20.94% of the total neural variance is in the direction of the first factor x 1 and the second factor, respectively. The distribution of the training dataset separated well on this axis finding a separation in the noise. We then projected testing data composed of the six remaining epochs not used to learn the direction z onto the learned direction z to verify that the separation was real, figure 8((f), middle and bottom). Finally, we illustrate the large separation of slow and fast speeds by showing a receiver operating characteristic plot, which shows how well we classify neural trajectories into the correct speeds as a function of decision threshold (see figure 8(g)).
The dotted curve represents the separation of slow and fast walking for the training neural trajectories, and the two black curves for the test neural trajectories. For the testing datasets, we can correctly separate more than 85% of the slow and fast walking trajectories by a linear classifier. We find that, though neural trajectories have largely similar features for different speeds, they have a small, but robust separation along the direction z.

PMd modulated during walking
Many textbooks describing the control of locomotion in mammals focus on the role of central pattern generators in the spinal cord (e.g., [50]). The role of cortex is often described as being used for skilled walking, such as walking over obstacles or on ladders (e.g., [51]). Studies using decerebrate cat models have shown that unskilled walking need not have cortical inputs to produce walking movements due to central pattern generators, but to the best of our knowledge no studies have effectively demonstrated the existence of central pattern generators in higher order primates with well developed corticospinal tracts [52]. Thus, we feel that treadmill locomotion is not only a reasonably accessible behavior to study because it is relatively stereotyped, but also a behavior that can lead to insights into motor function. If corticospinal control is more important in these higher order primates, it is unsurprising that in our data we found that most of the recorded channels in PMd encoded the phase of walking. Here again, we point out the ambiguity of the encoding. The neural signals in premotor cortex may arise from sensory regions, or they may be controlling the locomotion, or both. Although this ambiguity is unresolved, it is clear that the information about the phase of walking is present in threshold crossings in PMd.
On the other hand, if walking were truly controlled solely by the spinal cord, then it is unclear what purpose this modulation within motor cortex serves. One possible explanation comes from studies that have shown that successful skilled walking, such as walking on ladders and responding to changing terrain, requires cortical input. This means that cortex is required when changes to gait patterns are required. Therefore, modulation found within motor cortex may represent tracking of the phase of walking such that appropriate corrective movements can be instructed when intervention is required.

Implications for studies of voluntary movement
Scientific endeavors to understand how the brain controls voluntary movement have relied on using animal models that narrowly focus on a handful of behaviors. However, the neural responses to this limited set of behaviors is often complex and does not always follow simple, interpretable relationships. The hope for freely-moving models is to be able to find consistent relationships across many more sets of voluntary behaviors.
Studying the control of voluntary movement in this manner requires a new strategy for analyzing neural correlates. Instead of focusing on a repeated behavior for many trials, freely-moving animal models require recording many behaviors over an extended period of time, then using modern algorithms for large datasets to find consistent patterns within the data. To be analyzed efficiently, new statistical and machine learning techniques are required to study the neural correlates of behavior. One approach is to segment behavior into categorical epochs (e.g., reaches upwards) where behavior is similar, and then to analyze those behavioral epochs for common neural motifs (see supplementary figures 5 and 6). Conversely, we can use machine learning techniques to look for patterns in neural data and analyze the behavioral correlates of those patterns. Here, we present the types of data inquiries that can be used in freely-moving animal models for the treadmill monkey model. The treadmill monkey model has a number of advantages and disadvantages over in-cage and in-the-wild monkey models. The benefit of the treadmill monkey model is that similar behaviors (steps) are easy to collect in a short amount of time because the animal walks continuously. This simplifies the data storage problem because only a limited time is needed to collect sufficient behavior. In addition, computer vision for a task like walking is much easier than for tasks in the cage and in the wild where occlusion from cameras can make the computer vision task more difficult. For these reasons, we chose to limit the scope of this initial study to walking. However, to realize the full potential of freelymoving animal behavior, these systems need to be designed for and rigorously evaluated in multiple behavioral contexts. Though the complexity of developing in-cage and in-the-wild monkey models is higher, they have the benefit of lower complexity in gathering raw data. Data acquisition requires little intervention once the neural recording and behavioral monitoring systems are set up. The complexity is shifted from difficult recordings to difficult data analysis. However, incage and in-the-wild monkey will produce richer data leading to new insights in the control of voluntary movement.

Implications for BMIs
Given the inherent complexity of neural correlates of movements, it is difficult to determine how the neural population will respond to novel behaviors. For computer cursor BMIs, the space of novel behaviors is reduced because the cursor is confined to a 2D plane, and thus reaching monkey models have provided perhaps a sufficient amount of data to train prosthetic decoders-assuming that the person with paralysis is tetraplegic, and thus is approximately as immobile as the head-fixed (and thus shoulder-constrained) monkey model. However, to design decode algorithms for whole-arm BMIs capable of working with people who are ambulatory and have a single-arm amputation, more knowledge about neural correlates of complex reaches is needed. People who use wholearm prosthetics will likely need to change postures or move to new environments; the effect of which is unknown on neural correlates. Freely-moving monkey models afford us the opportunity to study these factors.

Conclusion
In this work, we have presented the design of a freely-moving monkey model. By using telemetered neural recording technologies and markerless computer vision techniques, we were able to study the neural correlates of quadrupedal walking for the first time. The techniques presented here build a foundation for freely-moving monkey models that future work can build on to study the neural correlates of many behaviors providing a comprehensive look at how the brain controls voluntary movement.