Abstract

Do motor commands encode movement independently or can they be represented in a reduced set of signals (i.e. synergies)? Motor encoding poses a computational and practical challenge because many muscles typically drive movement, and simultaneous electrophysiology recordings of all motor commands are typically not available. Moreover, during a single locomotor period (a stride or wingstroke) the variation in movement may have high dimensionality, even if only a few discrete signals activate the muscles. Here, we apply the method of partial least squares (PLS) to extract the encoded features of movement based on the cross-covariance of motor signals and movement. PLS simultaneously decomposes both datasets and identifies only the variation in movement that relates to the specific muscles of interest. We use this approach to explore how the main downstroke flight muscles of an insect, the hawkmoth Manduca sexta, encode torque during yaw turns. We simultaneously record muscle activity and turning torque in tethered flying moths experiencing wide-field visual stimuli. We ask whether this pair of muscles acts as a muscle synergy (a single linear combination of activity) consistent with their hypothesized function of producing a left-right power differential. Alternatively, each muscle might individually encode variation in movement. We show that PLS feature analysis produces an efficient reduction of dimensionality in torque variation within a wingstroke. At first, the two muscles appear to behave as a synergy when we consider only their wingstroke-averaged torque. However, when we consider the PLS features, the muscles reveal independent encoding of torque. Using these features we can predictably reconstruct the variation in torque corresponding to changes in muscle activation. PLS-based feature analysis provides a general two-sided dimensionality reduction that reveals encoding in high dimensional sensory or motor transformations.

Author Summary

In studying animal movement, one cannot always record all the motor commands an animal uses or know all the ways in which movement varies in response. A combined approach is necessary to find the relevant patterns: the changes in movement that correspond to changes in the recorded motor commands. Techniques exist to identify simple patterns in either the motor commands or the movements, but in this paper we develop an approach that identifies patterns in both simultaneously. We use this technique to understand how agile flying insects control aerial turns. The two main downstroke muscles of moths are thought to produce turns by creating a power difference between the left and right wings. The moth’s brain may only need to specify the difference in activation between the two muscles. We discover that moth’s brain actually has independent control over each muscle, and this separate control increases the moth’s ability to adjust turning within a single wingstroke. Our computational approach reveals sophisticated patterns of movement processing even in the small nervous systems of insects.

Introduction

To understand how motor spikes are transformed into action requires knowledge of how movement is encoded in these patterns of neuromuscular activation. Unfortunately, we cannot reliably predict the movement resulting from a particular motor signal simply from a muscle’s anatomy and static function (e.g. an “extensor”) [1—4]. The same pattern of activation to the same muscle, but in different dynamic contexts can even produce turning torques in opposite directions [2]. Moreover, both motor signals (the inputs to the motor transform) and movements (the outputs) are typically high dimensional and we may not be able to record all relevant motor signals electrophysiologically. Understanding how muscles work together to encode movement is therefore a computational challenge of both 1) dimensionality and 2) incomplete representation.

Dimensionality reduction techniques used in sensory neuroscience are typically one-sided, meaning that a high dimensional stimulus is reduced to describe variation in a single spiking neuron [5]. Similar use of one-sided dimensionality reduction in the motor transformation identifies patterns in high dimensional representations of movement encoded in the activity of individual muscles e.g. [6]. Conversely, the dimensionality of the neural signals can be high, as is the case for many central brain recordings, but in many experimental designs the representation of motor output is restricted to few dimensions, or even discrete states, e.g. [7]. Developing techniques that represent high dimensional movement encoded in high dimensional neural signals remains challenging [8], yet is ever more pressing as such datasets become the norm.

Muscle synergies, patterns of variation in activation across multiple muscles, are hypothesized to reduce the dimensionality of the motor commands and provide high level encoding of movement features [9,10]. However, not all variation in muscle activation might affect movement dynamics (and vice versa). If we consider all variation in a high dimensional description of movement as potentially relevant, we are likely to include variation that is not encoded by the muscles we are able to record from. This challenge exists for sensory encoding as well. Only certain features of a complex stimulus may be encoded by the neurons under consideration. Spike triggering is one way to extract only relevant variation. This method conditions the stimulus on the spiking of an individual neuron (or muscle) thereby limiting the reduced stimulus description to what is encoded in that spike. However, we can only align to a single neuron or discrete patterns of spiking across many neurons [5,6]. Relating activation of multiple muscles to a rich description of movement demands a reliable way to 1) reduce dimensionality when both input and output have multiple dimensions and 2) extract only the changes in movement that covary with the subset of muscles recorded.

