## Abstract

Several recent studies have used matrix factorization algorithms to assess the hypothesis that behaviors might be produced through the combination of a small number of muscle synergies. Although generally agreeing in their basic conclusions, these studies have used a range of different algorithms, making their interpretation and integration difficult. We therefore compared the performance of these different algorithms on both simulated and experimental data sets. We focused on the ability of these algorithms to identify the set of synergies underlying a data set. All data sets consisted of nonnegative values, reflecting the nonnegative data of muscle activation patterns. We found that the performance of principal component analysis (PCA) was generally lower than that of the other algorithms in identifying muscle synergies. Factor analysis (FA) with varimax rotation was better than PCA, and was generally at the same levels as independent component analysis (ICA) and nonnegative matrix factorization (NMF). ICA performed very well on data sets corrupted by constant variance Gaussian noise, but was impaired on data sets with signal-dependent noise and when synergy activation coefficients were correlated. Nonnegative matrix factorization (NMF) performed similarly to ICA and FA on data sets with signal-dependent noise and was generally robust across data sets. The best algorithms were ICA applied to the subspace defined by PCA (ICAPCA) and a version of probabilistic ICA with nonnegativity constraints (pICA). We also evaluated some commonly used criteria to identify the number of synergies underlying a data set, finding that only likelihood ratios based on factor analysis identified the correct number of synergies for data sets with signal-dependent noise in some cases. We then proposed an ad hoc procedure, finding that it was able to identify the correct number in a larger number of cases. Finally, we applied these methods to an experimentally obtained data set. The best performing algorithms (FA, ICA, NMF, ICAPCA, pICA) identified synergies very similar to one another. Based on these results, we discuss guidelines for using factorization algorithms to analyze muscle activation patterns. More generally, the ability of several algorithms to identify the correct muscle synergies and activation coefficients in simulated data, combined with their consistency when applied to physiological data sets, suggests that the muscle synergies found by a particular algorithm are not an artifact of that algorithm, but reflect basic aspects of the organization of muscle activation patterns underlying behaviors.

## INTRODUCTION

The fundamental task of the nervous system is to interpret and interact with a highly complex, multidimensional environment. Both aspects of this task potentially involve monitoring the state of many thousands of variables, considering either the many individual sensory receptors or individual motor units. One strategy for the nervous system to overcome this complexity might be to identify statistical regularities within the environment and then operate using these regularities rather than the individual variables of either sensation or action. In the sensory system, such regularities might correspond to “features” in the environment, such as edges and bars present in many visual scenes and apparently encoded by neurons in visual cortices (Hubel and Wiesel 1959). In the motor system, such regularities might reflect basic biomechanical properties of the skeletomotor apparatus that can be used synergistically in the performance of different tasks (Todorov and Ghahramani 2003; Tresch et al. 2002).

