Toward a unifying framework for the modeling and identification of motor primitives

A large body of evidence suggests that human and animal movements, despite their apparent complexity and flexibility, are remarkably structured. Quantitative analyses of various classes of motor behaviors consistently identify spatial and temporal features that are invariant across movements. Such invariant features have been observed at different levels of organization in the motor system, including the electromyographic, kinematic, and kinetic levels, and are thought to reflect fixed modules—named motor primitives—that the brain uses to simplify the construction of movement. However, motor primitives across space, time, and organization levels are often described with ad-hoc mathematical models that tend to be domain-specific. This, in turn, generates the need to use model-specific algorithms for the identification of both the motor primitives and additional model parameters. The lack of a comprehensive framework complicates the comparison and interpretation of the results obtained across different domains and studies. In this work, we take the first steps toward addressing these issues, by introducing a unifying framework for the modeling and identification of qualitatively different classes of motor primitives. Specifically, we show that a single model, the anechoic mixture model, subsumes many popular classes of motor primitive models. Moreover, we exploit the flexibility of the anechoic mixture model to develop a new class of identification algorithms based on the Fourier-based Anechoic Demixing Algorithm (FADA). We validate our framework by identifying eight qualitatively different classes of motor primitives from both simulated and experimental data. We show that, compared to established model-specific algorithms for the identification of motor primitives, our flexible framework reaches overall comparable and sometimes superior reconstruction performance. The identification framework is publicly released as a MATLAB toolbox (FADA-T, https://tinyurl.com/compsens) to facilitate the identification and comparison of different motor primitive models.

1. Introduction

The motor system controls a large number of degrees of freedom of the musculoskeletal system through a hierarchical architecture (Bernstein, 1947; d'Avella et al., 2015; Merel et al., 2019). It has been proposed that at the lower levels of the hierarchy, fixed modules—often referred to as motor primitives—reduce the complexity of the control problem, simplifying both motor control and learning (Tresch et al., 1999; Flash and Hochner, 2005; Bizzi et al., 2008; Bizzi and Ajemian, 2020; Cheung and Seki, 2021). Convincing evidence for the existence of motor primitives is provided by the consistent observation of temporal and spatial regularities at different levels of the control hierarchy. For example, such regularities have been reported at the level of the motor cortex (e.g., Overduin et al., 2015; Kadmon Harpaz et al., 2019), spinal interneurons (e.g., Hart and Giszter, 2010; Levine et al., 2014; Takei et al., 2017), motor neurons (e.g., d'Avella et al., 2003; Ivanenko et al., 2004; Torres-Oviedo et al., 2006), joint kinetics (e.g., Mussa-Ivaldi and Giszter, 1992; Santello and Soechting, 2000; Thomas et al., 2005), and joint kinematics (e.g., Santello et al., 1998; Kaminski, 2007; Chiovetto and Giese, 2013).

However, the heterogeneity of domains and observation levels has led to the development of a variety of computational models to explain how the observed regularities are plausibly generated by the underlying fixed modules. Such variety, in turn, has created the need to use different identification algorithms that sometimes have to be devised ex novo. For example, motor primitives at the electromyographic (EMG) level are typically extracted with non-Negative Matrix Factorization (NMF—e.g., Tresch et al., 1999; Ting and Macpherson, 2005; Godlove et al., 2016), while motor primitives at the kinematic level are often extracted with Principal Component Analysis (PCA—e.g., Santello et al., 1998; Kaminski, 2007; Chiovetto et al., 2012), Independent Component Analysis (ICA—e.g., Mori and Hoshino, 2002; Lambert-Shirzad and Van der Loos, 2017), or Factor Analysis (FA—e.g., Smith et al., 2006; Steinberg and Bock, 2013). Alternative approaches are also common: for example, ICA has also been successfully used to extract motor primitives at the EMG level (e.g., Hart and Giszter, 2004; Ivanenko et al., 2005; Dominici et al., 2011). More complex models, which introduce trial-dependent delays in the motor primitives activation (Omlor and Giese, 2006) and try to simultaneously capture both temporal and spatial primitives (d'Avella et al., 2003; Delis et al., 2014), require specialized identification algorithms. This multitude of mathematical models and identification algorithms complicates the comparison of the results from different studies: in the absence of a standardized framework for the definition and identification of motor primitives, it is hard to assess whether potential differences observed between studies are due to the use of a different model, identification algorithm, or genuine experimental manipulations.

To simplify such comparative analysis, we introduce here a new unifying framework to model and identify several popular classes of motor primitive models. Specifically, we first show that common models of spatial, temporal, and spatiotemporal modularity, with and without delays, can be considered as special cases of a single generative model: the anechoic mixture model. We then introduce a new class of identification algorithms, which we derived extending the Fourier-based Anechoic Demixing Algorithm (FADA—Chiovetto and Giese, 2013) to fit all the considered modularity models. Finally, we validate our framework by showing that it can robustly extract different classes of motor primitives from both simulated and experimental data with an accuracy that is comparable and sometimes superior to that achieved using model-specific identification methods.

2. Methods

2.1. Generative models of spatial and temporal regularities

This section provides a brief survey of the most common models of the modular organization of motor behavior. In general, such models explain the spatial and temporal invariances observed during movements as arising from spatial and temporal modules that are fixed across trials. Such modules are typically referred to as motor primitives or synergies. In the following, we will assume that the activity patterns of M degrees of freedom (DOFs) recorded during the execution of one of L different trials, over T time samples, are collected in an M by T matrix X^l, where l is the trial index. Depending on the model, it will sometimes be useful to refer to individual column vectors x^l(t), or to individual entries $x_m^l(t)$ of this matrix. Depending on the context, the DOFs represent different electromyographic (EMG), kinetic, or kinematic signals. Signals relative to different trials are considered to have a fixed duration T_s and to be sampled with a constant sampling frequency.

2.1.1. The spatial decomposition model

One classical definition of motor primitives is based on the observation, during rhythmic and goal-directed movements, of specific covariation patterns between different DOFs that are invariant across time and trials. Such fixed covariation patterns are typically interpreted as reflecting the coordinated recruitment of multiple muscles or joints. This type of model has been particularly successful at explaining regularities in EMG signals (Tresch et al., 1999; Ting and Macpherson, 2005; Torres-Oviedo et al., 2006). Consistent with these observations is the following generative model:

$\begin{array}{ll}{\text{x}}^l(t)={\displaystyle \sum _{p=1}^P}{\text{w}}_p·c_p^l(t)+residuals& (1)\end{array}$

In this equation, the vectors x^l(t) collect the values taken on by all the DOFs at time point t (assuming discrete time steps, 1 ≤ t ≤ T) in trial number l. The column vectors w_p capture the invariant spatial patterns and thus represent the motor primitives themselves. The number of primitives is P, and the scalars $c_p^l(t)$ are the time-dependent mixing weights of the primitives. Mixing weights (and residuals) can generally vary across trials. Importantly, the components of this model are often assumed to be non-negative [i.e., $c_p^l(t)\ge 0$ and w_p,m ≥ 0]. This assumption is particularly common when the model is used to explain EMG data, which typically consist of time series of non-negative signals [i.e., $x_m^l(t)\ge 0$, ∀m, l] reflecting the excitatory activity of underlying motoneurons (Farina et al., 2014). As this model is based on the invariant patterns in the spatial domain (i.e., in the DOF space), it is often referred to as spatial decomposition model.

2.1.2. The temporal decomposition model

An alternative definition of motor primitives is based on the observation of invariant covariation patterns across time, which are thought to represent the activity of latent temporal source functions s_p(t). Temporal components based on this definition have been identified in kinematic (Kaminski, 2007; Berret et al., 2009; Chiovetto and Giese, 2013), kinetic (Thomas et al., 2005), and EMG (Ivanenko et al., 2004, 2005; Chiovetto et al., 2010) space. The underlying generative model, which we will refer to as temporal decomposition model is defined by:

$\begin{array}{ll}{x}_m^l(t)={\displaystyle \sum _{p=1}^P}{c}_{mp}^l·s_p(t)+residuals& (2)\end{array}$

In this equation, $x_m^l(t)$ is the value of the m-th DOF at time t in trial number l, and the corresponding scalar mixing weights $c_{mp}^l$ change between different trials. The P temporal primitives s_p(t), however, are assumed to be invariant over trials. Both the spatial (1) and the temporal (2) decomposition models assume that the latent sources affect the activity patterns of all the different degrees of freedom simultaneously, that is, without any DOF-specific delays. For this reason, such models are sometimes referred to as synchronous decomposition models.

2.1.3. The temporal decomposition model with delays

An alternative to the synchronous decomposition models to explain temporal regularities has been proposed by Omlor and Giese (2006, 2007, 2011). This model allows for delayed activation of the temporal basis functions, where the delays can potentially vary across different primitives, DOF, and trials. This can be interpreted as reflecting, for example, delays between the activation of different muscles within the same temporal primitive. Mathematically, this model is characterized by the equations:

$\begin{array}{ll}{x}_m^l(t)={\displaystyle \sum _{p=1}^P}{c}_{mp}^l·s_p(t-{\tau }_{mp}^l)+residuals& (3)\end{array}$

Importantly, in this model, the time delays ${\tau }_{mp}^l$ and mixing weights $c_{mp}^l$ can vary over trials, while the basis functions s_p(t) are invariant, as in model (2). Like for model (1), inequality constraints can be imposed on the mixing weights and source functions of models (2) and (3) to account for the non-negativity of EMG signals.

2.1.4. The spatiotemporal decomposition model

The models discussed so far, defined by the Equations (1)–(3) can only account for regularities in the spatial or temporal domain, but not both. To deal with such a limitation, d'Avella and Tresch (2002), d'Avella et al. (2003, 2006) introduced the concept of spatiotemporal (or time-varying) primitives, which can be considered as latent spatiotemporal activity patterns that are invariant over trials. The resulting spatiotemporal decomposition model thus assumes that the activity patterns measured on the DOFs are generated by mixing such primitives with trial-specific weights. Such latent sources, or primitives, can also be shifted in time by a trial-specific delay. This results in the following generative model:

$\begin{array}{ll}{\text{x}}^l(t)={\displaystyle \sum _{p=1}^P}{c}_p^l·{\text{w}}_p(t-{\tau }_p^l)+residuals& (4)\end{array}$

Note that in this model, unlike in model (3), mixing weights $c_p^l$ and delays ${\tau }_p^l$ do not change across muscles. The vector-valued source functions w_p(t) are invariant across trials and represent the spatiotemporal primitives. Such primitives and the corresponding mixing weights have typically been assumed to be non-negative (d'Avella et al., 2003), although also models with unconstrained parameters have been applied to model phasic EMG activity (d'Avella et al., 2006).

2.1.5. The space-by-time decomposition model

An alternative approach to simultaneously model spatial and temporal regularities was introduced by Delis et al. (2014). This model, named space-by-time decomposition model, assumes the simultaneous existence of P_sp spatial and P_tp temporal primitives that, unlike in model (4), are not grouped in P spatio-temporal primitives, but are free to vary independently.

$\begin{array}{ll}{\text{x}}^l(t)={\displaystyle \sum _{p=1}^{P_{tp}}}{\displaystyle \sum _{q=1}^{P_{sp}}}{s}_p(t-{\tau }_{pq}^l)·c_{pq}^l·{\text{w}}_q+residuals& (5)\end{array}$

In this model, w_q and s_p(t) define the trial-independent spatial and temporal primitives as in models (1) and (2), while the mixing weights $c_{pq}^l$ and time delays ${\tau }_{pq}^l$ are trial-dependent. Since the model was originally designed to account for EMG data, all parameters are typically assumed to be non-negative, with the exception of the time delays.

2.2. The unifying anechoic mixture model

All previously discussed models can be derived as special cases of a single model: the anechoic mixture model. This type of model is popular in acoustics, where it is applied for modeling acoustic mixtures in reverberation-free rooms (Torkkola, 1996; Emile and Comon, 1998; Bofill, 2003; Yilmaz and Rickard, 2004). This model assumes typically a set of R recorded acoustic signals y_r(t) that are created by the superposition of U acoustic source functions f_u(t), where time-shifted versions of these source functions are linearly superposed with the mixing weights a_ru. The time shifts are given by the time delays τ_ru. This models the fact that for a reverberation-free room the signals from the acoustic sources arrive at the receiver with different time delays and attenuated amplitudes, which are dependent on the distances between the acoustic sources and the receivers. The corresponding generative model has the following form (for 1 ≤ r ≤ R):

$\begin{array}{ll}{y}_r(t)={\displaystyle \sum _{u=1}^U}{a}_{ru}·f_u(t-{\tau }_{ru})+residuals& (6)\end{array}$

2.3. The anechoic mixture model subsumes previous modular models of movement generation

In this section, we show that all the models of spatial, temporal, and spatiotemporal modularity discussed so far can be considered as a special case of the anechoic mixture model (6).

2.3.1. Derivation of the spatial decomposition model

Identifying the signals of type y_r(t) with the components of the vectors x^l(t), i.e., $y_r(t)=x_{m(r)}^{l(r)}(t)$ where the integer functions l(r) and m(r) define a one-to-one mapping between the m-th degree of freedom in trial l and the corresponding signal y_r(t) (with 1 ≤ r ≤ M · L), and constraining the time delays τ_ru to be zero, one obtains the model (1). Since in this model the weight vectors w_p are assumed to be invariant over trials, all mixing weights a_rp belonging to the same DOF and primitive number P have to be equal and independent of the trial number, so that a_rp = w_p,m(r), where the function m(r) returns the number of the DOF that belongs to index r independent of the trial number. The time-dependent mixing coefficients $c_p^l(t)$ of the model (1) correspond to the source functions f_u of the model (6), thus $f_u(t)=c_{p(u)}^{l(u)}(t)$ where here the index u runs over all combinations of the indices p and l, thus 1 ≤ u ≤ U = P · L and where the integer functions l(u) and p(u) establish mappings between the number of the source function in model (6) and the time-dependent mixing weights in model (1). Non-negativity constraints can be added for the model parameters a_rpand the functions f_u(t), e.g., for the modeling of EMG data.

2.3.2. Derivation of the temporal decomposition models

If one identifies the source functions in model (6) with the temporal primitive functions s_p(t), i.e., f_p(t) = s_p(t), 1 ≤ p ≤ P and again constrains the delays τ_ru to be zero, Equation (6) becomes equivalent to model (2). In this case, the mixing weights a_rp are equated to the mixing coefficients $c_{mp}^l$ in model (2), where the index r runs over all combinations of m and l, formally $a_{rp}=c_{m(r),p}^{l(r)}$, with appropriately chosen integer functions m(r) and l(r). Like for model (1), the components of the data vector have to be remapped over DOF and trials according to the relationship $y_r(t)=x_{m(r)}^{l(r)}$(t). Again, non-negativity constraints can be added for the parameters a_rp and to the source functions f.

Additionally, assuming that the delays of the anechoic model can take on any value (τ_ru ≠ 0), and equating the delays in model (3) according to the relationship ${\tau }_{rp}={\tau }_{m(r),p}^{l(r)}$, makes model (6) equivalent to model (3).

2.3.3. Derivation of the spatiotemporal decomposition model

Introducing individual sets of basis functions for the different DOFs, grouping them into vectors and equating the mixing weights and temporal delays for the components of each vector, transforms model (6) into the model (4). On the level of the time-dependent basis functions, this equivalence can be mathematically described by the equation f_u(t) = w_{p(u), m(u)}(t), where w_{p, m} corresponds to the component of the basis function vector w_p(t) that belongs to the m-th DOF, and where the integer functions m(u) and p(u) establish a one-to-one mapping between the indices of the basis functions in the two models and the number of the associated DOF. This assignment is independent of the trial index l. The index r in (6) runs over all combinations of DOF and trial numbers, thus 1 ≤ r ≤ M · L. The integer functions m(r) and l(r) assign the corresponding trial number and DOF to the index r in the model (6). Thus, the assignment equation for the data vector is again given by $y_r(t)=x_{m(r)}^{l(r)}(t)$ for the m-th DOF in the l-th trial. The requirement that all mixing weights and temporal delays belonging to the same basis function vector w_p are equal is equivalent to a set of equality constraints, which can be captured by the equation systems $a_{ru}=c_{p(r)}^{l(r)}$ and ${\tau }_{ru}={\tau }_{p(r)}^{l(r)}$. Again, non-negativity constraints can be added, if necessary.

2.3.4. Derivation of the space-by-time decomposition model

In order to establish equivalence with the model (5), the data vectors of the models are mapped onto each other according to the relationship $y_r(t)=x_{m(r)}^{l(r)}(t)$, where again l(r) and m(r) are integer mapping functions that assign the r-th element of the data vector of the model (6) to the m-th DOF of the data vector x^l for the l-th trial in (5) with 1 ≤ r ≤ M · L. Model (5) has a total of P_sp·P_tp temporal basis functions, where however the functional forms of the basis functions for different indices q (i.e., different spatial components) for the same p (i.e., same temporal component) just differ by time shifts. This is equivalent to an equality constraint for these functions, which can mathematically be characterized in the form f_u(t) = s_p(u)(t), with 1 ≤ u ≤ P_tp and the index functions p(u) and q(u) that map the index u in the model (6) onto the indices of the temporal and spatial primitive in (5). Since all indices with the same p(u) are mapped onto the same basis function s_p, the last equation specifies an equality constraint. With the same integer mapping functions, finally, also the relationship between the mixing weights can be established, which is given by the equation $a_{ru}=c_{p(u),q(u)}^{l(r)}·w_{q(r),m(u)}$, where w_q,m is the m-th element for the vector w_q. The last equation specifies a bilinear constraint for the weight parameters of the model (6). Using the same notation, the equivalence between the delays is established by the equation system ${\tau }_{ru}={\tau }_{p(u),q(u)}^{l(r)}$. A summary of the established equivalences between the general model (6) and the other models is given in Table 1.

TABLE 1

Table 1. Constraints that make the primitive models (1), (2), (3), (4), and (5) equivalent to the general anechoic model (6). See text for details.

2.4. FADA: An efficient algorithm for the identification of motor primitives within the unified framework

All algorithms for blind source separation require the identification of a large number of parameters. For example, a complete definition of model (6) requires the specification of: T · U parameters to represent the sources f_u(t), R · U parameters to represent the activation weights a_ru, and R · U parameters to represent the delays τ_ru. Thus, fitting a dataset with L trials and M degrees of freedom (i.e., with R = M · L), requires the estimation of a total of (T + 2M · L) · U parameters. In typical settings, the number of time samples T is much larger than the number of degrees of freedom M, sources U, and trials L. For applications in motor control, the relevant signals are subject to additional constraints, which can be exploited for the derivation of more efficient algorithms. Signals in motor control are typically smooth. This allows to reduce the complexity of the anechoic demixing problem considerably and to devise algorithms that are more robust than those developed for general purposes.

Here, we first present a robust algorithm for standard anechoic demixing, which can be used for the identification of the parameters associated with the unconstrained model (6). This algorithm, which relies on the representation of the signals in the frequency domain and is thus called Fourier-based Anechoic Demixing Algorithm (FADA), was introduced in a previous study (Chiovetto and Giese, 2013) to identify the temporal decomposition model (3). In the present work, we describe how this algorithm can be extended by inclusion of additional constraints to make it suitable for the identification of the parameters associated with all the most common models of motor modularity (1–5). The collection of algorithms for the estimation of all models of motor modularity is released as a public toolbox: the FADA toolbox (FADA-T, https://tinyurl.com/compsens).

The typical smoothness of the activity patterns (e.g., joint trajectories, EMG envelopes, etc.) recorded during the execution of body movements, implies that they can be well described by mixtures of smooth source functions (Chiovetto and Giese, 2013). These smooth source functions can, in turn, be well-approximated by truncated Fourier expansions, defined by only K nonzero complex Fourier coefficients, where K is typically far below the Nyquist limit (K ≪ T/2). Consequently, the number of parameters to identify to describe the original dataset drops to (K + 2M · L) · U. This decreases substantially the computational costs of the parameter estimation problem and makes it less prone to premature convergence to local minima.

Band-limited temporal signals y_r(t) and source functions f_u(t) can thus be approximated by truncated Fourier expansions of the form:

$\begin{array}{ll}{y}_r(t)={\displaystyle \sum _{k=-K}^K}{c}_{rk}{e}^{\frac{2\pi ikt}{T_s}}& (7)\end{array}$

and

$\begin{array}{ll}{f}_u(t-{\tau }_{ru})\cong {\displaystyle \sum _{k=-K}^K}{\nu }_{uk}{e}^{-ik{\tau }_{ru}}{e}^{\frac{2\pi ikt}{T_s}}& (8)\end{array}$

where c_rk and ν_uk are complex constants (i.e., $c_{rk}=\left|c_{rk}\right|e^{i{\phi }_{c_{rk}}}$ and ${\nu }_{uk}=\left|{\nu }_{uk}\right|e^{i{\phi }_{{\nu }_{uk}}}$, where i is the imaginary unit, and φ_crk and φ_νrk are real numbers). The positive integer K is determined by Shannon's theorem according to the limit frequency of the signals, and T_s is the temporal duration of the signal. The source separation algorithm tries to ensure that the source functions f_u(t) are uncorrelated over the distributions of the approximated signals. This implies $E\left\{f_u(t)·f_{u^{\prime }}(t^{\prime })\right\}=0$ for u ≠ u′ and any pair t ≠ t′. For the corresponding Fourier coefficients this implies $E\left\{{\nu }_{uk}·{\nu }_{u^{\prime }{k}^{\prime }}\right\}=0$ for u ≠ u′ and any pair k ≠ k′. Combining Equations (6), (7), and (8) we obtain by comparison of the terms for the same frequency

$\begin{array}{ll}{c}_{rk}={\displaystyle \sum _{u=1}^U}{a}_{ru}·{\nu }_{uk}{e}^{-ik{\tau }_{ru}}& (9)\end{array}$

From this follows with $E\left\{{\nu }_{uk}·{\nu }_{u^{\prime }{k}^{\prime }}^{*}\right\}=E\left\{|{\nu }_{uk}{|}^2\right\}·{\delta }_{u{u}^{\prime }}$ the equation:

$\begin{array}{l}|c_{rk}{|}^2=E\{\left|c_{rk}\right|\}\\ ={\displaystyle \underset{u=1}{\overset{U}{{{\displaystyle \sum }}^{\text{}}}}}{\displaystyle \sum _{u^{\prime }=1}^U{a}_{ru}{a}_{r{u}^{\prime }}E\{{\nu }_{uk}·{\nu }_{u^{\prime }{k}^{\prime }}^{*}\}{e}^{-ik({\tau }_{ru}-{\tau }_{r{u}^{\prime }})}}\\ ={\displaystyle \underset{u=1}{\overset{U}{{{\displaystyle \sum }}^{\text{}}}}}{a}_{ru}^{2}E\left\{\right|{\nu }_{uk}{|}^2\}\\ ={\displaystyle \underset{u=1}{\overset{U}{{{\displaystyle \sum }}^{\text{}}}}}\left|a_{ru}{|}^2\right|{\nu }_{uk}{|}^2& (10)\end{array}$

The symbol * indicates the conjugate of a complex number. The derivation of this equation replaces the expectations of the Fourier coefficients c_rk with their deterministic values and treats the source weights a_rk as deterministic trial-specific variables. This can be justified if these mixture weights are estimated separately from the sources in an EM-like procedure. Empirically, however, we obtain reasonable estimates of the model components based on Equation (10) also using other methods (see below). Since the signals f_u(t) and y_r(t) are real the corresponding Fourier coefficients fulfill $c_{rk}=c_{r,-k}^{*}$ and ${\nu }_{uk}={\nu }_{u,-k}^{*}$. Thus, the demixing problem needs to be solved only for parameters with k ≥ 0.

The previous considerations motivate the following iterative algorithm for the identification of the unknown parameters in model (6). After random initialization of the parameters to be estimated, the following steps are carried out iteratively until convergence:

1. Compute the absolute values of the coefficients c_rk and solve the following equations:

$\begin{array}{ll}|c_{rk}{|}^2={\displaystyle \sum _{u=1}^U}|a_{ru}{|}^2|{\nu }_{uk}{|}^2& (11)\end{array}$

with r = 0, 1, … R and k = 0, 1, … K. In our study we exploited non-negative ICA (Højen-Sørensen et al., 2002) to solve this equation. In the distributed version of the software, the Equation (10) can also be solved via non-negative matrix factorization (Lee and Seung, 1999, 2000).

2. Initialize for all pairs and iterate the following steps:

(a) Update the phases of the Fourier coefficients of the sources, defined as φ_νuk = angle(ν_uk) = arctan[Im(ν_uk)/Re(ν_uk)] by solving the following non-linear least square problem

$\begin{array}{ll}{\displaystyle \underset{\Phi }{\text{min}}}\parallel \text{C}-\text{Z}(\Phi ){\parallel }_F^2& (12)\end{array}$

where (C)_rk = c_rk, ${(\text{Z})}_{rk}={\displaystyle \sum _{u=1}^U{a}_{ru}{e}^{-ik{\tau }_{uk}}\left|{\nu }_{uk}\right|e^{i{\phi }_{{\nu }_{uk}}}}$ and ∥∥_F indicates the Frobenius norm. In order to avoid cluttered notation, for the function Z(.) only the arguments with relevance for the optimization are explicitly written.

(b) Keeping the identified source functions f_u(t) constant, identify for each signal y_r(t) the weights a_ru and delays τ_ru by minimization of the following cost functions:

$\begin{array}{ll}{\displaystyle \underset{{\text{a}}_r,{\tau }_r}{argmin}}\parallel y_r(t)-\text{f}{(t,{\tau }_r)}^{\prime }{\text{a}}_r{\parallel }_F^2& (13)\end{array}$

where ′ is the transposition operator. Optimization with respect to a_r and τ_r is feasible, assuming uncorrelatedness of the functions f_u and independence of the time delays (Swindlehurst, 1998). The column vector a_r concatenates all weights related to DOF r, i.e., ${\text{a}}_r={\left[a_{r1},\dots ,a_{rU}\right]}^{\prime }$. The vector function ${\text{f}}_r(t,{\tau }_r)={\left[f_1(t-{\tau }_{r1}),\dots ,f_U(t-{\tau }_{rU})\right]}^{\prime }$ concatenates the functions f_u, shifted by the time delays related to the r-th DOF.

The original version of the FADA algorithm was designed to solve the source separation problems without constraints. Additional constraints, such as the non-negativity of the parameters or additional equality constraints for the weights and delays can be easily added, due to the modular structure of the algorithm. A detailed description of how we imposed such constraints can be found in the Appendix. Importantly, the addition of such constraints allowed us to develop a unifying method to identify the parameters of all the considered models of motor modularity.

2.5. Validation of the proposed framework

To validate our framework, we used both biologically realistic simulated data and real experimental data. The simulated data included a total of 20 noisy realization from the models (1–5), for five different noise levels. Importantly, for models (1–3), we generated data from both the unconstrained and the non-negative variants (Figure 1). This allowed us to simulate kinematic-like and EMG-like signals from a total of eight different models. We used the data generated from these models to benchmark the ability of FADA-T to identify the ground-truth generative models against that of other standard model-specific identification methods. These methods included: the fast Independent Component Analysis (fastICA—Hyvärinen and Oja, 1997), the Non-negative Matrix Factorization algorithm (NMF—Lee and Seung, 1999) the Anechoic Demixing algorithm (AnDem—Omlor and Giese, 2011), the Shifted ICA (SICA—Mørup et al., 2007b), the Anechoic NMF algorithm (AnNMF—Omlor and Giese, 2011), the shifted NMF (sNMF—Mørup et al., 2007a), the spatiotemporal NMF (stNMF—d'Avella et al., 2003), and the sample-based Non-negative Matrix tri-Factorization algorithm (sNM3F—Delis et al., 2014).

FIGURE 1

Figure 1. Representative primitives generated to assess the identification performance of FADA-T. (A) Unconstrained spatial primitives, associated with model (1). (B) Non-negative spatial primitives, associated with model (1) and (5). (C) Unconstrained temporal primitives, associated with models (2) and (3). (D) Non-negative (EMG-like) temporal primitives, associated with models (2), (3), and (5). (E) Spatiotemporal non-negative primitives, associated with model (4).

The experimental data included a dataset of body joint kinematics recorded during the execution of emotional walks (Omlor and Giese, 2007; Roether et al., 2009; Endres et al., 2013), and a dataset of arm EMG signals recorded during the execution of reaching movements (d'Avella et al., 2006). We fitted the temporal decomposition model with delays (3) to the kinematic dataset, and the spatiotemporal decomposition model (5) to the EMG dataset.

To assess the performance of the identification algorithms, we considered two measures. Specifically, as a measure of the ability of the algorithms to retrieve a solution consistent with the observed data, we computed the fraction of explained variance R². As a measure of the ability of the algorithms to identify the ground-truth primitives, activation coefficients, and delays, we computed the normalized similarity S_N between the identified and ground-truth quantities. To normalize such measures, we computed a baseline similarity measure, which estimates the average similarity between random pairs of realizations of a single model. The normalized similarity S_N takes on values between zero and one, where zero indicates random estimates. Further details about the procedure we used to generate biologically realistic signals, the experimental dataset, the benchmark model-specific identification algorithms, and the similarity measures, can be found in the Appendix.

2.6. Statistical analyses

All tested measures were normally distributed according to a Chi-square goodness-of-fit test. Student's t-test was used to test whether the reconstruction and similarity measures were statistically different from the chance level. Group differences were statistically tested by two-way ANOVAs with Algorithm and Noise Level as factors. When appropriate, we performed post-hoc analysis with the Tukey-Kramer test. As a level of significance for the rejection of the null hypotheses we chose α = 0.05.

3. Results

3.1. Evaluation of algorithm performance on simulated data sets

To assess the identification performance of FADA-T, we generated EMG-like and kinematic-like ground-truth data, based on the mixture models defined by the Equations (1), (2), (3), (4), and (5). The aim of our comparison was to assess whether FADA-T could identify mixture parameters at least as well as other established model-specific unsupervised learning methods.

Figure 2A shows the average performance (±SD) of the FADA-T and the fastICA algorithms (Hyvärinen and Oja, 1997) on the identification of the parameters of the spatial decomposition model (1). The bar plots represent the reconstruction accuracy (R²) and the normalized similarity (S_N) between the ground-truth and the extracted primitives and weights, averaged across 20 realizations, for five different levels of signal-dependent noise. Asterisks indicate significant differences between algorithms, according to post-hoc testing. Qualitatively, both algorithms provided a good level of reconstruction accuracy and resulted in an accurate estimation of the original model parameters. Accuracy measures were typically larger than 0.5, and the similarity measures for the recovered primitives and weighting coefficients were always significantly larger than chance [t₍₁₉₎ > 9.93, p < 0.001]. The two-factor ANOVAs revealed a significant main effect of the Noise factor on both the reconstruction accuracy and the identification of the primitives [F_{(4, 190)} ≥ 5.08, p < 0.001], indicating a general decrease in performance for increasing levels of noise. Additionally, we found a significant main effect of the Algorithm factor on all the tested parameters [F_{(1, 190)} ≥ 11.92, p < 0.001]. The interaction between the two factors was also significant for the similarity of the primitives [F_{(4, 190)} ≥ 5.08, p < 0.001]. The post-hoc analysis revealed that the FADA and fastICA algorithms were equally able to retrieve the correct primitives and weights (p > 0.05), and that FADA-T was better able to reconstruct the simulated signals in the absence of noise (p < 0.05). However, fastICA tended to have higher reconstruction accuracy in the presence of noise (p < 0.05).

FIGURE 2

Figure 2. Identification of the spatial decomposition model. Performance (mean ± SD) of the Fourier-based Anechoic Demixing Algorithm Toolbox (FADA-T), fast Independent Component Analysis (factICA), and Non-negative Matrix Factorization (NMF) on artificial datasets generated with the spatial decomposition model and corrupted by increasing amounts of noise. (A) Unconstrained variant. From left to right: fraction of variance explained by the identified model, normalized similarities between original and identified primitives and activation weights. (B) Corresponding statistics for the model with positivity constraints. The * symbol indicates the statistically significant values of p < 0.05.

Figure 2B reports the identification performance of FADA-T on the spatial decomposition model (1) with non-negativity constraints. In this case, FADA-T was compared against the NMF algorithm (Lee and Seung, 1999), as fastICA does not provide a way to constrain parameters to be non-negative. Also in this case, both algorithms provided a good fit of the data and accurate estimates of the original primitives and mixture weights. Not surprisingly, performance of both algorithms degraded with increasing noise. ANOVAs indicated a significant main effect of the factor Algorithm on the similarity of the weighting coefficients [F_{(1, 190)} = 23.14, p < 0.001]. We also found a main effect of Noise Level on the R² and S_N measures [F_{(4, 190)} ≥ 20.85, p < 0.001]. The interaction of both factors was significant only for the R² measure [F_{(4, 190)} = 5.51,p < 0.001]. The post-hoc analysis showed that FADA-T had a higher reconstruction accuracy (p = 0.01) for the highest tested noise level (35%), and that FADA-T had a lower identification performance (p < 0.05) on the weight coefficients for the 25% noise level. Finally, all the measures were significantly above chance level [t₍₁₉₎ ≥ 6.85, p < 0.001].

Figure 3A shows the identification performance of FADA-T and fastICA on the identification of the parameters of the temporal decomposition model without delays (2). Qualitatively, reconstruction accuracy and primitive similarity were modulated by the noise level, while weight similarity was not. ANOVAs confirmed this trend revealing a significant main effect of Noise on reconstruction accuracy and primitive similarity [F_{(4, 190)} ≥ 9.67, p<0.001]. The interaction between Algorithm and Noise Level was significant only for R² [F_{(4, 190)} = 5.12, p < 0.001]. Post-hoc testing revealed that FADA-T and fastICA performed similarly (p>0.05), with the exception of the reconstruction quality for 35% noise level, where FADA-T outperformed fastICA.

FIGURE 3

Figure 3. Identification of the temporal decomposition model (2). Performance of the FADA-T, factICA, and NMF algorithms on artificial datasets generated with the temporal decomposition model. (A) Unconstrained variant; from left to right: fraction of variance explained by the identified model, normalized similarities between original and identified primitives and activation weights. (B) Corresponding statistics for the model with positivity constraints. The * symbol indicates the statistically significant values of p < 0.05.

Figure 3B shows the results of the comparison between FADA-T and NMF on the identification of the parameters of the temporal decomposition model without delays (2) with non-negativity constraints. The differences in performance between the two methods for the same noise levels were very small. Correspondingly, ANOVAs showed that the Algorithm factor had a significant main effect only on the reconstruction accuracy R² [F_{(1, 190)} = 25.99, p < 0.001], while the Noise factor had significant main effects on all three tested measures [F_{(4, 190)} ≥ 17.38p < 0.001]. Post-hoc testing revealed that, for the two highest noise levels, the NMF algorithm approximated the original data with significantly higher reconstruction accuracy (p < 0.05). Finally, all measures in Figure 3 were significantly above chance level [t₍₁₉₎ ≥ 8.86, p < 0.001].

Taken together, Figures 2, 3 show that, when applied to data generated with the synchronous models (1) and (2), FADA-T exhibited reconstruction performance overall comparable to those provided by the fastICA and NMF algorithms. In terms of the identification of the actual parameters, the differences between the tested algorithms were even smaller, with FADA-T underperforming only on the estimation of the weights of the constrained spatial decomposition model with 25% noise level.

In Figure 4, we show the performance of FADA-T on the identification of the temporal decomposition model with delays (3). For the unconstrained variant (Figure 4A), we compared FADA-T against the AnDem (Omlor and Giese, 2011) and the SICA (Mørup et al., 2007b) algorithms, while for the constrained variant (Figure 4B) we considered the AnNMF (Omlor and Giese, 2011) and the sNMF (Mørup et al., 2007a) algorithms. In addition to the similarity measures assessed for the models considered above, we also quantified the similarity between original and identified delays.

FIGURE 4

Figure 4. Identification of the temporal decomposition model with delays. Performance of FADA-T, Anechoic Demixing (AnDem), Shifted Independent Component Analysis (SICA), Anechoic Demixing with Non-negativity constraints (AnNMF), and shifted Non-negative Matrix Factorization (sNMF), on artificial datasets generated with the temporal decomposition model with delays. (A) Unconstrained variant. From top to bottom: fraction of variance explained by the identified model, normalized similarities between original and identified primitives, activation weights, and delays. (B) Corresponding statistics for the model with positivity constraints. The ** symbol represent two significant pairwise group differences (p < 0.05), each represented with a single asterisk (*).

Figure 4A shows that, overall, FADA-T reconstructed the signals and identified the parameters of the unconstrained anechoic mixture better than the benchmark algorithms, for all the noise levels. ANOVAs revealed a significant main effect of the Algorithm factor on all the considered measures [F_{(2, 190)} ≥ 210.5, p < 0.001]. Post-hoc analysis revealed that, compared to AnDem, FADA-T provided significantly higher reconstruction accuracy, primitive similarity, and weight similarity (p < 0.05). Additionally, compared to SICA, FADA-T had a higher reconstruction accuracy for all noise levels (p < 0.001), and a higher delay similarity for noise levels >15% (p < 0.001). All measures in Figure 4A were significantly different from chance [t₍₁₉₎ ≥ 3.23, p < 0.001].

Figure 4B shows that FADA-T tended to reconstruct the signals and identify the primitives of the non-negatives anechoic mixture model better than the benchmark algorithms, for all the noise levels. Specifically, we found a significant main effect of the factor Algorithm on R², primitive similarity, and delay similarity [F_{(2, 190)} ≥ 6.64, p < 0.05]. Post-hoc testing revealed that, compared to the AnNMF, FADA-T had higher reconstruction accuracy across all noise levels (p < 0.001), a higher primitive similarity for the 5, 15, 25, and 35% noise levels (p < 0.05), and a higher delay similarity for the 15% noise level (p < 0.05). On the other hand, FADA-T had a lower delay similarity than AnNMF for the noise levels 0% and 5%. Compared to sNMF, FADA-T had higher reconstruction accuracy across all noise levels (p < 0.05), and higher primitive similarity for the 0, 5, 25, and 35% noise levels (p < 0.05). All similarity measures in Figure 4B were significantly above chance level [t₍₁₉₎ ≥ 30.8, p < 0.01], except for the reconstruction accuracy provided by sNMF for the most noisy data sets [t₍₁₉₎ = 0.25, p = 0.80].

In Figure 5, we report the results on the identification of the spatiotemporal decomposition model (4). In this case, we compared FADA-T against the stNMF algorithm (d'Avella et al., 2003), developed specifically to fit this model. In this case, we found a significant main effect of the Algorithm factor on the reconstruction accuracy, weight similarity, and delay similarity [F_{(1, 190)} ≥ 13.34, p < 0.001], but not on the primitive similarity [F_{(1, 190)} = 0.4, p > 0.05]. We also found a significant Algorithm by Noise interaction on reconstruction quality and delay similarity [F_{(4, 190)} ≥ 2.84, p < 0.05]. Post-hoc testing revealed that FADA-T had a higher reconstruction accuracy than stNMF across all noise levels (p < 0.001); however, FADA-T had a lower delay similarity for the 35% noise level (p = 0.03). Finally, all the tested measures were significantly above chance level for all the noise levels [t₍₁₉₎ ≥ 11.78, p < 0.001].

FIGURE 5

Figure 5. Identification of the spatiotemporal decomposition model (4). Performance of FADA-T and spatiotemporal NMF (stNMF) on artificial datasets generated with the spatiotemporal decomposition model. The (top) panels report the fraction of explained variance (left) and the normalized similarities between original and identified spatiotemporal primitives (right). The (bottom) panels report the normalized similarities between the corresponding activation weights (left) and delays (right). The * symbol indicates the statistically significant values of p < 0.05.

Figure 6 summarizes the results on the identification of the space-by-time decomposition model (5). In this case, we compared FADA-T against the sNM3F algorithm (Delis et al., 2014), which was developed to fit this specific model. Qualitatively, FADA-T appears to perform better than sNM3F, especially in terms of reconstruction quality, weight similarity and delay similarity. The statistical analyses revealed: a main effect of the Algorithm factor on all the variables [F_{(1, 190)} ≥ 10.72, p < 0.001]; a main effect of Noise on reconstruction accuracy, weight similarity, and delay similarity [F_{(4, 190)} ≥ 2.74, p < 0.05]; a significant Algorithm by Noise interaction on reconstruction accuracy, spatial primitive similarity, and weight similarity [F_{(4, 190)} ≥ 2.46, p < 0.05]. Post-hoc testing showed that, in the presence of noise, FADA-T had a significantly better reconstruction accuracy than sNM3F (p < 0.001). Additionally, FADA-T had a higher temporal primitive similarity for the noise levels 0 and 35% (p < 0.01), and a comparable spatial primitive similarity (p > 0.05). Moreover, FADA-T outperformed sNM3F on the estimation of weights and delays, across all noise levels (p < 0.05). Finally, t-tests showed that FADA-T was always able to provide above chance estimates [t₍₁₉₎ ≥ 3.68, p < 0.01], while sNM3F provided chance-level weight estimates for the 5% noise level [t₍₁₉₎ = 1.91, p = 0.07].

FIGURE 6

Figure 6. Identification of the space-by-time decomposition model. Performance of FADA-T and the sample-based Non-negative Matrix tri-Factorization (sNM3F) algorithm on artificial datasets generated with the space-by-time decomposition model. The (top) panels report the fraction of explained variance (left) and the normalized similarities between ground-truth and estimated temporal (center) and spatial (right) primitives. The (bottom) panels report the normalized similarities between the corresponding activation weights (left) and delays (right). The * symbol indicates the statistically significant values of p < 0.05.

In summary, these results seem to indicate that, on simulated data, FADA-T performs comparably well to model-specific methods on the identification of synchronous mixture models (Figures 2, 3) and better on the identification of anechoic mixtures (Figures 4–6).

3.2. Evaluation of algorithm performance on real experimental data

In addition to validating FADA-T on synthesized data, we also assessed its performance on previously published real experimental data, comparing the primitives extracted by FADA-T with those identified with other techniques. The first experimental data set consisted of kinematic joint angle trajectories of the body joints of participants performing emotional walks. Trajectories represented a single gait cycle, resampled with 100 time steps (Roether et al., 2009; Endres et al., 2013). In this case, we tested FADA-T against the AnDem (Omlor and Giese, 2011) and the SICA (Mørup et al., 2007b) algorithms on the extraction of delayed temporal primitives (3). Figure 7A shows the fraction of explained variance (R²), as a function of the number of primitives. Consistently with previous studies (i.e., d'Avella et al., 2006; Berger et al., 2020), to select the number of primitives that allow a good trade-off between reconstruction accuracy and model complexity, we used the elbow method: we selected the minimum number of primitives from which a line could fit the R² curve well (i.e., with a mean squared error below 10⁻⁴—d'Avella et al., 2006; Berger et al., 2020). For all methods, this point was reached for N = 3, indicating that three anechoic components (Figure 7B) can approximate the experimental data set well. Interestingly, the primitives identified by the methods were almost identical (S ≥ 0.94).

FIGURE 7

Figure 7. Extraction of delayed temporal primitives from kinematic data of emotional walking. (A) Fraction of explained variance as a function of the number of extracted primitives, identified with FADA-T (cyan), AnDem (green), and SICA (purple). (B) Temporal primitives identified by the three algorithms for a model with three sources; legends indicate pairwise correlation values.

The second experimental dataset we used to validate FADA-T comprised EMG signals recorded from 16 different task-relevalant muscles during point-to-point arm reaching movements (d'Avella et al., 2006). In this case, we fitted the spatio-temporal decomposition model (4)—which was originally shown to capture important features of the dataset (d'Avella et al., 2006)—with both FADA-T and the stNMF algorithm. Figure 8 shows that the methods displayed very similar R² curves (first column), which leveled off at N = 5 (as assessed by the elbow method). Additionally, also in this case, the identified spatiotemporal primitives identified by the two methods were very similar (S ≥ 0.85).

FIGURE 8

Figure 8. Extraction of spatiotemporal primitives from EMG data during reaching movements. Primitives are extracted with the FADA-T (A) and stNMF (B), and are grouped according to their similarity. The leftmost panels show the fraction of explained variance as function of the number of extracted primitives. The remaining panels show the activation patterns of the spatiotemporal primitives extracted by the two algorithms.

Taken together, these results suggest that FADA-T performs as well as model-specific identification methods on both real kinematic (Figure 7) and muscle activity data (Figure 8).

4. Discussion

In this work, we have introduced a novel framework that allows to unify a number of common methods for the definition and identification of spatial, temporal, and spatiotemporal motor primitives. The framework harnesses the flexibility of the anechoic mixture model to capture qualitatively different classes of motor modularity models, and the robustness of the Fourier-based Anechoic Demixing Algorithm (FADA) to estimate the parameters reliably. We tested the framework on eight different model classes on both simulated and experimental data and showed that the reconstruction and identification performance were in most cases comparable to those achieved by established model-specific identification methods. As the experimental data comprised both electromyographic and kinematic spatio-temporal signals collected during the execution of qualitatively different motor behaviors, our results suggest that the framework constitutes a valid unsupervised method to identify spatial, temporal, and spatiotemporal regularities from signals extracted from different levels of the motor hierarchy. Our framework has thus the potential to facilitate the identification of motor modules, and to remove the potential confounding arising when comparing results obtained adopting different models and identification algorithms.

The algorithmic framework we used to identify the sources, delays, and mixing weights of the anechoic mixture model builds on the FADA algorithm, which we originally introduced in Chiovetto and Giese (2013) for the estimation of the temporal decomposition model with delays (3). In the present work, we have extended FADA to make it suitable for the identification of other popular models of spatial and spatiotemporal decomposition. This was accomplished by deriving classes of relevant constraints and by appropriately adapting individual optimization steps. The resulting classes of algorithms, which we have released publicly in the FADA toolbox (FADA-T), solve the well-known problem of over-determined anechoic demixing, where the number of signals to reconstruct outnumbers that of the latent source functions. A few algorithms have been proposes to address such a problem, including the Shifted Factor Analysis (Harshman et al., 2003), the Shifted Independent Component Analysis (Mørup et al., 2007b), and the Anechoic Demixing Algorithm (Omlor and Giese, 2011). However, these algorithms tend to be computationally expensive, as they do not sufficiently restrict the search space of the latent sources. Our identification framework, on the other hand, imposes a smoothness prior on the sources; this restricts the search space to band-limited source function, which speeds up the estimation process.

Influential studies have compared the performance of different existing identification methods for the estimation of motor primitives (e.g., Tresch et al., 2006; Lambert-Shirzad and Van der Loos, 2017; Ebied et al., 2018), or for the blind source decomposition of other biological signals (Virtanen et al., 2009; Cashero and Anderson, 2011; Erhardt et al., 2011). In contrast to these studies, here we have introduced a novel framework for the identification of motor primitives; moreover, we have shown that the framework allows the estimation of several classes of motor primitives and have benchmarked it against popular model-specific identification methods. Similarly to these studies, we have used both simulated and real experimental data to measure the identification performance.

A key element of the FADA algorithm is the mapping onto a finite Fourier basis. On the one hand, this strategy reduces remarkably the number of identified parameters in comparison to more general anechoic demixing methods (Mørup et al., 2007b; Omlor and Giese, 2011). On the other hand, this choice determines that only band-limited data can be adequately modeled. However, this is usually not an issue of major concern in fields such as motor control, where the typical activity patterns of interest (e.g., joint kinematics, EMGs, etc.) are naturally band-limited or artificially smoothed by filtering and trial averaging. An additional potential limitation of our work is that the current implementation of FADA-T includes some constrained optimization steps performed with gradient descent methods, which can potentially be time-consuming. Even though FADA-T performs such steps relying on built-in MATLAB functions, future work can potentially replace them with more efficient ones, harnessing the highly modular pipeline of FADA-T.

Empirically, we found that the imposed reduction of the dimensionality of the parameter space resulted in a more robust estimation of the primitives (even in the presence of substantial noise levels) and in a lower probability of convergence to local minima. Consequently, FADA-T performed better than other methods on the identification of anechoic mixture models (Figures 4–6). This seems to be due to a stronger ability to identify the correct weights and delays, especially in the cases of the space-by-time decomposition model (Figure 6), and the temporal decomposition model with delays (Figure 4A). Preliminary analyses, also suggest a better ability of FADA-T to deal with ambiguities in the estimation of delays and source functions, especially for sources with higher fundamental frequencies. Further studies, which are beyond the scope of this paper, will investigate this issue in more depth.

The only case where FADA-T showed consistently lower performance was on the identification of the unconstrained spatial decomposition model (Figure 2A), where fastICA had higher reconstruction accuracy in the presence of noise. We speculate that this happened because the smoothness prior imposed by FADA-T on the spatial domain was perhaps too restrictive. However, in spite of this problem, both the reconstruction and identification accuracy were sufficiently high. In addition to testing FADA-T on simulated datasets, we also considered the problem of identifying spatiotemporal sources from real experimental data of EMG and kinematic data sets, collected from human participants during the execution of goal-oriented and rhythmic motor tasks. In this case, our results show that FADA-T retrieved sources and mixing coefficients that were consistent with those obtained with other traditional techniques (cf. Figures 7, 8).

The central mathematical contribution of this article is the systematical analysis of the relationship between the different models of motor modularity, and, most importantly, their common derivation from the anechoic mixture model (6). This raises the question of how to identify the most appropriate modularity model for a given experimental dataset. Even though we have not proposed new solutions to this problem in this work, classical model selection methods can easily be applied to this context, including the Akaike method (Akaike, 1974) and the Bayesian Information Criterion (Schwarz, 1978). Alternative approaches such as Bayesian model selection (Bishop and Nasrabadi, 2006) can also be potentially applied. For this purpose, all tested models are embedded in a joint model space, and one marginalizes the prediction error (evidence) using an uninformative prior distribution over all possible model architectures. This procedure typically finds automatically a good balance between the goodness-of-fit and simplicity of the model. We have proposed an implementation of this idea for automatic model selection in Endres et al. (2013), where we approximated the resulting non-Gaussian distributions with a Laplace approximation. This allowed us to obtain an analytically tractable criterion to compare different demixing models, including those with time delays. Interestingly, a similar procedure can also allow to make inference about the most suitable smoothness priors for a given data set. A more recent AIC-based approach for the estimation of the number of motor primitives extracted with NMF is introduced by Ranaldi et al. (2021).

5. Conclusion

Experimental investigations over the last couple of decades have confirmed Bernstein's hypothesis (Bernstein, 1966) that the CNS simplifies the control of movement by relying on a modular organization of control (Flash and Hochner, 2005; Bizzi et al., 2008; Bizzi and Ajemian, 2020). The modules underlying such a control architecture have been defined in multiple ways (Tresch et al., 2006; Chiovetto et al., 2013; Giszter, 2015), and extracted by applying a variety of unsupervised learning algorithms to kinematic, kinetic, EMG, and neural data. In this work, we have introduced a unifying framework as a potential solution to this heterogeneity of approaches. FADA-T—the toolbox we have developed to provide a single estimation environment for the most common models of motor modularity—promises to facilitate the interpretation of the multimodal data recorded during the execution of body movements, by simplifying the process of identifying the most appropriate modularity models for the dataset at hand. We expect FADA-T to also stimulate the explorative adoption of the motor modularity models to other neuroscientific domains, which can potentially lead to the discovery of similar principles of hierarchical organization in other functional brain systems.

Data availability statement

The datasets analyzed in this study are available from the corresponding author on reasonable request.

Ethics statement

Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent from the patients/participants or patients/participants' legal guardian/next of kin was not required to participate in this study in accordance with the national legislation and the institutional requirements.

Author contributions

EC, AS, Ad'A, and MG have made a significant contribution to the design of the study, the development and implementation of the methods, the analysis of the results, and the writing of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the German Federal Ministry of Education and Research (BMBF FKZ 01GQ1704), the Human Frontiers Science Program (HFSP RGP0036/2016), the German Research Foundation (DFG GZ: KA 1258/15- 1), the European Research Council (ERC 2019-SyGRELEVANCE- 856495), and the European Union Horizon 2020 Programme (CogIMon H2020 ICT-23-2014/644727).

Acknowledgments

The authors thank Lars Omlor for his help during the first stage of this project and Dominik Endres for many helpful discussions. The authors also thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting AS.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2022.926345/full#supplementary-material

References

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Autom. Control 19, 716–723. doi: 10.1109/TAC.1974.1100705

CrossRef Full Text | Google Scholar

Berger, D. J., Masciullo, M., Molinari, M., Lacquaniti, F., and d'Avella, A. (2020). Does the cerebellum shape the spatiotemporal organization of muscle patterns? Insights from subjects with cerebellar ataxias. J. Neurophysiol. 123, 1691–1710. doi: 10.1152/jn.00657.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernstein, N. (1966). The co-ordination and regulation of movements. Pergamon Press.

Google Scholar

Bernstein, N. A. (1947). On the Construction of Movements. Moscow: Medgiz.

Berret, B., Bonnetblanc, F., Papaxanthis, C., and Pozzo, T. (2009). Modular control of pointing beyond arm's length. J. Neurosci. 29, 191–205. doi: 10.1523/JNEUROSCI.3426-08.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bishop, C. M., and Nasrabadi, N. M. (2006). Pattern Recognition and Machine Learning. Vol. 4. New York, NY: Springer.

Google Scholar

Bizzi, E., and Ajemian, R. (2020). From motor planning to execution: a sensorimotor loop perspective. J. Neurophysiol. 124, 1815–1823. doi: 10.1152/jn.00715.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Bizzi, E., Cheung, V. C., d'Avella, A., Saltiel, P., and Tresch, M. (2008). Combining modules for movement. Brain Res. Rev. 57, 125–133. doi: 10.1016/j.brainresrev.2007.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bofill, P. (2003). Underdetermined blind separation of delayed sound sources in the frequency domain. Neurocomputing 55, 627–641. doi: 10.1016/S0925-2312(02)00631-8

CrossRef Full Text | Google Scholar

Cashero, Z., and Anderson, C. W. (2011). “Comparison of EEG blind source separation techniques to improve the classification of p300 trials,” in 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Boston, MA: IEEE), 7183–7186. doi: 10.1109/IEMBS.2011.6091815

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheung, V. C., and Seki, K. (2021). Approaches to revealing the neural basis of muscle synergies: a review and a critique. J. Neurophysiol. 125, 1580–1597. doi: 10.1152/jn.00625.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiovetto, E., Berret, B., Delis, I., Panzeri, S., and Pozzo, T. (2013). Investigating reduction of dimensionality during single-joint elbow movements: a case study on muscle synergies. Front. Comput. Neurosci. 7:11. doi: 10.3389/fncom.2013.00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiovetto, E., Berret, B., and Pozzo, T. (2010). Tri-dimensional and triphasic muscle organization of whole-body pointing movements. Neuroscience 170, 1223–1238. doi: 10.1016/j.neuroscience.2010.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiovetto, E., and Giese, M. A. (2013). Kinematics of the coordination of pointing during locomotion. PLoS ONE 8:e79555. doi: 10.1371/journal.pone.0079555

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiovetto, E., Patané, L., and Pozzo, T. (2012). Variant and invariant features characterizing natural and reverse whole-body pointing movements. Exp. Brain Res. 218, 419–431. doi: 10.1007/s00221-012-3030-y

PubMed Abstract | CrossRef Full Text | Google Scholar

d'Avella, A., Giese, M., Ivanenko, Y. P., Schack, T., and Flash, T. (2015). Modularity in motor control: from muscle synergies to cognitive action representation. Front. Comput. Neurosci. 9:126. doi: 10.3389/fncom.2015.00126

PubMed Abstract | CrossRef Full Text | Google Scholar

d'Avella, A., Portone, A., Fernandez, L., and Lacquaniti, F. (2006). Control of fast-reaching movements by muscle synergy combinations. J. Neurosci. 26, 7791–7810. doi: 10.1523/JNEUROSCI.0830-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

d'Avella, A., Saltiel, P., and Bizzi, E. (2003). Combinations of muscle synergies in the construction of a natural motor behavior. Nat. Neurosci. 6, 300–308. doi: 10.1038/nn1010

PubMed Abstract | CrossRef Full Text | Google Scholar

d'Avella, A., and Tresch, M. C. (2002). “Modularity in the motor system: decomposition of muscle patterns as combinations of time-varying synergies,” in Advances in Neural Information Processing Systems, Vol. 1 (Vancouver, BC), 141–148.

Google Scholar

Delis, I., Panzeri, S., Pozzo, T., and Berret, B. (2014). A unifying model of concurrent spatial and temporal modularity in muscle activity. J. Neurophysiol. 111, 675–693. doi: 10.1152/jn.00245.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominici, N., Ivanenko, Y. P., Cappellini, G., d'Avella, A., Mondì, V., Cicchese, M., et al. (2011). Locomotor primitives in newborn babies and their development. Science 334, 997–999. doi: 10.1126/science.1210617

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebied, A., Kinney-Lang, E., Spyrou, L., and Escudero, J. (2018). Evaluation of matrix factorisation approaches for muscle synergy extraction. Med. Eng. Phys. 57, 51–60. doi: 10.1016/j.medengphy.2018.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Emile, B., and Comon, P. (1998). Estimation of time delays between unknown colored signals. Signal Process. 69, 93–100. doi: 10.1016/S0165-1684(98)00061-9

CrossRef Full Text | Google Scholar

Endres, D., Chiovetto, E., and Giese, M. A. (2013). Model selection for the extraction of movement primitives. Front. Comput. Neurosci. 7:185. doi: 10.3389/fncom.2013.00185

PubMed Abstract | CrossRef Full Text | Google Scholar

Erhardt, E. B., Rachakonda, S., Bedrick, E. J., Allen, E. A., Adali, T., and Calhoun, V. D. (2011). Comparison of multi-subject ICA methods for analysis of fMRI data. Hum. Brain Mapp. 32, 2075–2095. doi: 10.1002/hbm.21170

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Jiang, N., Rehbaum, H., Holobar, A., Graimann, B., Dietl, H., and Aszmann, O. C. (2014). The extraction of neural information from the surface EMG for the control of upper-limb prostheses: emerging avenues and challenges. IEEE Trans. Neural Syst. Rehabil. Eng. 22, 797–809. doi: 10.1109/TNSRE.2014.2305111

PubMed Abstract | CrossRef Full Text | Google Scholar

Flash, T., and Hochner, B. (2005). Motor primitives in vertebrates and invertebrates. Curr. Opin. Neurobiol. 15, 660–666. doi: 10.1016/j.conb.2005.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Giszter, S. F. (2015). Motor primitives-new data and future questions. Curr. Opin. Neurobiol. 33, 156–165. doi: 10.1016/j.conb.2015.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Godlove, J., Gulati, T., Dichter, B. K., Chang, E., and Ganguly, K. (2016). Muscle synergies after stroke are correlated with perilesional high gamma. Ann. Clin. Transl. Neurol. 3, 956–961. doi: 10.1002/acn3.368

PubMed Abstract | CrossRef Full Text | Google Scholar

Harshman, R. A., Hong, S., and Lundy, M. E. (2003). Shifted factor analysis-part I: models and properties. J. Chemometr. 17, 363–378. doi: 10.1002/cem.808

CrossRef Full Text | Google Scholar

Hart, C. B., and Giszter, S. F. (2004). Modular premotor drives and unit bursts as primitives for frog motor behaviors. J. Neurosci. 24, 5269–5282. doi: 10.1523/JNEUROSCI.5626-03.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hart, C. B., and Giszter, S. F. (2010). A neural basis for motor primitives in the spinal cord. J. Neurosci. 30, 1322–1336. doi: 10.1523/JNEUROSCI.5894-08.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Højen-Sørensen, P. A. F. R., Winther, O., and Hansen, L. K. (2002). Mean-field approaches to independent component analysis. Neural Comput. 14, 889–918. doi: 10.1162/089976602317319009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hyvärinen, A., and Oja, E. (1997). A fast fixed-point algorithm for independent component analysis. Neural Comput. 9, 1483–1492. doi: 10.1162/neco.1997.9.7.1483

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Cappellini, G., Dominici, N., Poppele, R. E., and Lacquaniti, F. (2005). Coordination of locomotion with voluntary movements in humans. J. Neurosci. 25, 7238–7253. doi: 10.1523/JNEUROSCI.1327-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2004). Five basic muscle activation patterns account for muscle activity during human locomotion. J. Physiol. 556, 267–282. doi: 10.1113/jphysiol.2003.057174

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadmon Harpaz, N., Ungarish, D., Hatsopoulos, N. G., and Flash, T. (2019). Movement decomposition in the primary motor cortex. Cereb. Cortex 29, 1619–1633. doi: 10.1093/cercor/bhy060

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaminski, T. (2007). The coupling between upper and lower extremity synergies during whole body reaching. Gait Posture 26, 256–262. doi: 10.1016/j.gaitpost.2006.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert-Shirzad, N., and Van der Loos, H. M. (2017). On identifying kinematic and muscle synergies: a comparison of matrix factorization methods using experimental data from the healthy population. J. Neurophysiol. 117, 290–302. doi: 10.1152/jn.00435.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, D., and Seung, H. S. (2000). “Algorithms for non-negative matrix factorization,” in Advances in Neural Information Processing Systems, Vol. 13, eds T. Leen, T. Dietterich, and V. Tresp (Denver, CO: MIT Press).

PubMed Abstract | Google Scholar

Lee, D. D., and Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature 401, 788–791. doi: 10.1038/44565

PubMed Abstract | CrossRef Full Text | Google Scholar

Levine, A. J., Hinckley, C. A., Hilde, K. L., Driscoll, S. P., Poon, T. H., Montgomery, J. M., et al. (2014). Identification of a cellular node for motor control pathways. Nat. Neurosci. 17, 586–593. doi: 10.1038/nn.3675

PubMed Abstract | CrossRef Full Text | Google Scholar

Merel, J., Botvinick, M., and Wayne, G. (2019). Hierarchical motor control in mammals and machines. Nat. Commun. 10, 1–12. doi: 10.1038/s41467-019-13239-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Mori, H., and Hoshino, J. (2002). “Independent component analysis and synthesis of human motion,” in 2002 IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 4 (Orlando, FL), IV-3564. doi: 10.1109/ICASSP.2002.5745425

CrossRef Full Text | Google Scholar

Mørup, M., Madsen, K., and Hansen, L. (2007a). “Shifted non-negative matrix factorization,” in 2007 IEEE Workshop on Machine Learning for Signal Processing (Thessaloniki: IEEE), 139–144. doi: 10.1109/MLSP.2007.4414296

CrossRef Full Text | Google Scholar

Mørup, M., Madsen, K. H., and Hansen, L. K. (2007b). “Shifted independent component analysis,” in International Conference on Independent Component Analysis and Signal Separation (London: Springer), 89–96. doi: 10.1007/978-3-540-74494-8_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Mussa-Ivaldi, F. A., and Giszter, S. F. (1992). Vector field approximation: a computational paradigm for motor control and learning. Biol. Cybern. 67, 491–500. doi: 10.1007/BF00198756

PubMed Abstract | CrossRef Full Text | Google Scholar

Omlor, L., and Giese, M. (2006). “Blind source separation for over-determined delayed mixtures,” in Advances in Neural Information Processing Systems, Vol. 19 (Vancouver, BC), 1049–1056.

Google Scholar

Omlor, L., and Giese, M. A. (2007). Extraction of spatio-temporal primitives of emotional body expressions. Neurocomputing 70, 1938–1942. doi: 10.1016/j.neucom.2006.10.100

CrossRef Full Text | Google Scholar

Omlor, L., and Giese, M. A. (2011). Anechoic blind source separation using Wigner marginals. J. Mach. Learn. Res. 12, 1111–1148. doi: 10.5555/1953048.2021037

CrossRef Full Text | Google Scholar

Overduin, S. A., d'Avella, A., Roh, J., Carmena, J. M., and Bizzi, E. (2015). Representation of muscle synergies in the primate brain. J. Neurosci. 35, 12615–12624. doi: 10.1523/JNEUROSCI.4302-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranaldi, S., De Marchis, C., Severini, G., and Conforto, S. (2021). An objective, information-based approach for selecting the number of muscle synergies to be extracted via non-negative matrix factorization. IEEE Trans. Neural Syst. Rehabil. Eng. 29, 2676–2683. doi: 10.1109/TNSRE.2021.3134763

PubMed Abstract | CrossRef Full Text | Google Scholar

Roether, C. L., Omlor, L., Christensen, A., and Giese, M. A. (2009). Critical features for the perception of emotion from gait. J. Vis. 9:15. doi: 10.1167/9.6.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Santello, M., Flanders, M., and Soechting, J. F. (1998). Postural hand synergies for tool use. J. Neurosci. 18, 10105–10115. doi: 10.1523/JNEUROSCI.18-23-10105.1998

CrossRef Full Text | Google Scholar

Santello, M., and Soechting, J. F. (2000). Force synergies for multifingered grasping. Exp. Brain Res. 133, 457–467. doi: 10.1007/s002210000420

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwarz, G. (1978). Estimating the dimension of a model. Ann. Stat. 6, 461–464. doi: 10.1214/aos/1176344136

CrossRef Full Text | Google Scholar

Smith, C., Gilleard, W. L., Hammond, J., and Brooks, L. O. (2006). The application of an exploratory factor analysis to investigate the inter-relationships amongst joint movement during performance of a football skill. J. Sports Sci. Med. 5, 417–524.

PubMed Abstract | Google Scholar

Steinberg, F., and Bock, O. (2013). Influence of cognitive functions and behavioral context on grasping kinematics. Exp. Brain Res. 225, 387–397. doi: 10.1007/s00221-012-3379-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Swindlehurst, A. L. (1998). Time delay and spatial signature estimation using known asynchronous signals. IEEE Trans. Signal Process. 46, 449–462. doi: 10.1109/78.655429

CrossRef Full Text | Google Scholar

Takei, T., Confais, J., Tomatsu, S., Oya, T., and Seki, K. (2017). Neural basis for hand muscle synergies in the primate spinal cord. Proc. Natl. Acad. Sci. U.S.A. 114, 8643–8648. doi: 10.1073/pnas.1704328114

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, J. S., Corcos, D. M., and Hasan, Z. (2005). Kinematic and kinetic constraints on arm, trunk, and leg segments in target-reaching movements. J. Neurophysiol. 93, 352–364. doi: 10.1152/jn.00582.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ting, L. H., and Macpherson, J. M. (2005). A limited set of muscle synergies for force control during a postural task. J. Neurophysiol. 93, 609–613. doi: 10.1152/jn.00681.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Torkkola, K. (1996). “Blind separation of convolved sources based on information maximization,” in Neural Networks for Signal Processing VI. Proceedings of the 1996 IEEE Signal Processing Society Workshop (Kyoto: IEEE), 423–432. doi: 10.1109/NNSP.1996.548372

PubMed Abstract | CrossRef Full Text | Google Scholar

Torres-Oviedo, G., Macpherson, J. M., and Ting, L. H. (2006). Muscle synergy organization is robust across a variety of postural perturbations. J. Neurophysiol. 96, 1530–1546. doi: 10.1152/jn.00810.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Tresch, M. C., Cheung, V. C. K., and d'Avella, A. (2006). Matrix factorization algorithms for the identification of muscle synergies: evaluation on simulated and experimental data sets. J. Neurophysiol. 95, 2199–2212. doi: 10.1152/jn.00222.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Tresch, M. C., Saltiel, P., and Bizzi, E. (1999). The construction of movement by the spinal cord. Nat. Neurosci. 2, 162–167. doi: 10.1038/5721

PubMed Abstract | CrossRef Full Text | Google Scholar

Virtanen, J., Noponen, T. E., and Meriläinen, P. (2009). Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals. J. Biomed. Opt. 14:054032. doi: 10.1117/1.3253323

PubMed Abstract | CrossRef Full Text | Google Scholar

Yilmaz, O., and Rickard, S. (2004). Blind separation of speech mixtures via time-frequency masking. IEEE Trans. Signal Process. 52, 1830–1847. doi: 10.1109/TSP.2004.828896

CrossRef Full Text | Google Scholar