By two-sided, we mean that PLS uses variation in both input and output to reduce dimensionality, addressing the two challenges above. Using this approach, we analyze the movement encoded in the flight muscles of the hawkmoth, Manduca sexta, to test a muscle synergy hypothesis for flight muscle coordination [13,14].

In their most general sense, muscle synergies are linear combinations of variables describing muscle activation that capture variation with fewer dimensions than the complete set of variables [9]. To avoid confusion, this use of muscle synergy differs from “infor-mation synergy” where two signals jointly provide more information than the sum of their individual contributions [21]. It is also different from the terminology of synergistic (vs. antagonistic) muscles, which refers to muscles acting on the same joint to produce movement in the same (vs. opposite) direction.

In invertebrates some synergies exist simply due to anatomy because individual motor neurons can innervate multiple muscles (e.g. [22]). Even when innervation is separate, motor units can fire in very tight synchrony, acting like a simple synergy. In Manduca, each of the main downstroke muscles (dorsolongitudinal muscles or DLMs) has five discrete subunits (Fig 1A), each innervated by a separate motor neuron (one per subunit) [13,23]. However, the whole muscle fires as one combined motor unit because the timing of the motor neurons is very precisely synchronized [23]. In fact, each DLM is driven by only a single muscle potential during each wingstroke (Fig 1B; [13,14]). Activation of the DLMs therefore varies only in the timing of the spike rather than magnitude (e.g. number of spikes). During straight flight the timing of the left and right DLM are also very precisely synchronized (spikes occur within < 1 ms of each other), but during turning the timing is modulated over an ~8 ms window [14].

Nonetheless, such small shifts produce large changes in mechanical power output [14,26]. Prior to this discovery, control of the wings was thought to be the exclusive domain of the small steering muscles that trim and tension the wingstroke [13].

Three reasons why these muscles could act as a synergy are: 1) They are activated at the midpoint of a steep but monotonic region of the power-phase curve, meaning that small changes in timing produce large, but nearly linear, changes in power output [14,26]. Therefore the timing difference between the left and right muscle’s spike translates into a power differential [14]. 2) Variation in the timing of activation is also highly correlated between the DLMs, meaning that little independent variation exists that could encode movement separately in each muscle [14]. 3) The two muscles both attach to the same large exoskeletal plate that deforms to drive wing motions, resulting in mechanical coupling of their action [24]. As we will see, the synergy hypothesis is supported when we consider wingstroke-averaged turning torque, but fails to describe the variation in motor features captured by the PLS analysis.

The coordination of the moth’s downstroke muscles is a very simple synergy hypothesis: that one variable describing the combination of these two muscles’ activities has as much predictive power as considering the two muscles independently. Manduca flight muscle is also a good system for assessing the PLS approach because the timing of activation of multiple muscles translates into a continuously varying, high dimensional pattern of torque throughout the wingstroke. Our goal is to identify the few relevant dimensions of torque that are encoded by changes in the DLMs and then to use these to assess the synergy hypothesis. Even though we focus on a specific hypothesis about an insect’s flight motor program, we use this system to show how PLS feature analysis may be generally applied to produce a data-driven decomposition of two simultaneously measured datasets.

Materials and Methods

Experimental data

It demonstrates a strong visual tracking (optomotor) response to oscillating wide-field optical patterns [27]. This behavior is ecologically relevant because moths must hover and feed while visually tracking flowers’ movement [28,29]. Flight in hawkmoths is powered by a left-right pair of DLMs and a pair of upstroke muscles (the dorsoventral muscles or DVMs; Fig 1A) [24]. However, a suite of several smaller steering muscles further trim or tension the movement of the wings, and contribute to the control of turning [13,30]. It is therefore still challenging to isolate the turning dynamics encoded in the DLMs’ activation alone.