One common hypothesis is that these regularities in the motor system might be represented as “muscle synergies,” each of which specifies a particular balance of muscle activations, with movements produced through the coordinated activation of these synergies (Lee 1984; Macpherson 1991; Tresch et al. 2002). As a result of producing movements through synergies, the number of degrees of freedom needed to be coordinated is substantially reduced (Tresch et al. 1999, 2002). Note that this definition of a “muscle synergy,” although common in the literature, differs from the “functional synergies” characterized by methods such uncontrolled manifold analyses (e.g., Krishnamoorthy et al. 2003). Instead, the definition of a synergy used here is identical to the “muscle modes” used in those studies. Evidence in favor of such a hypothesis has come from several studies, each examining the muscle activation patterns observed in a particular behavior and attempting to describe these activations in terms of a small number of muscle synergies (Cheung et al. 2005; d'Avella et al. 2003; Hart and Giszter 2004; Ivanenko et al. 2004; Jacobs and Macpherson 1996; Krishnamoorthy et al. 2003; Merkle et al. 1998; Olree and Vaughan 1995; Patla 1985; Ting and Macpherson 2005; Tresch et al. 1999; Weiss and Flanders 2004). In general, these studies have shown that many behaviors, ranging from withdrawal reflexes in the frog to hand manipulation in humans, can be well described in terms of combinations of a small number of muscle synergies.

Although generally supporting this hypothesis, these studies have used different methods to provide this support and to identify the synergies that might underlie behaviors. In general, these studies have used factorization algorithms, in which the observed data are modeled as a linear combination of a small set of basis vectors. Although similar to one another in their basic models, the algorithms used in these studies are based on different assumptions (Attias 1999; Basilevsky 1994; Dayan and Abbott 2001; Hyvärinen and Oja 2000; Roweis and Ghahramani 1999). Because these different assumptions will affect the performance of each algorithm, it is difficult to integrate these results. Further, although the performance of many of these algorithms is well characterized on standard data sets, their performance on data sets with properties expected for muscle activation patterns has not been systematically assessed. In particular, the fact that muscle activations, as measured in the firing rates of motoneurons or in electromyograms, are restricted to nonnegative values might impose a very strong constraint on the performance of different algorithms (Donoho and Stodden 2004; Oja and Plumbley 2004). Other characteristics of muscle activations such as the noise structure of the data (Schmidt et al. 1979; Sutton and Sykes 1967; van Beers et al. 2004) might also affect the performance of these algorithms.

For these reasons we assessed the ability of a number of different factorization algorithms to identify muscle synergies and their activation coefficients underlying a set of muscle activation patterns. We performed this analysis both on simulated data sets, for which the properties of the synergies and activations were known, and on experimental data sets, collected from behaving animals. The properties of the simulated data sets examined here were chosen so as to evaluate properties expected in physiological data sets. In addition, we examined methods to identify the correct number of synergies responsible for producing different data sets. Although we consider these issues in the particular context of the analysis of muscle activation patterns, these analyses are also more generally relevant to the application of these algorithms to nonnegative data sets in different contexts.

Preliminary results of these analyses were presented previously (Tresch et al. 2004).

## METHODS

### Generation of simulated data

We first tested the performance of each algorithm on data sets with known statistical properties. In all cases, each data set was constrained to contain only nonnegative data. Simulated data were generated as a weighted combination of basis vectors (1) where *d⃗* is an M-dimensional data vector, *w⃗ _{i}* is the

*i*th of K basis vectors also of M dimensions,

*c*is the scalar activation coefficient for the

_{i}*i*th basis vector,

*ε⃗*is an M-dimensional noise vector, and

*y⃗*=

*g*(

*x⃗*) is a thresholding function such that

*y*= 0 for

_{i}*x*< 0 and

_{i}*y*=

_{i}*x*for

_{i}*x*≥ 0. This thresholding function

_{i}*g*(

*x⃗*) enforces the nonnegative constraint. Each basis vector was scaled to have a vector norm of 1. The variables in

*Eq. 1*can be referred to differently in different contexts. In a physiological context,

*d⃗*represents the observed EMG activity for

*M*recorded muscles, each

*w⃗*represents a muscle synergy related to the synaptic weights from premotor neurons to different motoneuronal pools, each

_{i}*c*represents the synergy activation coefficient or firing frequency recruiting a synergy, and

_{i}*g*(

*x⃗*) is roughly related to the thresholding function of motoneurons. In the context of principal component analysis,

*d⃗*represents the data to be decomposed, the

*w⃗*values constitute the principal components, and the

_{i}*c*values are the component scores. Here, we will refer to the

_{i}*c*as activation coefficients and the

_{i}*w⃗*values as synergies or basis vectors. Each of these variables was manipulated within the simulated data sets to examine the performance of different algorithms. Although many different types of data sets could be examined, we focus on those that are relevant to experimental data sets.

_{i}### Simulated data sets

In all simulated data sets the number of basis vectors was four, the number of data dimensions was 12, and each data set consisted of 1,000 data points. For each type of data, 25 different individual data sets were generated according to *Eq. 1*, randomly selecting basis vectors, activation coefficients, and noise for each set.

##### W EXPONENTIAL, C EXPONENTIAL, SIGNAL-DEPENDENT NOISE.

In this data set, basis vectors and activation coefficients were drawn from exponential distributions with a mean of 10. These distributions were chosen because they were roughly similar to the distributions observed in previous analyses of experimental data. The data generated from the *w⃗ _{i}* and

*c*values were then corrupted by signal-dependent noise, with SD σ =

_{i}*a*·

*g*(

*∑*), where

_{i=1}^{N}c_{i}w⃗_{i}*a*is the slope of the relationship between the SD and the noiseless data value (Schmidt et al. 1979; Sutton and Sykes 1967; van Beers et al. 2004).

For this data set and those described below, we systematically varied the magnitude of the noise in the data set. The magnitude of this noise was calculated as 1 − *R*^{2}, where *R*^{2} is the coefficient of determination representing the percentage of variance in the noise-corrupted data set explained by linearly combining the original basis vectors and activation coefficients. We chose the slope parameter *a* so as to generate data sets with 1 − *R*^{2} levels equal to (0.05 0.15 0.25 0.35 0.45). Given that the relationship between noise magnitude and the slope of the signal dependency *a* is not straightforward, and that this relationship changes with the characteristics of each data set, we determined the slopes for each data set separately.

The properties of all data sets described below were the same as those of this data set type (i.e., exponential synergy values and activation coefficients with signal-dependent noise) unless otherwise noted. These properties were chosen as defaults because physiological signals are commonly considered to be sparsely distributed (Baddeley et al. 1997; Bell and Sejnowski 1997; Olshausen and Field 1996; Vinje and Gallant 2000) and corrupted by signal-dependent noise (Schmidt et al. 1979; Sutton and Sykes 1967; van Beers et al. 2004).

##### W EXPONENTIAL, C EXPONENTIAL, GAUSSIAN NOISE.

The properties of this data set were the same as the above default data set but with noise drawn from a constant-variance Gaussian distribution.

##### W EXPONENTIAL, C EXPONENTIAL, SIGNAL-DEPENDENT NOISE, DEPENDENCIES BETWEEN C VALUES.

For this data set, we introduced a linear dependency between the activation coefficients of basis vectors. Dependencies between activation coefficients were introduced using the following equation: *c*_{2} = *kc*_{1} + ε_{c}, where *c*_{1} and *c*_{2} are the activation coefficients of the two basis vectors, *k* is the slope of the dependency between them, and ε_{c} is noise drawn from a Gaussian distribution (variance = 9). This noise was necessary because if the two activations were perfectly correlated, any algorithm could at best identify only the correlated basis vector, *kw⃗ _{2}* +

*w⃗*, as a single vector. To ensure that the distributions of

_{1}*c*

_{1}and

*c*

_{2}were similar, for half of the activation coefficients the above equation was reversed such that

*c*

_{1}= (1/

*k*)

*c*

_{2}+ ε

_{c}. Any negative activation coefficients were set to zero. We introduced either a single dependency (with a slope of

*k*= 1; see Fig. 4

*A*) or two dependencies (with slopes of

*k*= 0.5 and

*k*= 2; see Fig. 4

*C*) between two of the four basis vectors generating a data set. The other two generating basis vectors were uncoupled.

##### W WITH NEGATIVE AND POSITIVE COMPONENTS, C EXPONENTIAL, SIGNAL-DEPENDENT NOISE.

Basis vectors were generated from exponential distributions with mean 1. Small components (<0.6) within these synergies were replaced by values resampled from an exponential distribution with the same mean and then negated, so that on average approximately half of the synergy values were negative. We generated data sets with these synergies with either no offset value or with a tonic offset value added to each muscle, mimicking a tonic activation of muscles. We chose a tonic offset so that 10% of the data values were negative and therefore thresholded (*Eq. 1*). Without the offset, an average of 49% of the data values were thresholded. Noise values were selected simultaneously with the offset to ensure 10% of values were thresholded and that the desired *R*^{2} value was achieved.

##### W EXPONENTIAL, C FROM TRUNCATED LAPLACIAN/TRUNCATED GAUSSIAN/OFFSET GAUSSIAN, SIGNAL-DEPENDENT NOISE.

For these data sets the values of the C components were drawn from different distributions. For truncated Laplacians, only nonnegative values generated from a Laplacian distribution (Dayan and Abbott 2001) with scale parameter of 10 were used for the activation coefficients. For truncated Gaussians, only nonnegative values generated from a Gaussian distribution with variance 10 were used. For offset Gaussians, values were generated from a Gaussian distribution with variance 10 and the distribution offset by its minimum so that all values were nonnegative.

### Experimental data

We also assessed the performance of each algorithm on EMG recordings obtained in behaving frogs. Muscle activation patterns were recorded during withdrawal reflexes in spinalized frogs (Tresch et al. 1999). This data set was previously used to examine the production of behaviors through combination of muscle synergies, using a variant of the nonnegative matrix factorization (NMF) algorithm (Lee and Seung 1999).

### Factorization algorithms

We compared a number of different matrix factorization algorithms. Each of these algorithms models the data according to which is similar to *Eq. 1*, except that the thresholding function, *g*(*x⃗*), is absent. Although sharing this basic model, each algorithm makes different assumptions about the properties of the elements in this equation (Attias 1999; Basilevsky 1994; Dayan and Abbott 2001; Hyvärinen and Oja 2000; Roweis and Ghahramani 1999; see discussion). Principal component analysis (PCA) was performed using the princomp function in Matlab. Note that this function uses the data covariance to identify components. In some cases, we found that the identification was slightly better using the correlation matrix, although the differences were small. Maximum likelihood factor analysis with varimax rotation (FA) was performed using the factoran function in Matlab, and wls function (weighted least squares) to estimate the factor scores. Note that the number of basis vectors that can be assessed using factor analysis is limited by the number of degrees of freedom in the model; that is, [(M − K)^{2} − (M + K)] > 0 (see Basilevsky 1994). In the data sets used here with 12 dimensions, a maximum of seven basis vectors could be used. Nonnegative matrix factorization (NMF) was implemented using the matrix multiplication update rules based on the Euclidean distance objective function described by Lee and Seung (2001). The results of this algorithm were very similar to those obtained using the gradient descent algorithm we used previously (Tresch et al. 1999) but the Lee and Seung algorithm is much more efficient. Infomax independent component analysis (ICA) was performed using the function runica in the EEGLAB package (v4.1; Bell and Sejnowski 1995; Makeig et al. 1996; http://sccn.ucsd.edu/eeglab/). We examined the performance of ICA both with the number of identified sources equal to the full number of dimensions of the data set (referred to as ICA here) and after reducing the dimensionality of the data with PCA (ICAPCA). In this latter case, the data reconstructed using the first four principal components were subjected to an ICA decomposition (note that this is essentially a post hoc factor rotation). Finally, we used a version of probabilistic ICA (pICA) that, in addition to allowing estimates of the noise for each data dimension (analogous to FA), also allowed the basis vectors and activation coefficients to be constrained to be nonnegative (Højen-Sørensen et al. 2002; Kolenda et al. 2002; http://mole.imm.dtu.dk/toolbox/menu.html). For PCA, FA, ICA, and ICAPCA, the data were zero meaned before decomposition, as required by the formulations of these algorithms.

A number of these algorithms have been used to identify muscle synergies underlying different behaviors: PCA (Krishnamoorthy et al. 2003; Olree and Vaughan 1995; Soechting and Lacquaniti 1989; Weiss and Flanders 2004), FA (Ivanenko et al. 2003, 2004; Merkle et al. 1998; Sabatini et al. 2002), ICA (Hart and Giszter 2004; Kargo and Nitz 2003), and NMF or related nonnegative factorization algorithms (d'Avella et al. 2005; Saltiel et al. 2001; Tresch et al. 1999).

In all cases, we examined these algorithms using four basis vectors to reconstruct each data set. For PCA, we took the four basis vectors that explained the largest amount of variance in the data set. For infomax ICA with dimensionality reduction (ICAPCA), we used these PCA reduced-dimension data and applied ICA to them. For infomax ICA without dimensionality reduction, we took the first four bases with the largest amount of projected mean variance (Makeig et al. 1997). For all other algorithms (FA, NMF, and pICA), the number of bases was a parameter of the model used to fit the data set and was set to four.

The algorithms described above, in addition to varying in their underlying assumptions of the data structure (see discussion for comparisons), also differ in their computational efficiency. Operating in Matlab 6.13 running on a Pentium IV (2.8 GHz, 1 Gb RAM), identification of four synergies took 0.02 s for PCA, 0.08 s for FA, 2.61 s for ICA, 0.25 s for ICAPCA, 0.05 s for NMF, and 70.17 s for pICA.

### Assessment of algorithm performance

We assessed three different aspects of algorithm performance, quantifying the ability of each algorithm to identify *1*) the subspace spanned by the original generating basis vectors, *2*) the original basis vectors, and *3*) the original activation coefficients.