The data used here are from the same moths used previously and under conditions elaborated on in [14]. In brief, bipolar tungsten electrodes inserted through the scutum on the dorsal surface of each moth recorded from the DLMs (Fig 1B). We tethered moths to a custom optical torque-meter that produced an output voltage dependent on the yaw (left-right) turning torque of the animal (Fig 1C) [14]. Following at least a minute of warmup shivering, moths produced full-amplitude wingstrokes spontaneously or when we elicited flight with a light touch to the neck region. The visual stimulus consisted of a sinusoidal grating of light and dark bars with a spatial frequency of 0.05 cycles degree'l, oscillating sinusoidally at 1 Hz. Because wingstroke frequency is much faster (~25 Hz), this slow variation in optic flow magnitude produced individual wingstrokes spanning a wide range of average yaw torque. We recorded electromyograms (EMGs) from the left and right DLMs and detected spikes using simple threshold crossing. Invertebrate EMGs usually afford resolution of individual muscle potentials, or spikes, whereas vertebrate recordings typically record from many motor units simultaneously, obscuring individual spikes [9,10]. Trials were analyzed further only if the right DLM’s average spike rate exceeded 18 spikes sec'l, corresponding to the lowest flight flapping frequency. To extract the torque Within each wingstroke, we had to decouple the internal dynamics of the torquemeter from the forces applied by the animal. We modeled the torquemeter as a forced, damped rotational oscillator. Where gb is the angle of rotation, I is the moment of inertia, C is the torsional damping coefficient, and K‘ is the torsional stiffness. We fit the spring, damping, and inertial parameters using a series of calibration trials following [31].

Analytical approach & synergies

During periodic movement, the motor signals form a matrix U of k timing (or phase) variables {u 1, . . ., uk}, characterizing neural or muscular action potentials (hereafter “spikes”) for each of N periods of movement (e.g. wingstrokes). During rhythmic movement with a characteristic period T, the motor output matrix is an N X 19 matrix, M, where b is the number of samples of a set of movement variables over T.

These variables were each centered and scaled by their variances. The N X 500, M movement matrix was composed of 50 ms (500 samples) long waveforms, {T1 .. . . 7500}, based on a typical tethered wingbeat frequency of 20 Hz and a 10 kHz sampling of torque. The ensemble of wingstrokes were centered and scaled by their overall variance, s. Our approach to extracting motor features follows four steps (Fig 1D—1H): 1) Alignment of the torque waveforms, 2) dimensionality reduction of the movement matrix, M, to extract a reduced basis, termed motor features, 3) synergy testing and 4) reconstruction of the torque waveforms from the motor signals, U. By comparing how well torque was encoded by different combinations of the motor signals, U, we asked if these signals act as a synergy or independently encode information about movement. In the independence model, the activation of each muscle {tL, tR} contributes significantly to predicting torque, and these two variables (the two columns of U) cannot be reduced to a single variable.

The first is based on an a priori physiological hypothesis that the timing difference At between the muscles’ spikes could translate directly into a power differential between the muscles [14,25]. We refer to this as the differential synergy model. The second synergy model is entirely data-driven: we extract the first principal component (tPCA—the PCA synergy model) of the two motor variables {tL, tR}. This closely mirrors the existing methods of muscle synergy calculation in the literature [9]. However, we use PCA instead of the sign-dependent nonnegative matrix factorization (NMF) technique [32] because timing can shift positively or negatively. We call this the empirical synergy model. Finally, we test a redundancy model, in which the variation in torque is equally well explained by only one muscle’s activation.

Using this same feature basis, we then test whether either of the synergy models can fully explain the variation in the feature basis or if the independence model (full U) is needed (Fig 1G). This avoids circularity because we identify the torque features contingent on U, but all the synergy and independence models are subsets or combinations of U. We test the different synergy, Spike-triggered ensemble time (ms) C Anima| J D Animal L extra peak doubled peak

Motor-spike-triggered ensembles and turning deciles. We aligned the torque waveforms to either the zero phase crossing (A) or the timing of the right DLM’s spike (B). We then divided the data into deciles ordered by average torque. The decile averages (C, D) and the interpolated contour surfaces (below) from two animals (“J” and “L”) show the range of torque variation within wingstrokes. Deciles are ordered from the greatest leftward (“L”) to rightward (”R”) torque. Grey transects highlight distinctive features mentioned in the text. Because wingstrokes are triggered on the right muscle’s spikes, we do not expect these patterns to be left-right symmetric. independence, and redundancy models by how well they explain variation in the projection of torque onto the feature basis (the PLS “scores”) and also how well the competing models perform in reconstruction of the entire torque waveform (Fig 1H).

Alignment of torque waveforms

We first produced an alignment that did not depend on the spikes directly. We computed the Hilbert transform of the torque signal, and filtered around the dominant 20 Hz wingstroke frequency (3—35 Hz bandpass, 8th order Type II Chebychev). The Hilbert transform returns a periodic function whose value estimates the phase of the original time series data and has been used for both gait [33] and rat whisking analyses [34]. After transforming the raw voltage signal (Fig 1C) into actual torque (Fig 1E), each torque segment was aligned to the zero phase crossing of each wing-stroke (“phase-triggered”). As an alternative, we also aligned the torque to the timing of the right DLM’s action potential (“spike-triggered”; Fig 1E). In both alignments, the resulting wingstrokes were assembled into the movement matrix, M (Fig 2A and 2B).

Dimensionality reduction & feature extraction

We first applied standard principal component analysis (PCA) by computing the eigendecomposition of the covariance matrix of M alone for both the phase and spike-trig-gered ensembles. Note that this is different from the PCA used on the motor signals, U, to create the empirical synergy model. We compared the one-sided PCA analysis of M with a two-sided, cross-covariance decomposition based on the partial least squares (PLS) method (Fig 1F). We explored how effectively these reduced descriptions of the torque output (termed motor features) captured variation in the motor signals {tL, tR} via standard regression of the features onto the signals.

PLS

PLS regression is one of several methods based on Wold’s original Projection onto Latent Structures [11], which are used widely in chemometrics [35]. PLS regression uses an iterative approach that greedily extracts features in one dataset (here M) that maximally predict the remaining variance in another dataset (here U). The set of features identified from PLS therefore are not guaranteed to maximize a global statistical property of the data [36], although each subsequent feature is optimal for that iteration. This approach has been shown to have good predictability in empirical datasets ranging from neural imaging analyses [37], ecology [38], geometric morphometrics [39,40], and paleontology [41]. We implement the faster Statis-tically-Inspired Modification of PLS (SIMPLS) developed by de long [35] rather than the original NonLinear Iterative PLS (NIPALS) [11] that requires iterative optimization steps.

This also allows the number of relevant motor output features to exceed the rank of U. The weights are the left and right singular values (SVs) from an SVD of the cross covariance matrix of M and U. Projecting these into M and U gives the loadings, which are the features and the coefficients of these projections are the scores. The SIMPLS algorithm, applied to Uand M, has the following steps: motor signals (neuromuscular) motor output (torque) Predictor var. —> M

Compute 1'1 and cl, the leading right and left singular value decompositions of So. These are the “weight” vectors of M and Urespectively. 3. From these weight vectors, the motor output “scores” vector k1 is computed as: These scores are then normalized: 4. The scores represent the amount of the leading feature represented in each wingstrokeThe motor output “loadings” vector 1)], defined as: is the first motor feature itself and is analogous to a principal component or eigenvector. 6. This leading feature is projected out of the current So matrix, which is “deflated” to obtain a new cross-covariance matrix S 1: This is the same as first estimating the regression coefficient (91 and then projecting out the newly obtained feature from M and U: and 81 can then be reconstructed from the residual M1 and U1 matrices: 7. Steps 2—5 are iterated to pull out the next set of scores, weights, and loadings {ki, di, ri, Ci, pi, qi, vi}. However, each new v,- is made orthogonal to all prior (normalized from Equation 9)

. ., k; 1}, so that only the independent variation explained in U is considered for each feature beyond the first:

However, CCA is symmetric, uses only a single cross-covariance decomposition, and only extracts a number of features up to the smaller of the rank of the inputs or the rank of the outputs (here only two). In contrast, PLS is a greedy, iterative algorithm that captures the unique contributions of successive features in describing the motor signals While preserving the asymmetry through a regression step [12]. It can therefore produce a number of features up to the number of variables describing the movement, given by the dimensionality of M.

Synergy model testing

We quantify model performance using the predicted residual sum of squares (PRESS). This form of leave-one-out cross-validation sequentially withholds each individual wingstroke from the analysis and then predicts the withheld wingstroke. If the two-variable independence model has greater predictive power than the one-variable synergy or redundancy models, then each muscle contributes significantly to the decoding of torque. If the independence model does not significantly improve PRESS, then at least one of the reduced models are favored. Note that the independence model cannot explain less variance than the synergy or redundancy models because these two alternatives are subsets of the independence model with one fewer free parameters. We used ANOVAs and paired tests to compare across models using each animal as a separate observation for model testing.