To compare two sets of basis vectors (one used to generate the data set, one identified by an algorithm), we first matched vectors between the two sets. To perform this matching, we first calculated all dot products between the vectors from one set and the vectors from the other set. The two vectors in the pair with the largest absolute value of dot product were matched to one another. In cases for which the largest magnitude dot product was negative, the sign of the basis vector identified by the algorithm was reversed. The dot products between the remaining unmatched vectors were then calculated and the process repeated until all vectors were matched.

We assessed the ability of each algorithm to identify the subspace spanned by the original generating basis vectors. A set of K basis vectors defines a K-dimensional subspace out of the original M-dimensional space spanned by the muscle activations and, in principle, many different sets of K basis vectors can span this same subspace. To assess the similarity between original and identified subspaces, we calculated the principal angles (Golub and van Loan 1983) between the subspaces defined by the generating basis vectors and the basis vectors identified by each algorithm. If two subspaces are identical, then the K identified principal angles between them will be zero.

We also assessed the ability of each algorithm to identify the original basis vectors that generated the data sets. The similarity between two sets of basis vectors was taken as the average of the correlation coefficients between each of the matched basis vectors. This measure is relatively insensitive to the biases introduced by the nonnegative constraint of the NMF and pICA algorithms used here, although the results described here were generally similar when dot product or Euclidean distance similarity measures were used.

We also assessed the ability of each algorithm to identify the activation coefficients used to generate the data sets, using the correlation between the activation coefficients of two matched synergies. The ordering of activation coefficients was taken from the ordering used to match the basis vectors, as described above. These correlation coefficients were averaged across pairs within each data set.

We assessed the levels of these similarity measures expected by chance, assessed separately for each algorithm and each data set. To calculate the baseline level of these similarity measures for a particular data set, we first randomly generated 24 sets of vectors and coefficients according to the same properties used to generate the original data set. We then calculated the similarity of the identified vectors and coefficients to each of the 24 randomly generated basis vectors and activation coefficients. This process was repeated for each of the 25 repetitions for each type of data, resulting in 600 (25 × 24) values representing the case in which there was no relationship between identified and original basis vectors and activations other than that attributed to properties of the data set and biases of the algorithm. These similarities were then averaged. The baseline was then used to normalize the observed similarity value, according to *d _{s}^{norm}* = (

*d*−

_{s}*d*)/(

_{b}*d*

_{max}−

*d*), where

_{b}*d*is the normalized similarity measure,

_{s}^{norm}*d*is the nonnormalized similarity measure,

_{s}*d*is the baseline similarity measure, and

_{b}*d*

_{max}is the value for the similarity measure corresponding to perfect similarity. For the correlation coefficient,

*d*

_{max}was 1. For the subspace similarity,

*d*

_{max}was 0. These normalized measures vary between 0 and 1, with a value of 1 corresponding to maximum similarity and a value of 0 corresponding to a chance level of similarity.

For simplicity of presentation, we report only the results of applying statistical tests to these similarity measures at a single noise level, with a noise level of 15%. This level represents a case with moderate noise, similar to those expected in experimental data sets and for which most algorithms were below ceiling performance. To assess whether similarity measures for these algorithms were different from chance values, we performed two-way ANOVAs (similarity vs. algorithm and observed/chance values). In all cases, except where noted, the performance of each algorithm was better than chance for each of these three similarity measures (*P* < 0.05). To compare the performance of algorithms on a particular data set, we performed one-way ANOVAs (normalized similarity vs. algorithm). To compare the performance of algorithms between different data sets, we performed two-way ANOVAs (normalized similarity vs. algorithm and data set). Post hoc comparisons were performed using Tukey honestly significant difference corrections. Significance levels were set at *P* < 0.05.

### Choosing the number of basis vectors