Torque waveform reconstruction

In addition to comparing how well motor signals encoded the features, we also tested how well we could reconstruct the wingstroke torques from the features (Fig 1H). The torque waveforms, M’, can be approximated from PLS features as a sum of the average motor output (the spike-triggered average—MSTA) and the features weighted by their scores:

The final score and loading vectors for the motor signal UmatriX, D = {d 1, . . ., dn} and Q = {(11, . . ., qn} define the dimensions of UmaXimally predicted in a least squares sense by each feature, pi. This approach is similar to reverse reconstruction of the stimulus given a set of spikes [5].

To demonstrate the efficacy of Wingstroke-averaged methods, we next considered the STA With just the average torque over the whole Wingstroke, < T >, reconstructed from the synergy and independence models and added as a constant offset (the mean torque models). This is the same as assuming that only Wingstroke-to-wingstroke changes are encoded in the motor signals. We then reconstructed the torque using Equation 15, Which includes the STA and the first two motor features:

Finally, we reconstructed these waveforms from the motor signals themselves by predicting the features’ scores (k 1 and k2 in Eq. 17) from the motor signals. These reconstructions were based on using the independence, redundancy, differential synergy, and empirical synergy representations of Uto generate the k’s (via the regression equations).

We also considered the residual, or unexplained, variance in the model, calculated as 1—r2, where r is the correlation coefficient. RMSE is sensitive to small, but systematic deviations (e.g. offsets) whereas uneXplained variance is sensitive to small phase shifts. Throughout the analyses we used repeated measures ANOVAs with each animal contributing an observation for each decile of turning and for each model (i.e. model and decile were each factors). Since we were primarily interested in comparing to the best model (the actual two feature waveform or the independence model), we used Hsu’s “multiple comparisons to the best” (MCB) test [43]rather than a Tukey comparison of all pairs. In paired comparisons we use paired t-tests. We confirmed that non-paramet-ric tests (Kruskal-Wallis and pairwise Wilcoxon tests) did not affect our conclusions. Statistical tests were performed in Matlab (Mathworks, Natick, MA, USA) with Hsu’s MCB tests performed in IMP (SAS Institute, Cary, NC, USA).

Cross-validation was repeated one thousand times. We also reconstructed each individual wingstroke’s torque, rather than the decile averages. In this case, the maximum reconstruction performance is likely to be limited because the motor commands from the main flight muscles should only predict a portion of the overall variation in the wingstroke. However, the ability of these motor commands to predict the scores (k,) of each PLS feature (pi) should remain high because these include only variation in torque corresponding to flight muscle variation.

Results

Aligned wingstrokes reveal patterns of torque variation

After alignment via phase or spike-triggering, we separated the resulting ensembles of wingstrokes into deciles ordered from left to right turns by mean torque. We found substantial variation, confirming our ability to visually induce a range of motor outputs (Fig 2C and 2D). While mean torque varied smoothly across deciles, there were changes in shape of the wingstroke, including phase shifts. These were stronger in some animals than in others (Fig 2D). In animal I (Fig 2C), the amplitude of the torque around ventral wingstroke reversal (~25 ms) varied and a secondary peak arose in the middle of the upstroke (~37 ms). In animal L (Fig 2D), a similar double peak formed between 50 and 80% of the wingstroke cycle and there was a prominent phase shift. These patterns were consistent for both phase and spike-triggered ensembles.

Dimensionality reduction identifies two significant features

In PCA, the spike-triggered waveforms required fewer features to reach the same explanatory power as the phase-triggered waveforms (Fig 3A and 3B). This is presumably because the spike-triggered ensemble contains more implicit information about the timing variables. PCA features are ranked in order of their ability to describe variation in the motor output. This ordering does not correspond to each feature’s ability to predict motor timing variables: some higher-ranked PCA

Identified features support the independence model. We plotted the timecourse (loadings) of each significant PLS feature (A). The decile-averaged score for each feature varied from one extreme of turning to the other (B). We determined how well different subsets of the timing variables (the independence, synergy, and redundancy models) predicted eitherthe mean torque (C) orthe feature scores (D, E) with a hold-one-out PRESS statistic. The box plots represent the mean and quartiles, with whiskers encompassing all non-outliers. Horizontal lines indicate statistically separate groups.

Some lower-rank PCA features also describe variability in the waveform ensemble that is not correlated with spike timing. As a result, the cumulative sum of the variance explained by the PCA features did not have a constant plateau (Fig 3C).

The spike and phase-triggered alignments performed comparably under PLS decomposition (Fig 3D and 3E) presumably because the spike timing information is now incorporated in the dimensionality reduction. However, with spike triggering, there was smoother accumulation of variance concentrated in the first features (Fig 3F). Including At along with tL and tR in the original U matrix did not affect the number of motor features extracted from torque (Fig 3G). The cumulative variance explained in timing variables does not reach 100% because features iteratively maximize the cross-covariance; the procedure aims to capture variance explained by the timing variables, not all variance in the output. The amount of variance explained also differs across animals. To combine all individuals, we scaled variation explained to the variance captured with 10 features (Fig 3H and 31). We did the same for the random features extracted from resampling the torque for each animal. The first two features explain significantly more variation than chance in the timing of both the left and the right muscle. While the variance generally drops off rapidly after the first feature, the second feature was important in some animals and so we retain it. Our conclusions about synergies do not depend on the inclusion of the second feature.

Features generally had four periods of oscillation per wingstroke, consistent with the mean torque waveform, or spike-triggered average (“STA”; Fig 4A). The score of the first feature correlated with the degree of turning from left to rightmost. The second feature’s score had a maximum at intermediate torque values, corresponding to straight flight, and decreased during extreme left and right turns (Fig 4B).

Varying the amplitude of the first feature has the greatest effect on torque at the beginning and midpoint of the wingstroke (Fig 4A). These points correspond to ventral and dorsal wingstroke reversals, critical moments for flight control [44]. Furthermore, the features oscillate at approximately four times the wingstroke frequency, which is consistent with the patterns observed in the aerodynamics of a robotic Drosophila wing model [45].

PLS motor features reject synergy models

We compared the variation in the torque captured by reduced combinations of the motor signals (the synergy or redundancy models) with that predicted by the two signals together (the independence model). We predicted mean torque and the identified torque features through regression and cross-validation (PRESS). When we predicted only the mean torque, but not the torque features, the one-variable synergy models performed as well as the two-variable independence models

However, when we considered either of the two PLS features, we found that the independence model significantly outperformed the redundancy or synergy models (p < 0.05 in all cases; Fig 4D and 4E). This was true even for the first feature alone, which indicates that our conclusions do not rely on the iterated extraction of features—even the first feature is encoded in both tL and tR.

Torque reconstruction improves with independently driven motor features

The first two PLS features produced reconstructions of torque that closely matched the decile means (Figs 5A, 5B and 6A), describing 95% of the variance and 70% of the residual variance after correlating the STA to the waveforms (Fig 6B). If we predicted the amount of each feature (its score) from the motor signals, the reconstruction was necessarily worse than if we used the measured feature scores (P < 0.001; Figs 5, 7A); compare cyan and green), but the resulting reconstructions were still significantly better than those based on mean torque alone (P < 0.003; Fig 7A, 7B and 7D) or using the STA without any added features (P < 0.0004; Fig 6).