In the analyses described above, we assumed that the correct number of basis vectors was known before applying each algorithm. However, in most experimental situations this is not the case. We therefore also examined methods to estimate the number of basis vectors. This problem of model order selection has been addressed by many researchers and a large number of different criteria have been proposed to choose the correct number of basis vectors (Akaike 1987; Basilevsky 1994; Minka 2000; Zucchini 2000). In general, these criteria each attempt to determine the number of basis vectors that captures the structural variation in a data set, such that any additional bases describe only unstructured noise in the data. This number is often indicated by observing changes in the slope in the plot of the criterion versus the number of bases. Although a thorough comparison of all these methods is beyond the scope of this paper, we did try several different methods, including scree plots for PCA, Bartlett's test for subsets of PCA components (barttest in Matlab), likelihood ratios for each algorithm (Basilevsky 1994), the projected variance for ICA components, the Akaike and Bayesian information criteria (AIC and BIC) for each algorithm, and the Laplacian information criterion (LIC) suggested by Minka for PCA (Minka 2000) (http://research.microsoft.com/∼minka/).

Criteria such as likelihood ratios and information criteria require likelihood calculations. For factor analysis and pICA the likelihood can be calculated straightforwardly. For the other algorithms, we followed the approach of Tipping and Bishop (1997) described for PCA. This approach cannot be used for the straight ICA model because its likelihood is defined for only the full M-dimensional data set: i.e., we can calculate the likelihood of the model identified by ICA with K = M basis vectors but not with any K < M. For ICAPCA, the likelihood can be calculated as the product of the likelihood of the ICA model applied to the K-dimensional subspace defined by PCA and the likelihood for the residual error modeled as isotropic Gaussian noise (Hansen et al. 2001). For NMF, we used the PCA formulation, estimating the observed covariance matrix as the combination of the first K basis vectors defined by NMF with isotropic residual noise. We found that pICA had difficulties identifying models with K = M bases; in this case, we used the likelihood of the model identified with K = 8 bases to be the maximal likelihood model in the likelihood ratio test, and the degrees of freedom used in the test were adjusted accordingly.

In many cases the correct number of basis vectors was indicated in the variance and likelihood curves as a change in slope (see results, Fig. 6), even though the above methods failed to identify the correct number. We thus developed an ad hoc procedure based on these explained variance and likelihood curves. For a particular data set, we applied each of the algorithms described above, varying the number of bases used in the algorithm from one to eight. We then calculated the amount of variance and the likelihood of the data set explained by each algorithm with a particular number of bases. To determine the point at which the resulting curves changes slope, we calculated the curvature for every adjacent three points in the curve (i.e., first of points 1, 2, 3; then of points 2, 3, 4; and so on). The variance and likelihood curves were first range normalized to vary between 0 and 1. We then calculated the curvature for these points and found the first resulting curvature (starting from points 6, 7, 8 and proceeding to smaller numbers), which was >0.075 and for which the corresponding three-point curve was concave (i.e., the center of the curve was located below the points in the curve). This threshold of 0.075 was arbitrarily set for all methods to give robust estimates across different factorization algorithms.

## RESULTS

### Simulated data sets

Figure 1 shows an example of the basis vectors identified by each algorithm on a simulated data set. This data set was generated with the components of the basis vectors and the activation coefficients drawn from exponential distributions and with data corrupted by signal-dependent noise (5%). The *first row* in the figure shows the four 12-dimensional basis vectors used to generate this data set. The subsequent rows show the basis vectors identified by each algorithm applied to this data set. As can be seen in the figure, FA, NMF, ICAPCA, and pICA each identified the correct bases with a high degree of fidelity. ICA identified three of the four bases very well, but misidentified the second basis as a variant of the first. The second basis vector was identified by ICA as its sixth basis vector (as ordered by their projected variance; not shown). The vectors identified by PCA, on the other hand, were different to a greater degree from, although not completely unrelated to, the original bases. The results illustrated in this figure capture the basic results of this study: that each method was capable of identifying the basis vectors underlying a set of nonnegative data, although PCA performed worse than the other methods.

Figure 2, *A–C* summarizes the performance of each of these algorithms on data sets corrupted by signal-dependent noise and with synergy values and activation coefficients drawn from exponential distributions. Performance was assessed by examining the ability of each algorithm to correctly identify the generating synergies (Fig. 2*A*), to identify the activation coefficients (Fig. 2*B*), and to identify the subspace spanned by the generating synergies (Fig. 2*C*). Each bar in these plots represents the normalized similarity averaged across 25 data sets for each of the five noise levels examined here. As suggested by the examples shown in Fig. 1, several algorithms were able to identify the correct synergies with good fidelity for data sets corrupted with small amounts of noise (leftmost bars for each algorithm in Fig. 2*A*). In particular, ICAPCA and pICA performed very well, followed by FA, NMF and ICA, then PCA (Fig. 2*A*). Not surprisingly, all algorithms were degraded with increasing noise, although notably FA was less affected than others. However, even with the highest levels of noise examined here (45%), performance of several algorithms was >0.6 (FA, ICA, NMF, ICAPCA, pICA) (Fig. 2*A*), suggesting that even with considerable noise, it is possible to estimate the synergies underlying a data set. These observations confirm the observations of Fig. 1, suggesting that most algorithms (FA, ICA, NMF, ICAPCA, and pICA) are capable of identifying a set of synergies underlying data sets.

To assess and compare the performance of these algorithms more systematically and to simplify presentation of results, we performed statistical tests on data sets with a moderate amount of noise (15%). Statistical results from this noise level were expected to be representative because performance of most algorithms with this noise level was below the ceiling, and this level is also similar to levels observed experimentally (Cheung et al. 2005).

Using these tests we confirmed that PCA performed lower than the other methods in identifying the original basis vectors as well as in identifying the activation coefficients. For synergy similarity measures, PCA was significantly below all other algorithms (*P* < 0.05). For activation coefficients, PCA was significantly lower than all but ICA (*P* < 0.05). However, PCA performed very well at identifying the correct subspace, performing at levels indistinguishable from the best of the other methods (FA, NMF, ICAPCA, pICA; *P* > 0.05) and better than ICA (*P* < 0.05).

FA performed better than PCA, worse that pICA, and was indistinguishable from the other algorithms in identifying synergies. On the other hand, it performed at levels indistinguishable from NMF, ICAPCA, and pICA (*P* > 0.05) and better than PCA and ICA in identifying activation coefficients (*P* < 0.05). In identifying the subspace, FA was better than ICA and NMF, and was not significantly different from the other methods (*P* > 0.05).

ICA performed better than PCA in identifying synergies, was indistinguishable from FA and NMF, and was below ICAPCA and pICA (*P* < 0.05). In identifying activation coefficients, ICA was lower than all methods but PCA. For the subspace identification, ICA was below all methods (*P* < 0.05).

NMF performed intermediately on all measures. It was better than PCA, indistinguishable from ICA and FA, and lower than ICAPCA and pICA in identifying muscle synergies (*P* < 0.05). In identifying activation coefficients it was better than PCA and ICA. It was better than ICA but below FA in identifying the subspace (*P* < 0.05).

Finally, ICAPCA and pICA consistently performed the best among the algorithms evaluated here across all three measures. In no case did any algorithm perform better than either ICAPCA or pICA on these measures (*P* < 0.05). Moreover, the performance of these two algorithms was indistinguishable from one another across these measures (*P* > 0.05).

In summary, the results of Figs. 1 and 2, *A*–*C* suggest that when applied to nonnegative data sets with signal-dependent noise *1*) PCA performs generally lower than the other methods, although it performs well at identifying the correct subspace; *2*) FA, ICA, and NMF perform at intermediate levels in identifying synergies, although ICA is impaired at identifying the correct subspace and activation coefficients; and 3) pICA and ICAPCA are consistently the best-performing algorithms.

In subsequent analyses, we focus on the normalized synergy similarity as a measure of algorithm performance because synergy identification has generally been the motivation for applying these algorithms to muscle activation data sets.

### Effects of noise structure: Gaussian versus signal-dependent noise

We compared the performance of each algorithm on data sets corrupted either by constant variance Gaussian noise or by signal-dependent noise. Although each algorithm examined here assumes that the data are corrupted by constant-variance Gaussian noise (see discussion), several studies have shown that muscle activations are corrupted by signal-dependent noise (Schmidt et al. 1979; Sutton and Sykes 1967; van Beers et al. 2004). We therefore assessed whether this assumption differentially affected the performance of these algorithms.

Only ICA was substantially affected by this change in noise structure, with its performance on data sets corrupted by signal-dependent noise (Fig. 2*A*) being considerably worse than its performance with Gaussian noise (Fig. 2*D*). In particular, with Gaussian noise, ICA was better than all algorithms except FA, ICAPCA, and pICA, from which it was indistinguishable, whereas with signal-dependent noise it was better than PCA, indistinguishable from NMF and FA, but below ICAPCA and pICA. The other algorithms were not significantly affected by the type of noise within the data set. ICA was also substantially impaired in identifying the correct subspace and activation coefficients with signal-dependent noise compared with its performance with Gaussian noise.

### Synergies with unequal variances

In the preceding data sets, each synergy contributed an equal amount of variance to the data set. However, in many data sets, synergies are expected to contribute unequal variances, often showing an exponential falloff in variance, as illustrated in Fig. 3*A*. We therefore assessed how each algorithm performed in identifying synergies that contributed differing amounts of variance. We simulated synergies that on average contributed 76, 18, 5, and 1% of the variance in the data set (Fig. 3*A*). As can be seen in Fig. 3*B*, there was a general trend that the performance of each algorithm decreased for synergies that contributed smaller amounts of variance to the data set. PCA showed a significantly lower performance on synergies 2–4 when compared with synergy 1 (*P* < 0.05), indicating a strong sensitivity of this algorithm to the amount of variance contributed by different synergies. The performance of ICA dropped significantly between synergies 2 and 3 (*P* < 0.05), at the point when the amount of variance contributed by the synergy was less than the noise level in the data set. FA, NMF, ICAPCA, and pICA were impaired only on the synergy that contributed the smallest amount of variance. Note also that the performance of pICA, and to a lesser extent FA, although degraded compared with their performance on the other synergies, was still quite high on the synergy with the smallest variance. These observations suggest that, although all algorithms are impaired in identifying synergies contributing small amounts of variance, some information can still be obtained about these synergies using FA and pICA.

### Introducing activation coefficient dependencies

Several of the algorithms used here assume that the data set is generated by activation coefficients with either no correlation (PCA, FA) or with complete statistical independence (ICA, ICAPCA, pICA) between them (Basilevsky 1994; Dayan and Abbott 2001; Hyvärinen and Oja 2000). NMF makes no explicit assumption about correlation. Such correlations between synergy activation sources, however, have been observed in physiological data sets (Cheung et al. 2005; Saltiel et al. 2001) and there is no a priori reason not to expect such correlations in muscle data sets. We therefore examined the effect of introducing correlations between activation coefficients on the performance of each algorithm.

Figure 4*A* shows the correlation (slope = 1) introduced between the first two synergies. The other two of the four synergies had no coupling between them. Figure 4*B* (*left*) shows the normalized synergy similarity between the coupled original synergies and the best-matching synergies identified by each method. Figure 4*B* (*right*) shows the same similarity, but for the uncoupled synergies within the same data set. As shown in this figure, all methods except PCA were impaired relative to their performance on data sets with no coupling (compare with Fig. 2*A*; *P* < 0.05). FA performed worse than all methods except PCA (*P* < 0.05), whereas the other methods were indistinguishable from one another. The performance of FA, NMF, ICAPCA, and pICA was better on the remaining, uncoupled synergies than on the coupled synergies (*P* < 0.05), and when compared with their performance on data sets with no coupled synergies (i.e., those illustrated in Fig. 2*A*), only ICA was impaired on these uncoupled synergies (*P* < 0.05). Similar results were observed for the identification of activation coefficients, although the impairment of the performance of ICA on both coupled and uncoupled synergies was greater than that of the other methods.

With the single correlation between synergies illustrated in Fig. 4*A*, each algorithm will tend to identify the coupled synergy (represented by *w⃗ _{2}* +

*w⃗*) as a single synergy within the data set. In fact, if there were no noise in the coupling of these two synergies, the data would effectively lie within a three-dimensional subspace, rather than the four-dimensional subspace spanned by the four synergies. It is thus not altogether surprising that algorithms were impaired on the correlated synergies. To examine this issue further, we introduced an additional dependency between the two synergies, as illustrated in Fig. 4

_{1}*C*. Again, the activation of the first two synergies was coupled, but this time with two different slopes, so that half of the activations was coupled with a slope of 0.5, and the other half were coupled with a slope of 2. The introduction of the second correlation will again cause the data to lie within a four-dimensional subspace and there should now be more information available to correctly identify the original synergies. Figure 4

*D*(

*left*) shows that the performance of algorithms in identifying each of the coupled synergies was considerably improved after introducing this second correlation (

*P*< 0.05), and were now indistinguishable from synergies identified for data sets with no correlations (

*P*> 0.05). Similarly, for most algorithms, their performance on the uncoupled synergies (Fig. 4

*D*,

*right*) was also improved and was not different from their performance on the data sets with no correlations (

*P*> 0.05). The notable exception to this good performance, however, was ICA, which performed substantially worse with two dependencies than it did with no dependencies (

*P*< 0.05). Inspection of the solutions found by ICA indicated that in these cases ICA identified the “coupled” synergies created by the two correlations shown in Fig. 4

*C*(i.e., 2

*w⃗*+

_{2}*w⃗*and 0.5

_{1}*w⃗*+

_{2}*w⃗*), rather than the remaining, uncoupled synergies. Of the 25 repetitions with randomly generated data sets examined with ICA, five repetitions identified both of the “coupled” synergies within the first four synergies, 17 identified one coupled synergy and one correct uncoupled synergy, and three identified both of the uncoupled synergies. It thus appears that the poorer performance of ICA for these data sets reflects a tendency for the correct, separate synergies to be identified by ICA with lower projected mean variances than the summated synergies. As expected by its degraded performance on synergy identification, a decrement of ICA's performance was also seen when assessing activation coefficient similarity.

_{1}### Muscle synergies with negative components