The independence model accurately reconstructed 91% of the average torque waveform across all turning deciles (Fig 7B) and 53% of the residual variance after accounting for the STA. Wingstrokes during more extreme turns were less completely reconstructed (Fig 7A and 7B), but this is not surprising given that there are other muscles involved in turning whose activity we did not record. If we restrict the reconstruction to the torque projected into the two-feature subspace, the motor signals are even more accurate and consistent across deciles (Fig 5). Using the RMSE waveforms for all the deciles, we first confirmed that adding more features beyond the first two did not improve reconstruction (Fig 7C). Reconstruction based on the independence model (using both tR and tL to predict the feature scores in Equation 16) outperformed both synergy models and the redundancy model in minimizing the RMS error in the Error ratio feature-based average torque feature-based average torque models models models models

Cross-validation and sensitivity. To combine reconstruction performance into a single error metric we normalized the RMS of the errorfor each decile (averaged across animals) to the RMSE of the two feature reconstruction (*** indicates P<0.001 for all paired conditions). In addition to the standard reconstruction (A) we repeated the feature extraction and model testing with At included in the original U matrix (B) and with the phase-triggered rather than spike-triggered torque waveforms (C). Finally, we cross-validated (D) the reconstructions performing 1000 replicates with 70% of the wingstrokes as a training set and 30% withheld for a test set. Model abbreviations: 2F—two features, SAt—differential synergy, SpC—empirical (PCA) synergy, R—redundancy, l—independence, S< T > mean torque differential synergy, l< T >—mean torque independence. torque—even when controlling for the differences across deciles (effect of model in two-factor ANOVA with decile; P < 0.05 in all cases; Fig 7A and 7D). As before, we only rejected the synergy models when considering PLS features because when considering the mean torque alone, the independence and synergy models were equivalent (two-factor ANOVA; P > 0.1; Fig 7A and 7E). To summarize reconstruction performance across all deciles into a single performance measure, we took the mean ratio between the decile RMSE of each model and the RMSE of the two-feature reconstruction, which is the best two-dimensional reconstruction possible: This allows for paired tests across all animals and Wingstrokes (Fig 8A; paired t-tests comparing feature-based synergy models to independence P < 0.001 in all cases; average-torque based: P = 0.94). The conclusions held even if the most extreme left and right turning deciles