Another important distinction between the algorithms described above is whether they constrain the identified basis vectors to be nonnegative (NMF, pICA) or whether they allow both positive and negative components in the identified basis vectors (PCA, FA, ICA, ICAPCA). In the context of muscle activation patterns, such negative values in the basis vectors would reflect inhibitory components within the muscle synergy underlying the behavior. Such inhibitory components are obviously relevant to the production of behaviors and have been proposed to be a critical aspect of the function of spinal modules (Jordan 1991). On the other hand, the fact that the analyzed data are nonnegative potentially imposes constraints on the solutions found by any algorithm and negative values might therefore be fundamentally difficult to identify: negative values can never be observed directly in the data but only by their interactions with the positive components in other basis vectors. We therefore introduced such negative values into the generating synergies and assessed whether any algorithm was able to use these interactions to identify the correct negative values in the muscle synergies.

Figure 5*A* (*left*) shows the similarity between the generating and identified synergies for each algorithm, considering only the values identified for the negative components in the generating synergies. Note that for this measure of similarity the correlation coefficient was calculated for all negative components taken together, rather than averaging across the correlation coefficients obtained for each of the four synergies separately. This single correlation was taken because there could be cases for which there were only a few negative values in the generating synergies, and the resulting individual correlation measures would be consequently poorly characterized. As can be seen in the figure, no algorithm was able to identify the negative components in the generating synergies: PCA, FA, ICA, and ICAPCA were not better than chance levels. The positive correlations shown in the figure for NMF and pICA, even though these algorithms are explicitly constrained to identify strictly nonnegative synergy values, reflect the tendency of these algorithms to identify synergy values closer to zero when values of generating synergies are more negative. As a weaker similarity measure we also assessed how often each algorithm identified a negative value when the value of the originating synergy component was also negative. With this measure all methods except NMF and pICA gave some information about negative values, with ICA and ICAPCA performing the best (Fig. 5*B*, *left*). The performance of ICAPCA was significantly better than that of PCA and FA (*P* < 0.05) and was not different from ICA. ICA was not different from any of the other methods. Each algorithm was also impaired in identifying the correct subspace for these data sets (Fig. 5*C*) and, in fact, each was significantly worse than for data sets with no negative values (*P* < 0.05).

In the data sets used in these analyses, an average of 49% of the data values were thresholded (i.e., values of X = WC, which were negative). However, in muscle activation patterns obtained experimentally, the level of muscle activation can have a tonic offset level. Such offsets might make it possible to identify inhibition more readily because this positive offset will help the muscle activations avoid the thresholding nonlinearity. We therefore assessed the ability of algorithms to identify the negative values of synergies when the data are offset by a constant value. We chose the offset value so that 10% of data values were thresholded. As seen in Fig. 5*A* (*right*), the reduction in number of thresholded values increased the performance of each of these algorithms in identifying the negative values of synergies. This increased ability was also seen for the weaker measure of similarity assessing how often the algorithm identified negative values when the correct value was also negative, as shown in Fig. 5*B* (*right*). ICAPCA in particular performed well in identifying the negative values of synergies for these data sets.

Figure 5*D* shows the ability of these algorithms to identify the positive values of synergies when the generating synergies contained negative values. With the exception of NMF and PCA, each algorithm was impaired in identifying these positive components (*P* < 0.05). However, the similarities were still quite high, ranging between 0.6 and 0.8 for most algorithms. This suggests that, although these algorithms are generally impaired in determining inhibitory components of synergies, they are still capable of conveying substantial information about positive synergy components.

### Activation coefficients with different distributions

Several of the algorithms make explicit assumptions about the distributions from which the activation coefficients are drawn. We found here that the performance of each of these algorithms was unchanged as long as similarly sparse distributions were used to generate activation coefficients. We found that the performance of each algorithm on data sets generated with either truncated Laplacian distributions or truncated Gaussians were not different from their performance on data sets generated with activation coefficients drawn from exponential distributions (data not shown). However, when data sets were generated with activation coefficients generated from an offset Gaussian distribution, then all ICA-based methods (ICA, ICAPCA, and pICA) were impaired, as expected, based on the assumption of ICA that the sources' distribution is non-Gaussian (see discussion).

### Choosing the correct number of synergies

In all the analyses described above, we have assumed that the correct number of basis vectors was known. However, in real data sets this number must be estimated before the application of each algorithm. We therefore attempted to identify criteria that might be able to select the number of bases generating a data set consisting of nonnegative data.

As described in methods, many criteria have been proposed for this problem of model order selection. We found that for data sets generated with Gaussian noise, many of these criteria identified the correct number of bases in the simulated data for a large number of cases. In particular, Bartlett's test chose the correct number of synergies for 14/25 data sets with constant-variance Gaussian noise (noise level = 0.15), the AIC and BIC were correct in 24/25 cases using the ICAPCA algorithm, and the LIC for PCA identified the correct number in 20/25 data sets. Similarly, the likelihood ratio test using FA and ICAPCA identified the correct number in 25/25 cases; pICA, 15/25 cases; and NMF, only 3/25 cases. These results suggest that several criteria are capable of identifying the correct number of synergies for nonnegative data sets corrupted by normally distributed noise. However, for data sets generated with signal-dependent noise, we found that none of these criteria performed well: Bartlett's test, the AIC, BIC, and LIC all identified the number of synergies as being ≥10. For the likelihood ratio test, only FA was able to identify the correct number in 6/25 cases; the others identified more than eight dimensions for all data sets. In 12 additional cases, the likelihood ratio test using FA identified only one extra dimension. This lack of robust performance for most measures was true even though in many cases an indication of the correct number of synergies could be seen in the plots of either the explained variance or, more usually, of the log-likelihood versus the number of synergies.

Figure 6 shows an example of the explained variance (Fig. 6*A*) and log-likelihood (Fig. 6*B*) plots found using different algorithms when applied to one particular data set, corrupted by signal-dependent noise (15%). Note that the FA model could be identified with models with only up to seven synergies because more synergies resulted in models with more free parameters than degrees of freedom (see methods). Also note that the explained variance calculated for ICAPCA was identical to that calculated for PCA, and ICAPCA was therefore omitted from Fig. 6*A*. On the other hand, the likelihood for ICA cannot be calculated for different numbers of synergies, and ICA was therefore omitted from Fig. 6*B*. It can be seen that the explained variance curves (Fig. 6*A*) did not generally indicate a clear number of synergies: there was no sharp change in the slope of the curves at any particular number of synergies, although each curve (with the exception of ICA) approached a plateau with increasing numbers of synergies. In contrast, many of the log-likelihood curves showed a sharp change in slope. In several of the log-likelihood plots in Fig. 6*B* there is a sharp change in the slope of the curve at the point corresponding to the model fit with four basis vectors. Although the change in slope for the likelihood plots is sharper for some methods than others, this change indicates that it might be possible to extract information about the correct number of synergies from these curves that was not captured by the previously described criteria.

In particular, the results shown in Fig. 6 suggest the use of an ad hoc procedure to determine the number of synergies underlying a particular data set, based on local estimates of the curvature of the explained variance and log-likelihood curves (such as those shown in Fig. 6). Figure 7*A* shows the number of instances for which this procedure identified different numbers of synergies for data sets generated with signal-dependent noise (noise magnitude of 15%), based either on the explained variance curves (*top*) or the log-likelihood curves (*bottom*) for each algorithm. The six bars shown for each number of synergies in the figure correspond to the number of times that number of synergies was chosen using the relevant curve from each of the algorithms. As can be seen in the figure, this ad hoc procedure identified the correct number of synergies (*n* = 4) in the large majority of cases, especially when using likelihood-based measures. These results suggest that the correct number of synergies can be estimated in data sets corrupted by signal-dependent noise using such an ad hoc procedure, especially when based on log-likelihood curves.

We also examined the consequences of choosing the incorrect number of synergies. In particular, we assessed whether synergies identified in a set with *N* synergies were preserved in a set with *N* + 1 synergies (d'Avella et al. 2003). If synergies are preserved as the number of synergies are increased, then incorrectly estimating the number of synergies should not drastically alter the conclusions obtained using these algorithms. Figure 7*B* shows the average similarity between the synergies found in a set of *N* synergies and the synergies found in a set of *N* + 1. Note that for PCA and ICA, the similarity will always be 1 because the solutions found by these algorithms do not depend on the number of synergies chosen. For the other algorithms, however, it can be seen that when the number of synergies chosen is close to the correct number, synergies are well preserved between synergy sets with either one more or one less synergy. This result suggests that a slightly incorrect estimate of the number of synergies does not lead to a drastically incorrect estimate of the underlying synergies, but that features of the estimated synergies are preserved.

### Application to real data: withdrawal reflexes in the spinalized frog

To illustrate the potential use of these algorithms given the above observations on simulated data, we applied these methods to experimental data produced during withdrawal reflexes in the spinalized frog (Tresch et al. 1999). We first estimated the number of synergies. From plots similar to those illustrated in Fig. 6, we found that several variance (PCA, NMF, pICA) and log-likelihood curves (PCA, NMF, pICA) had a sharp change in slope at four synergies. Application of the ad hoc procedure to these cases confirmed this qualitative observation, identifying the number of muscle synergies to be four. The information criteria and likelihood ratio test, however, identified at least six dimensions for this data set. The FA model could not be estimated for numbers of synergies greater than five because nine muscles were used in these analyses, and it was thus difficult to estimate the number of synergies based on the curves from FA. The explained variance curve for ICA did not show a clear change in slope, similar to the curve shown in Fig. 6*A*, and although there was some indication in the likelihood curve for a change in slope at four muscle synergies, this change was not as clear as in the other methods.

Based on these observations, we applied each of the factorization algorithms to this data set using four synergies. Figure 8 shows the four muscle synergies identified by each algorithm on this data set. As in Fig. 1, the synergies identified by FA, NMF, ICAPCA, and pICA are all very similar to one another. These synergies are also all very similar to those reported previously, extracted using a nonnegative gradient descent factorization procedure (Tresch et al. 1999). Synergies 1, 2, and 4 identified by ICA are also very similar to the synergies identified by these methods, but ICA identified a different third synergy. This synergy appears to be a combination of the first and second synergies. Interestingly, the third synergy identified by the other algorithms was found by ICA as the fifth synergy, as ordered by the projected variance of each synergy. Finally, although PCA identified synergies that are generally similar to those identified by the other algorithms, there were still clear differences between the two sets of synergies. Note also that ICAPCA identified a negative ST activation in the second, third, and fourth synergies, potentially suggesting the presence of inhibition to this muscle in these synergies.

Several observations suggest that the synergies found by FA, NMF, ICAPCA, and pICA are likely to be the best estimates of the muscle synergies underlying this data set. First, these methods were generally the most robust algorithms across the range of simulated data sets studied in previous sections, especially for data sets with physiologically plausible properties such as signal-dependent noise or the presence of dependencies between activation coefficients. Finally, ICA was the only algorithm not to identify the third synergy for this data set, even though it did identify it as the fifth synergy for this data set, a pattern similar to that seen in the example of Fig. 1. This inconsistent synergy might in part reflect the impairment of ICA on data sets with signal-dependent noise. Additionally, as mentioned earlier, the third synergy identified by ICA appears to be a combination of other synergies. This type of combination was observed in the simulated data sets when dependencies between synergies were introduced (Fig. 4) and, in fact, the activation coefficients of the first and second synergies found by other algorithms (e.g., pICA) were weakly correlated (*r* = 0.31, *P* < 0.001). These observations suggest that performance of ICA in this case might have been affected by both the presence of signal-dependent noise in this data set and of correlations between the activations of different synergies. Based on these considerations, we conclude that the best estimates of the synergies underlying this data set would be those identified by FA, NMF, ICAPCA, and pICA.

## DISCUSSION

We evaluated the ability of a number of factorization algorithms to identify the basis vectors and activation coefficients underlying both simulated and experimentally obtained data sets. Taken together, these results facilitate interpretation of results of studies using different factorization algorithms, and suggest guidelines for the use of these algorithms. Importantly, the similarity of the synergies identified by these different algorithms suggests that these analyses do not produce arbitrary fits to the data but capture basic features underlying muscle activation patterns. The present study therefore suggests that factorization algorithms for the identification of the muscle synergies underlying muscle activation patterns are useful tools for examining the organization of motor behaviors.

### Assumptions underlying different factorization algorithms

In general, the results described here can be attributed to fundamental assumptions underlying each algorithm. These issues have been discussed elsewhere and we refer the reader to those studies for a more thorough discussion than that presented here (Attias 1999; Basilevsky 1994; Dayan and Abbott 2001; Hyvärinen and Oja 2000; Roweis and Ghahramani 1999). We also note that each of these algorithms, although having specific assumptions on the statistical properties of the signals, do not incorporate knowledge of the mechanical actions of muscles, such as a division of muscles into agonist/antagonist pairs. Although it might be possible and advantageous to incorporate such information into similar analyses as those presented here, at present it is not clear how this would be done.

Several of the algorithms described here differ in their assumptions on two issues: on the distributions of activation coefficients and on the noise within the data set. For the first assumption, PCA and FA explicitly assume a Gaussian distribution of activation coefficients, whereas ICA and pICA assume that the distribution is non-Gaussian. Gaussian distributions imply that the solutions found by PCA and FA can be arbitrarily rotated to produce new solutions that explain the variance in the data equally well and that have the same likelihood. To overcome this ambiguity, criteria such as the varimax criterion used here are needed as additional constraints to specify a particular rotation (Basilevsky 1994). In this context, it is interesting to note that the ICAPCA algorithm used here can be considered as a post hoc rotation of the solution found by PCA, the post hoc criterion being the independence between synergy activations. Moreover, the assumption of Gaussian or non-Gaussian distributions leads respectively to the assumption of either no first moment dependencies (PCA and FA) or complete independence between synergy activations (ICA) (Attias 1999; Basilevsky 1994; Hyvärinen and Oja 2000), potentially explaining the degraded performance of ICA when correlations were introduced between activation coefficients. Although the ICA algorithms used here explicitly assume particular distributions for activation coefficients, several studies have suggested that the particular distribution is not critical to the performance of these algorithms, as long as the distribution in the data set is as sparse as the assumed distribution (Dayan and Abbott 2001; Hyvärinen and Oja 2000). This condition was confirmed here by our simulation results from the nonnegative data sets, finding little difference in the performance of these algorithms on data sets generated with exponential, truncated Laplacian, or truncated Gaussian distributions of activation coefficients, but a more substantial decrement when the activations were drawn from offset Gaussian distributions.

There are no explicit assumptions about the distributions of activation coefficients for NMF, other than that they be nonnegative. The performance of NMF was relatively robust when the activation coefficients were drawn from different distributions. It is likely that the robust performance of NMF is to a large extent explained by the strong constraints imposed by its assumption of nonnegativity (Donoho and Stodden 2004; Oja and Plumbley 2004).

The other assumption that these algorithms differ on is the characteristics of the noise within the data set. Most algorithms assume that the data are corrupted by constant-variance Gaussian noise. PCA and ICA assume that this noise variance approaches zero; that is, they are deterministic models (Dayan and Abbott 2001; Roweis and Ghahramani 1999). The NMF update rules used here can similarly be interpreted as assuming Gaussian noise with variance approaching zero, although they can also be recast as using other noise models (Lee and Seung 1999, 2001; Cheung and Tresch, unpublished observations). FA and pICA, on the other hand, assume that the data are corrupted by nonzero noise; they are probabilistic models (Attias 1999; Dayan and Abbott 2001; Højen-Sørensen et al. 2002; Roweis and Ghahramani 1999).

The decrement of ICA for data sets with signal-dependent noise might be attributable to these assumptions about the noise structure. The relative robustness of NMF, ICAPCA, and pICA to differences in the noise structure was somewhat surprising to us, given their assumptions about noise described above. The constraints imposed by the nonnegativity condition of the NMF and pICA algorithms used here, as described above, might at least in part explain the robustness of these algorithms (Donoho and Stodden 2004; Oja and Plumbley 2004). Reasons for the robustness of ICAPCA, on the other hand, are less clear to us, although it is interesting to note that the subspace defined by PCA, on which the ICA algorithm was applied, was very close to the subspace spanned by the original basis vectors. The solution found by ICAPCA was therefore constrained to lie within the correct subspace, which might in part explain the good performance of ICAPCA. We also note that the relative ability of FA and pICA to convey information about synergies that contributed small amounts of variance to the data might reflect the fact that, alone among the algorithms tested here, these algorithms account for the noise variation differentially from structural variation in the data. As such, they are better able to find structural covariation in data sets even when the magnitude of such covariation is comparable to noise levels.

### Determining the number of synergies

We also evaluated here procedures for determining the number of synergies underlying a particular data set consisting of nonnegative values. For nonnegative data sets corrupted with Gaussian noise, a large number of model selection criteria, including information criteria and likelihood ratio tests, were capable of identifying the correct number of synergies well. Of all previously proposed criteria, only the likelihood ratio test using FA was partially able to identify the correct number of synergies in data sets corrupted with signal-dependent noise (6/25 cases). Although most of these criteria failed, an ad hoc procedure based on detecting a change in slope for either the explained variance or log-likelihood curves identified the correct number of synergies in a large fraction of cases. Although this procedure is ad hoc, it does capture the basic aspect used by several criteria in determining the correct number of synergies by searching directly for a change in slope. Its robustness in situations other than those assessed here, however, is not clear and will require further examination. It is therefore clear that this problem of model order selection remains a difficult one, and development of more principled and robust criteria for determining the number of synergies underlying a data set is an important topic of future research.

### Application of factorization algorithms to real data sets

We were surprised by the consistency between the synergies identified by many of the algorithms when applied to the experimental data set. The fact that the most consistent synergies were found by the algorithms that performed the best on simulated data sets, especially on those data sets intended to mimic physiological properties, suggests that those synergies are the best estimates of the synergies responsible for withdrawal reflexes in the frog. These observations were reassuring to us that the results previously obtained applying NMF and ICA to frog behaviors and FA to human locomotion and other behaviors most likely do not reflect the function of these particular algorithms but, instead, reflect basic aspects of the organization of behavior. As stated above, these results suggest that factorization algorithms can be profitably used to examine the production of movements through the combination of muscle synergies. A similar consistency of some of these methods (FA, NMF, ICA) across data sets has recently been reported (Ivanenko et al. 2005).

### Recommendations for application of factorization algorithms

Over all data sets examined here the best-performing algorithms were ICAPCA and pICA, with pICA consistently better than ICAPCA by a small but insignificant margin. pICA is also considerably more computationally intense than ICAPCA, requiring an order of magnitude more time to converge for a particular data set. Further, ICAPCA was the best algorithm for providing any information about inhibition in the identified muscle synergies. It would therefore seem that ICAPCA provides the best trade-off between computational efficiency and accuracy of results.

Instead of using only one algorithm, however, it would seem useful in the exploratory stage of data analyses to apply many of the algorithms described here and examine the range of solutions obtained. The solutions obtained with different methods illustrated in Figs. 1 and 8 give a good sense of the potential range of synergies underlying a data set, and the similarity between solutions, especially those found by the best-performing algorithms, is very reassuring that these are good estimates of the “correct” muscle synergies. Although it remains unclear whether this data set was in fact produced through a combination of synergies, this consistency suggests that the evaluation of this hypothesis is not critically dependent on the algorithm used.

Similarly, it would seem helpful to use several of these algorithms in determining the number of synergies underlying a data set. As illustrated in Fig. 6, although explained variance and log-likelihood curves for one algorithm might not clearly indicate a particular number of synergies, the corresponding curves for another algorithm often did.

In conclusion, we have demonstrated the ability of several factorization algorithms to robustly identify the muscle synergies and their activation coefficients within nonnegative data sets. Although these different algorithms are based on different sets of assumptions, they were surprisingly similar in the solutions they found. We therefore conclude that these algorithms can be profitably used to identify hypothetical muscle synergies underlying experimental data sets. Despite their good performance, it remains clear that the results of these analyses provide only a best estimate of the muscle synergies underlying a data set. Determining whether a particular behavior is produced through the combination of muscle synergies requires confirmation from other, additional experiments.

## GRANTS

This work was supported by the Chyn Doug Shiah Memorial Fellowship and the Schoemaker Foundation Fellowship to VCK Cheung, and National Institute of Neurological Disorders and Stroke Grants NS-09343 and NS-39865.

## Acknowledgments

We are grateful to E. Bizzi, R. Ajemian, T. Doutriaux, S. Seung and E. N. Brown for discussions of this material, and to the University of California, San Diego and Danish Technical University groups for posting proprietary algorithms on the web.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 by the American Physiological Society