A100 \° .53; 80 '52 60 L-I—D

Reconstructing individual wingstrokes. RMS error (A) and residual variance (B) measures of the average performance for reconstructing individual wingstrokes show necessarily greater error, but consistent patterns across synergy and independence models. Colors and models correspond to those of Fig 7. Data are from reconstructing all of the individual wingstrokes from animal L (n = 298 wingstrokes) with error bars being s.e.m.

To test if the rejection of the synergy hypothesis was due to statistical bias introduced by not including a linear combination of tR and tL in the original UmatriX of motor signals, we repeated the entire analysis including At in U. The independence model still outperformed synergy models in the feature based cases (Fig 8B; feature based: P < 0.0004; average torque based: P = 0.98).

In that case, spike triggering could introduce a spike-timing dependent signal that might bias our data. It is unlikely that variation due to changes in the muscle’s timing would be phase locked independent of muscle timing, but to check for this bias we repeated the feature extraction and reconstruction using the phase-triggered torque waveforms rather than spike-triggered ensembles. The results were the same as before (Fig 8C; feature based:

First, we tested that the results were robust to cross-validation (Fig 8D). Second, we reconstructed the torque of individual wing-strokes rather than decile averages. RMS errors were naturally higher (Fig 9A), but the reconstructed waveforms were highly correlated with the measured torque (Fig 9B) and the synergy models were still rejected in feature-based analyses (paired t-test compared to independence model; all P < 10—7).

Discussion

PLS features support independent rather than synergy encoding in the DLMs

Reconstruction of decile averages (Fig 7C), cross-validated predictions (Fig 8D) and individual wing-strokes (Fig 9) are all improved by incorporating PLS features. We found that retaining the timing of both DLM activations captured more of the variance in torque than compressing these timing variables into a single linear combination, or synergy (Fig 4D and 4E). The independence model also more accurately decodes within-wingstroke torque by reconstruction (Figs 7D, 7E and 8A).

The synergy models are sufficient to account for the torque averaged over the wingstroke (Figs 4C, 7D and 8). When considering turning torque only from wing-stroke to wingstroke, the difference in timing (the differential synergy) may be sufficient to describe dynamics. However, independent variation in the two muscles is necessary to encode within-wingstroke dynamics. This does not imply that each muscle’s activation is orthogonal. These within-wingstroke dynamics are likely critical to flight control. A recent Floquet analysis treating the periodic dynamics of flapping flight rather than just the wingstroke-averaged forces demonstrated that ignoring within-wingstroke dynamics alters the stable and unstable modes of moth flight [46].

Most synergy analyses consider a large number of muscles, each described by a single, time-varying activation variable [9,10,19]. This study connects individual spike-level encoding of movement with muscle synergy questions. Additionally, while it is common to identify the force vectors or movements that correspond to muscle synergies [19,47], the extraction of synergies is usually an independent dimensionality reduction. Considering only the variation in muscle activation does not address whether or not this variation is relevant for movement generation. Here, we decompose a high-dimensional representation of torque to identify a relevant basis dependent on the variation in the muscles’ activity. We assess synergy performance against this basis.

We found that a muscle representation that retains the independent timing of both muscle activations captures more behavioral variation in torque output than three specific synergy alternatives [14]. If recordings from more muscles were included, more complex synergies might exist. For example, some of the information encoded in the timing of the left DLM might also be shared with a steering muscle’s activity. This is always a concern with encoding studies, which alone do not establish causality. Fortunately, in this case we do know that the DLMs are causally involved in producing turning torque [14]. More importantly, while we cannot say if a larger set of recordings might not express synergies, we do know that the variation in the downstroke muscles must contribute to at least two separate synergies—they cannot be expressed in a single linear dimension even with more muscles included.

However, PLS provides an alternative way of constructing synergies. We obtain not only the motor features, P, but also a matrix of equal dimension Q, whose vectors, qi, represent a nonorthogonal set of synergies for the motor signals, U (Table 1, Equation 5). This approach scales well as dimensionality increases and could also be advantageous for analyzing population encoding of complex sensory stimuli. In our case, reanalysis using the first vector q1 instead of the first PC of the DLMs’ activations does not change our conclusions.

It is possible that nonlinear analyses that take into account higher-order correlations may reach different conclusions. For example, information-based methods provide an alternative to covariance approaches that satisfy some of these same goals. The technique of maximally informative dimensions [48] seeks a feature basis which captures the most mutual information between input and output, Where the output typically is the occurrence of a single spike [49]. This approach can likely be generalized to discover complex output symbols that preserve mutual information [50]. Such techniques can handle naturalistic stimuli and minimize a priori assumptions about the statistical structure of the data, but they typically require estimating the marginal (or conditional) probability distribution of the outputs with respect to the inputs. This is very challenging When the dimensionality of both the input and output becomes large (although see [51]). In contrast, PLS readily scales to a large number of inputs and outputs because it relies only on estimating the cross-covariance structure in the data.

Motor features in the context of insect flight

DLMs can cause significant turning variation because small changes in their timing produce large changes in power [14]. However, other muscles almost certainly play a role in turning as well. In particular, small steering muscles modulate wingstroke angle (e.g. the 3rd axillary muscle) and demonstrate correlated phases of activation [13,30,52]. The role of steering muscles is one reason why the PLS approach is critical. PLS extracts features from the movement matriX, M, only if they covary with motor signals, U, and should ignore variation that is due to steering muscles but not encoded in the DLMs. Therefore, while a mechanistic understanding of turning control is not yet complete, our results do indicate that each DLM independently encodes information about the optomotor response.

Optomotor sensory input is bilateral, at least in flies, because each eye possesses distinct populations of left and right directionally selective cells [53]. However, this information could be centralized into a single descending command setting the timing difference between the flight muscles. The fact that we reject a reduced representation for the pair of DLMs indicates that visual information can separately modify the descending commands to each of the two DLMs, in addition to its known feedback to steering muscles

PLS feature analysis captures task relevant variation when multiple neuromuscular signals encode movement

By ranking features based on the cross-covariance, the PLS approach only captures motor variation relevant to the changing motor signals, While to some extent excluding components that are unrelated [12,35]. These factors account for the improved performance of PLS compared to PCA (Fig 3) and should only improve further as the dimensionality of both the motor signals and of movement increases. Our analysis considering only two muscles is the minimal case in which PLS could outperform a spike-triggered PCA. Even in this case, the abilities of PLS to deal with the challenges of multiple high dimensional datasets and incomplete representation are necessary to discover independent encoding.

Here, the approach improves predictability of wingstroke variability and our understanding of the motor program. In contrast to encoding in the central nervous system [54—56] and sensory systems, there have been few applications of dimensionality reduction approaches to encoding in the peripheral motor system. Just as sensory encoding reveals the patterns of stimuli to which neurons respond, the encoding of movement in the activation of multiple muscles reveals the structure of the motor program. New approaches like PLS-based feature analysis set the stage for understanding how the peripheral nervous system represents locomotor control.

Acknowledgments

We are grateful to Jacob Lockey and Stephanie Sundier for assisting in data collection. We thank Zane Aldworth, Iessica Fox, Jean-Michel Mongeau, and Robert Full for valuable discussions.

Author Contributions

Performed the experiments: SS. Analyzed the data: SS TLD ALF. Contributed reagents/materials/analysis tools: SS TLD ALF. Wrote the paper: SS TLD ALF.