The trial-to-trial variability within the neural commands differed depending on both the MSA method and the method of quantifying variability (**Table 3**). For both cosine similarity and cross-correlation, negative effect sizes indicate greater trialto-trial variability among RES, compared to nRES. For the cosine similarity metric with two synergies and constant SVs, there were three medium effects and no large or very large effects (**Figure 5**). Two synergies with varying SVs produced large effects for all fifteen combinations of EMG and output normalization. Three synergies with constant SVs revealed five medium effects and two large effects. Three synergies with varying SVs revealed five medium and no large or very large effects. For the cross-correlation comparison with constant SVs, neither two nor three synergies produced notable effects (**Figure 6**). Two synergies with varying SVs revealed six medium effects and six large effects (**Table 3**). Three synergies with varying SVs revealed three medium effects and one large effect. In general, the same patterns

#### TABLE 2 | Notable effects in synergy vector variability.

SVs, synergy vectors; MagPer, magnitude per trial; SV Mag, synergy vector magnitude; SV Max, synergy vector maximum value; NC Max, neural command maximum value; MaxPer, maximum value per trial; UnitPer, unit variance per trial; MaxOver, maximum value over all trials; UnitOver, unit variance over all trials.

aEffect sizes calculated using Hedge's g. Positive values in SV variability indicate greater variability in RES.

FIGURE 3 | Absolute effect sizes at baseline, measured by variance accounted for (VAF). Effect sizes for two synergies with constant synergy vectors (SVs, A), two synergies with varying SVs (B), three synergies with constant SVs (C), and three synergies with varying SVs (D). EMG normalization methods include: magnitude per trial (MagPer), maximum value over all trials (MaxOver), maximum value per trial (MaxPer), unit variance per trial (UnitPer), and unit variance over all trials (UnitOver). Color bar represents effect size values.

were present across output normalization methods but two synergies with varying SVs produced large effects most frequently.

represents effect size values.

To validate our findings, we compared results from the original nRES to those from five additional nRES identified from the larger data set. Example synergies for a validation nRES are shown in **Figure 2** (bottom row), which appear qualitatively similar to the example nRES in the same figure. Two synergies were sufficient to achieve 90% VAF in all validation nRES across all methods, and there were no differences in VAF between the original nRES and the validation group. Because two synergies with varying SVs produced the most consistent effects in the initial RES/nRES analysis, validation focused on large effects within this subset of results. In this case, large effects in the initial RES/nRES analysis coupled with small or no effects in the nRES validation analysis are desirable findings, indicating that the method has the capacity to differentiate between response groups. For synergy vector variability, MagPer EMG normalization produced large effects in the initial analysis and no effects when comparing nRES to nRES TABLE 3 | Notable effects in neural command variability.

SVs, synergy vectors; MagPer, magnitude per trial; SV Mag, synergy vector magnitude; SV Max, synergy vector maximum; NC Max, neural command maximum; MaxOver, maximum value over all trials; MaxPer, maximum value per trial; UnitPer, unit variance per trial; UnitOver, unit variance over all trials; –, no notable effect.

aEffect sizes calculated using Hedge's g. Negative values in NC cosine similarity and NC cross-correlation indicate greater variability in RES.

bNeural command similarity compared by trial-to-trial cosine similarity.

cNeural command similarity compared by trial-to-trial cross-correlation.

(**Table 4**). For cosine similarity within the neural commands, MagPer EMG normalization produced medium effects when comparing between nRES and validation nRES across all output normalization methods. All other EMG methods that produced large effects from the initial analysis produced small or no effects in the validation analysis. For NC cross-correlation, all comparisons within the validation analysis revealed small or no effects.

## DISCUSSION

This study demonstrates that the methodological choices made during MSA have a significant impact on the outcome. Several, but not all, MSA methods were able to differentiate between two groups of pathologic individuals with no distinguishing clinical differences, but presumed intrinsic physiologic differences. Two synergies with varying SVs across trials produced the highest frequency and greatest magnitude of effects. When assessing SV variability, MagPer EMG normalization was the most salient method within this analysis. However, when comparing NCs by cosine similarity or cross-correlation, any of the EMG normalization methods other than MagPer produced differences between RES and nRES that were validated in a secondary analysis. MSA output normalization had no notable influence on results. In general, our results were highly sensitive to changes in MSA methodology, illustrating the need for careful methodological consideration and disclosure when conducting MSA or comparing with results found in the literature.

## Methodological Considerations

We do not believe that our choice of one signal decomposition method (i.e., NNMF) limits the ability to detect group differences within our results. We chose to use NNMF for this MSA for three reasons: muscle activation signals are inherently non-negative, other methods require assumptions such as orthogonality and independence of parameters (Ting and Chvatal, 2010), and NNMF is a commonly used technique, allowing for greater generalizability of results. Previous work using MSA illustrates that NNMF produces synergy vectors and neural commands that are highly correlated with those produced by other computational methods such as PCA and ICA (Ivanenko et al., 2005; Cappellini et al., 2006). NNMF is particularly robust to differences in data distribution and noise across data sets (Tresch et al., 2006). NNMF tends to define a solution subspace where the synergies can be found, building non-negative components (i.e., SVs and NCs) together to reconstruct the original signals (Ting and Chvatal, 2010). This approach is appropriate for the physiological interpretation of the muscle synergy data. There are several variations of the NNMF algorithm and it is possible that changing the algorithm could change MSA results. However, the current consensus in the literature is that the decomposition method does not drastically change the structure of the extracted synergies (Ivanenko et al., 2005, 2006; Cappellini et al., 2006; Tresch et al., 2006). Ultimately, the choice of decomposition method should be tailored to the intended application of the results.

It is common for investigators to specify a minimum VAF cutoff of 90% for identification of the number of synergies that can be used to reconstruct the original signals

(Clark et al., 2010; Frère and Hug, 2012; Roh et al., 2012). If this common cutoff was applied to the current analysis, two synergies would have been more than adequate to reconstruct the original EMG walking signals with reasonable accuracy. If a cutoff of 95% VAF for all methods and all subjects was applied, all but one subject would have required three synergies, with the remaining subject requiring four. The decision to use either two or three synergies simplifies the analysis and prevents potential over-fitting problems. This result is also consistent with findings that individuals post-stroke require fewer synergies to reconstruct their movements (Clark et al., 2010; Cheung et al., 2012). However, three varying synergies were also able to detect differences between groups for some methods.

Decisions beyond decomposition method and number of synergies, such as fixing or varying SVs and EMG normalization method, exhibited a profound influence on the ability to detect group differences within our sample. Allowing SVs to vary from trial-to-trial, rather than fixing them across trials, produced more frequent and larger effect sizes, although this choice is less common in the literature (Ivanenko et al., 2005; Oliveira et al., 2014). MagPer EMG normalization produced large differences between groups for SV variability. When comparing NCs with cosine similarity, all methods of output normalization and EMG normalization (except MagPer) produced group differences that were validated with further analyses. When applying a cross-correlation analysis, MaxPer and UnitPer EMG normalization were the best performing methods across all output normalization approaches. EMG normalization is a common step in MSA, however one utilized without much justification for the choice of method. The impact of normalization strategies is shown when one investigator attempts to replicate the analysis of another. Two research groups performed the same experiment with one added condition, but found different results, which could be attributable in part to differences in MSA methods (de Rugy et al., 2013; Berger and d'Avella, 2014). Our analysis indicates that choice of normalization methods and the choice of constant or variable SVs across trials have a strong effect on the outcomes, thus representing a non-trivial choice among MSA parameters.

It is conceivable that variation in the selection of muscles included in the analysis could impact results (Steele et al., 2013). Addition of more muscles or muscles that could more accurately capture the differences between RES and nRES would likely increase the resolution to detect group differences. These data, and other synergy analysis data, should be interpreted within the context of the muscles included in the analysis.

## Detecting Physiological Differences

Differences in synergy vector and neural command variability between RES and nRES were greatest using two synergies with varying SVs. The largest difference between the two groups was found in the SVs when EMG was normalized to its maximum

TABLE 4 | Comparison of initial results for two synergies with varying synergy vectors with non-responder validation analysis.

RES, responders; nRES, non-responders; SV, synergy vector; NC, neural command; MagPer, magnitude per trial; SV Mag, synergy vector magnitude; SV Max, synergy vector maximum; NC Max, neural command maximum; MaxOver, maximum value over all trials; MaxPer, maximum value per trial; UnitPer, unit variance per trial; UnitOver, unit variance over all trials. aEffect sizes calculated using Hedge's g. Positive values in SV variability and negative values in NC cosine similarity and NC cross-correlation indicate greater variability in RES.

value per trial, while EMG normalization was less important when detecting differences in the neural commands. The effect sizes indicate that RES had more variability in their SVs than did nRES. Previous findings indicate that the SVs, but not the NCs, reveal differences between individuals post-stroke and healthy controls (Gizzi et al., 2011). Our results differ slightly, indicating that variability in the SVs and NCs can be used to quantify differences between groups after stroke. However, we compared two groups, both with neural pathology. Our most salient results could not have been produced if the SVs were held constant, an assumption that many investigators have used in their analyses (Clark et al., 2010; Ting and Chvatal, 2010). Differences in trial-to-trial neural command similarity between RES and nRES were also greatest when two varying synergy vectors were used. Similar to our findings in the SVs, these effects indicate that RES exhibited more trial-to-trial variability in the NCs than nRES. This characteristic could be due to a greater flexibility of commands that can be combined to produce a richer variety of movement patterns (Bernstein, 1967), which could be relevant to identifying one's capacity for recovery or treatment response. When NCs were compared using cross-correlation instead of cosine similarity, the results were similar; however, the effects were often blunted. This difference in effect magnitude provides further indication regarding the importance of all aspects of MSA with respect to the final results.

## Recommendations for Reporting Results

MSA has been popularized due to its ability to reduce the dimensionality of complex patterns of muscle activity (Ivanenko et al., 2006). However, the rapid popularization of a technique, coupled with the lack of validation of each MSA methodological choice, presents implementation challenges. In publications utilizing synergy analysis, some investigators provide more than adequate disclosure of the methods employed; however, others provide little to no detail, making comparison of results a challenge. We therefore propose a list of decisions that should be reported explicitly in manuscripts performing MSA:

The present comparison of MSA methods is not intended to suggest there is any one "best" method for all future applications of MSA. Rather, our results illustrate and emphasize the vast differences produced by variations in methodological choices. Ultimately, the choice of analysis methodology should be tailored to the application and research design. To facilitate understanding and reproducibility we recommend disclosure and justification of methods. Such openness between investigators will move research forward, improving the likelihood and timeliness of a research impact involving MSA.

#### CONCLUSIONS

The main goal of MSA, when applied to pathologic populations, is to better understand the intrinsic physiologic characteristics reflected in muscle activity. MSA, with specific focus on trial-to-trial variability, has the potential to provide insight regarding neural strategies that could be relevant to human performance and rehabilitation. Notably, our analysis revealed large differences between response groups with only 10 subjects. Clinical assessments typically require large samples, lack sensitivity, and in this case, fail to differentiate subpopulations that would later respond to an intervention. If information from MSA can be successfully employed as a predictor of recovery or intervention response on an individual basis, this information could improve outcomes of neurorehabilitation. The analysis must, however, be performed properly, with careful selection and justification of methods.

#### REFERENCES

#### AUTHOR CONTRIBUTIONS

BF, CP, and MP conceived and designed experiments; MP, TM, and BF contributed analysis tools; CB and MP performed experiments; CB analyzed data, prepared figures, and drafted manuscript; CB, BF, and CP interpreted results of experiments; CB, BF, and CP revised manuscript; CB, MP, TM, BF, and CP approved final version of manuscript and agreed to be accountable for all aspects of the work.

#### FUNDING

This research was supported by the Department of Veterans Affairs, Rehabilitation Research & Development Service [Project #B29792R (CP) and Research Career Scientist Award #F7823S (CP)] and NSF CBET #1159735 (BF, CP). CB receives support from a University of Florida Graduate School Fellowship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

#### ACKNOWLEDGMENTS

We thank the members of the Computational Neuromechanics and Neural Control of Movement labs for their support, especially Nick Bianco and Andrew Meyer for methodological consultation, and Eric Walker and Virginia Little for helpful comments on previous versions of this manuscript. A portion of this work has been presented in abstract form at the International Society of Electrophysiology and Kinesiology XX Congress in July 2014 in Rome, Italy.

stroke survivors. J. Neurophysiol. 109, 768–781. doi: 10.1152/jn.006 70.2012

**Conflict of Interest Statement:** 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.

Copyright © 2017 Banks, Pai, McGuirk, Fregly and Patten. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Mobility as the Purpose of Postural Control

Charlotte Le Mouel\* and Romain Brette

Institut National de la Santé et de la Recherche Médicale, Centre National de la Recherche Scientifique, Institut de la Vision, Sorbonne Universités, Université Pierre et Marie Curie, Univ. Paris 06, Paris, France

Counteracting the destabilizing force of gravity is usually considered to be the main purpose of postural control. However, from the consideration of the mechanical requirements for movement, we argue that posture is adjusted in view of providing impetus for movement. Thus, we show that the posture that is usually adopted in quiet standing in fact allows torque for potential movement. Moreover, when performing a movement—either voluntarily or in response to an external perturbation—we show that the postural adjustments are organized both spatially and temporally so as to provide the required torque for the movement. Thus, when movement is performed skillfully, the force of gravity is not counteracted but actually used to provide impetus to movement. This ability to move one's weight so as to exploit the torque of gravity seems to be dependent on development and skill learning, and is impaired in aging.

Keywords: neuromechanics, posture, balance, motor control, movement

## INTRODUCTION

#### Edited by:

Matthew Tresch, Northwestern University, United States

#### Reviewed by:

Robert H. Lee, Emory University, United States Daya Shankar Gupta, Camden County College, United States

#### \*Correspondence: Charlotte Le Mouel<EMAIL_ADDRESS>
Received: 23 February 2017 Accepted: 05 July 2017 Published: 27 July 2017

#### Citation:

Le Mouel C and Brette R (2017) Mobility as the Purpose of Postural Control. Front. Comput. Neurosci. 11:67. doi: 10.3389/fncom.2017.00067 The position of the center of mass (CoM) is adjusted by the central nervous system during quiet standing (Winter et al., 1998; Sasagawa et al., 2009), in reaction to perturbations (Horak and Nashner, 1986), and in voluntary movement (Cordo and Nashner, 1982; Pedotti et al., 1989; Lee et al., 1990). The traditional theory is that the purpose of this postural control is to immobilize the CoM despite movement and external perturbations (Nashner et al., 1989; Massion et al., 2004; Horak, 2006; Bouisset and Do, 2008). We will refer to this theory as the immobility theory. The underlying assumption is that, because of gravity, standing is unstable. Therefore, if the CoM is displaced from its equilibrium position, then the displacement must be counteracted by postural adjustments, so as to return the CoM to its equilibrium position, otherwise the person will inevitably fall. As argued by Hasan (2005), this notion stems from an analysis of how linear systems respond to perturbations: in linear systems, if deviations from the unique equilibrium position are not corrected, then they grow exponentially. Balance (the ability to prevent falling), is therefore assumed to be equivalent to stabilization, in the strict sense of immobilizing the CoM at a unique equilibrium position by counteracting any displacement away from this position. From this assumption, it follows that moving poses a threat to balance, since any voluntary movement might displace the CoM. This theory has motivated a large body of experiments, performed over the last thirty years, in which a subject is asked to perform a movement of the upper body, while their muscle activity is being recorded (Cordo and Nashner, 1982; Crenna et al., 1987; Pedotti et al., 1989; Lee et al., 1990). In these experiments, a change in the contraction of the lower leg muscles is systematically observed, and this change often precedes the contraction of the upper body muscles. This is interpreted by saying that movement of the upper body might displace the CoM, and must therefore be counteracted by the contraction of the lower leg muscles so as to immobilize the CoM despite movement.

We will argue however that the equivalence between balance and immobilization does not hold for human postural control, and that these postural responses should be understood as providing the impetus for the movement. We will indeed show that during quiet standing, voluntary movement, and in reaction to perturbations, the position of the CoM is not immobilized at a unique equilibrium position, but on the contrary adjusted so as to use the torque of one's own weight, either to counteract external forces so as to maintain balance, or to provide impetus for voluntary movement. We therefore develop an alternative to the immobility theory. We propose that the purpose of postural control is mobility, the ability to produce appropriate impetus by adjusting the position of the CoM. We will refer to this theory as the mobility theory.

We will first show that the posture which is typically adopted in quiet standing allows for one's weight to be used to provide impetus to potential movement, and that when the direction of the movement to be performed can be anticipated, the position of the CoM during stance is shifted in that direction. Secondly, we will show that, during voluntary movement, postural adjustments which are traditionally thought of as immobilizing the CoM despite movement should on the contrary be interpreted as displacing the CoM at the initiation of the movement, so that one's own weight can be used to provide impetus to the movement. Finally, we will show that this ability to use displace one's weight, rather than immobilize it, plays a crucial role when balance is upset by external forces.

#### ADJUSTMENT OF POSTURE DURING STANCE

#### The Standing Posture Allows For Mobility The Standing Posture Requires Tonic Muscular Contraction

When someone is asked to stand quietly, without further instructions, they typically maintain their CoM vertically aligned with the middle of the foot, a few centimeters forwards of the ankle joint (Schieppati et al., 1994). However, when requested to do so, a young, healthy person can maintain their CoM at positions up to 40% of their foot length forwards of its typical position, and up to 20% backwards (Schieppati et al., 1994). There is therefore no unique equilibrium position for the CoM in quiet standing, since a young, healthy person can maintain a range of standing postures without this posing a threat to balance.

If the position of the CoM were controlled only in view of counteracting the torque of one's weight, then it would be most appropriate to place it vertically above the ankles, such that weight would exert no torque (**Figure 1A**). This position can indeed be maintained with minimum lower leg muscle contraction (Schieppati et al., 1994). However, when no instructions are given, subjects maintain their CoM vertically aligned with the middle of the foot, a few centimeters forwards of the ankle joint (**Figure 1B**), so that the weight exerts a forwards torque. In order to maintain this posture, an equivalent backwards torque must be exerted by the ground reaction force (see Section 1.1 in Appendix). As developed in the Appendix (Section 1.2), the torque of the ground reaction force is determined by the contraction of the lower leg muscles. Indeed, if we consider the forces acting on the foot, the weight of the body, carried by the skeleton, is applied at the ankle and therefore exerts no torque. The ground prevents the foot from turning, therefore the ground reaction torque instantly opposes the torque exerted by the lower leg muscles onto the foot (**Figure 2**). Maintaining a standing posture with the CoM forwards of the ankles therefore requires tonic contraction of the calf muscles (**Figure 1B**, Schieppati et al., 1994). The normal standing posture is therefore not the most economical in terms of muscular contraction.

#### The Standing Posture Allows Torque for Movement

Why would subjects actively maintain their CoM forwards of the ankles in quiet standing if this is not efficient? We suggest that this allows them to use their own weight for initiating forwards movements. Forwards torque for movement can only be induced by the external forces: the person's weight and the ground reaction force. As we have shown (Appendix Section 1.2) the ground reaction torque instantly follows the torque of the lower leg muscles. However, this torque is limited. Indeed, as long as the person neither jumps up nor collapses, the ground reaction force has the same magnitude as the person's weight. Its torque is therefore the product of the weight, and the distance between the ankle and the point of application of the ground reaction force, called the center of pressure and noted CoP. Thus, contracting the calf muscles (gastrocnemius and soleus) shifts the CoP forwards of the ankle (**Figure 2A**), and contracting the shin muscle (tibialis anterior) shifts the CoP backwards of the ankle (**Figure 2C**), but the CoP can only move within the limited range of the foot (see Appendix Section 1.3 for further detail).

The net torque is proportional to the distance between the CoM and the CoP. Whereas, the CoP moves instantly when the forces exerted by the muscles change, but can only move within the limited range of the foot, the position of the CoM on the other hand, does not change instantly when the forces exerted by the muscles change. This first requires the sum of the external forces to accelerate the CoM. Displacements of the CoM therefore occur more slowly than displacements of the CoP (as seen for example in Burleigh et al., 1994). Thus, the initial

then the CoP is below the ankle. (C) The torque around the ankles exerted by the shin muscle onto the foot is instantly compensated for by a backwards

net torque that can be produced, either for opposing external perturbation forces or for voluntary movement, is limited by the initial position of the CoM (see Appendix Section 1.4 for further detail).

When initiating fast forwards movements, either starting to walk (Burleigh et al., 1994) or movements performed with the feet in place such as leaning forwards (Crenna et al., 1987) or rising onto one's toes (Nardone and Schieppati, 1988), the CoP is first brought toward the heel by inhibiting the calf muscle contraction and contracting the shin muscle (Crenna et al., 1987; Nardone and Schieppati, 1988; Burleigh et al., 1994). If the CoM were initially above the ankle, this would produce little initial forwards torque (**Figure 3A**), whereas with the CoM forwards of the ankle this produces larger torque (**Figure 3B**). Maintaining the CoM forwards of the ankle thus allows one's own weight to be used for initiating forwards movement. Maintaining the CoM in the middle of the foot allows for either forwards or backwards initial torque to be induced by changes in the forces of the lower leg muscles.

shift of the CoP.

## The Standing Posture Is Actively Maintained

This position of the CoM is precisely and actively maintained on a short timescale, with small adjustments of the CoP in quiet standing serving to immobilize the CoM at this position (Winter et al., 1998). Moreover, the tonic contraction of the calf muscles is adjusted when standing on different slopes so as to maintain the CoM aligned with the middle of the foot (**Figure 4A**, Sasagawa et al., 2009). This precise positioning is also maintained at the longer timescales of growth and aging. Indeed, the curvature of the spine and trunk increases with aging (red line in **Figure 4B**, Schwab et al., 2006), and the position of the CoM is maintained across people with different trunk curvatures by shifting the position of the pelvis relative to the heels (**Figure 4B**, Schwab et al., 2006; Lafage et al., 2008).

Moreover, this forwards position of the CoM emerges with skill learning. Thus, Clément and Rézette (1985) observed acrobats at various competitive levels performing handstands. All the acrobats were able to maintain their balance in the upside-down posture, however they did so in different ways. The acrobats at lower competitive levels maintained their mean CoP a few millimeters forwards of their wrist; they could therefore maintain their posture with very little tonic contraction in the arm muscles (**Figure 4C**, left). The acrobats at higher competitive levels maintained their mean CoP more forwards of their wrists, with the acrobat at the highest level maintaining his mean CoP 3 cm forwards of his wrists; this posture requires tonic contraction of the wrist extensors (**Figure 4C**, right).

Thus, the standing posture is actively adjusted so as maintain the CoM above the middle of the foot (and above the middle of the hand in handstands). Contrary to the immobility theory, this position is not a unique equilibrium point, since a variety of standing postures can be maintained without this leading to a loss of balance. According to the mobility theory, this position is maintained because it allows for torque of the appropriate direction to be produced at short notice, even when this direction cannot be anticipated. This may be useful both for opposing external perturbations and for initiating voluntary movements.

## The Standing Posture Is Adjusted in Anticipation of Movement

When the direction of the appropriate torque can be anticipated, the mobility theory predicts that the CoM would be displaced in that direction in anticipation of the movement. Such a shift can indeed be induced experimentally, either by challenging someone's balance in a predictable direction, or by indicating in advance the direction of a voluntary movement to be performed.

Someone's balance can be challenged by having them stand facing the edge of the platform they are on. According to the immobility theory, this should lead, if anything, to an even more stringent immobilization of the CoM at its equilibrium position,

FIGURE 4 | Adaptation of the position of the CoM. (A) The tonic calf muscle contraction decreases when going from a slope with the toes down (left panel), to a flat slope (middle panel), to a slope with the toes up (right panel) such that position of the CoM is maintained vertically aligned with the middle of the feet. (B) People of different trunk curvatures maintain their pelvis at different distances from the heel line (vertical line above the heel), such that the CoM line (vertical line passing through the CoM) is at the same distance from the heel line. (C) Left panel: Acrobats at lower competitive levels maintain their CoP and CoM aligned with their wrist without tonic contraction of their wrist extensors. Right panel: Acrobats at higher competitive levels maintain their CoP and CoM forwards of their wrist, through tonic contraction of their wrist extensors.

Le Mouel and Brette Poised for Movement

but what is observed is that the CoM is shifted slightly backwards (**Figure 5A**, Carpenter et al., 2001). This is in accordance with the mobility theory, since it increases the person's capacity for producing backwards torque, in the eventuality that they might be subjected to a forwards push. In the experiment, the person's balance was not challenged beyond placing them in front of a drop, which might explain why the shift in CoM position was rather small (less than a centimeter).

Another way of challenging someone's balance is to have them stand on a platform (**Figure 5B**) which is then translated backwards (**Figure 5C**). The person ends up with their CoM in a forward position relative to the feet. A commonly observed response to such a translation is to straighten up (Welch and Ting, 2014). This requires backwards torque, however their capacity for producing backwards torque is limited by the forwards position of their CoM (**Figure 5C**). If such a perturbation is repeated, then over a few trials, the person adjusts their quiet standing posture by shifting their CoM backwards by a few centimeters (**Figure 5D**, Welch and Ting, 2014). This is again in contradiction with the immobility theory, but in accordance with the mobility theory, since the backwards shift of the CoM increases the person's capacity to produce backwards torque for straightening up (**Figure 5E**). When the platform is repeatedly translated forwards, then the person shifts their CoM forwards (Welch and Ting, 2014).

The mobility theory predicts that the position of the CoM in quiet standing would also be shifted if the direction in which a voluntary movement to be performed could be anticipated. This occurs at the start of a race: in sprinting, the initial forwards acceleration is crucial in winning the race. Consistently with the mobility theory, the CoM in the starting position is shifted even beyond the toes by several tens of centimeters (Slawinski et al., 2010). This is achieved by placing the hands on the ground and having the hands carry some of the weight (**Figure 5F**). This ability to use one's own weight to produce torque for movement again seems to depend on skill learning. Indeed, in elite sprinters, the CoM is shifted 5 centimeters further forwards than for welltrained sprinters (Slawinski et al., 2010).

#### Summary

Thus, when the direction of the appropriate torque to be produced cannot be anticipated, the CoM is positioned at the middle of the feet, in a position which allows for both forwards and backwards torque to be produced. When the direction of the torque to be produced can be anticipated, then the standing posture is adjusted by shifting the CoM in that direction. This adaptation of the standing posture in view of movement seems to be dependent on learning.

## ADJUSTMENT OF POSTURE DURING VOLUNTARY MOVEMENT

According to the immobility theory, when a voluntary movement is being performed, postural control serves to immobilize the CoM despite the movement or the perturbation. The mobility theory predicts, on the contrary, that the position of the CoM

FIGURE 5 | Adjustment of the position of the CoM. (A) When a person stands facing a slope, they shift their CoM slightly backwards. (B,C) When someone stands normally (B) and the platform they stand on is shifted backwards, their CoM ends up far forwards of the ankle joints, which limits the net backwards torque for straightening up (C). (D,E) When a backwards perturbation is repeated, the person shifts their CoM backwards in quiet standing (D), which increases the net backwards torque for straightening up after the perturbation (E). (F) In the posture adopted before a sprint, the CoM is placed far forwards of the feet by having the arms carry some of the weight.

is adjusted so as to use the torque of weight for movement. It therefore predicts that muscular contractions are temporally organized so as to accelerate the CoM at the initiation of the movement in the appropriate direction for producing torque for movement.

## Initiation of Voluntary Movement

#### Pulling on a Handle

When someone pulls on a handle placed in front of them, the contraction of the arm muscles is preceded then accompanied by the contraction of the calf muscles (Cordo and Nashner, 1982; Lee et al., 1990). Cordo and Nashner (1982) suggest that this contraction of the calf muscles allows for the CoM to be immobilized despite the movement. However, in order for the CoM to be immobilized, the ground reaction torque would have to exactly compensate for the handle reaction torque throughout the movement, and this would notably require the calf and arm muscle contractions to be simultaneous (as in **Figure 6A**). On the contrary, the initial contraction of the calf muscles which is observed (Cordo and Nashner, 1982) accelerates the CoM backwards (**Figures 6B,C**, further details are provided in Appendix Section 2); and when the person is asked to pull harder on the handle, this initial period lasts longer, the calf muscle activation is stronger, and the initial backwards acceleration of the CoM is larger (Lee et al., 1990). This is in accordance with the mobility theory, since initially accelerating the CoM backwards allows one's own weight to be used to assist the movement (**Figure 6C**).

#### Leaning the Trunk

numbers 1 to 2.

When someone leans the trunk forwards, the contraction of the abdominal muscles is preceded then accompanied by the inhibition of calf muscle contraction and the contraction of the shin muscle (**Figures 7A–C**, Crenna et al., 1987). The CoM could in theory be immobilized if the shin and abdominal muscle contractions were simultaneous, such that the forwards acceleration of the CoM induced by the shin muscle contraction would compensate for the backwards acceleration of the CoM induced by the abdominals contraction (further details are provided in Appendix Section 2), as suggested by Alexandrov et al. (2001). However, these authors report an initial backwards displacement of the CoP (**Figure 7A**), followed by a forwards displacement of the CoM (**Figure 7B**), in accordance with the sequential muscular contraction observed by Crenna et al. (1987). This contradicts the immobility theory, but concords with the mobility theory's predictions.

Thus, postural responses should be considered as an integral part of the movement itself, since they provide the torque for the movement, first by shifting the CoP and secondly by accelerating the CoM through sequential muscle contraction (a more complete explanation can be found in Appendix Section 2).

#### Gait Initiation

Bouisset and Do (2008) distinguish between two types of anticipatory postural adjustments. For voluntary movements without a change in the basis of support, such as raising the arm, they provide a very classical interpretion for the displacement of the CoM which precedes the displacement of the arm. They present it as a counterperturbation whose purpose is to "counterbalance the disturbance to postural equilibrium due to the intentional forthcoming movement" (Bouisset and Zattara, 1981). However, for voluntary movements involving a change of the basis of support, such as walking, or rising onto one's toes,

Frontiers in Computational Neuroscience | www.frontiersin.org

CoM backwards, allowing for a larger net backwards torque to be exerted during the handle pull (C). The sequence of activation of the muscles is indicated by the

FIGURE 7 | Leaning the trunk. When control subjects perform fast forwards leaning, the initial contraction of the shin muscle (A) accelerates the CoM forwards, thus allowing for more net forwards torque during the subsequent contraction of the ventral muscles (B), which enables the person to lean the trunk (C). When control subjects perform fast backwards leaning, the dorsal muscles contract simultaneously (D), which increases backwards rotational (Continued)

#### FIGURE 7 | Continued

momentum without translating the CoM (E). When gymnasts perform fast backwards leaning, the initial contraction of the calf muscles (F) accelerates the CoM backwards, thus allowing for more net backwards torque during the subsequent contraction of the dorsal muscles (G) which enables the gymnast to lean the trunk (H). The sequence of activation of the muscles is indicated by the numbers 1 to 3.

they present anticipatory postural adjustments as a perturbation involved in "body weight transfer" (Do et al., 1991).

We propose that in movements with or without a change in the basis of support, anticipatory postural adjustments play the same role of moving the CoM in order to provide impetus for movement. Indeed, the changes in posture which precede walking are organized in the same way as those which precede pulling on a handle or leaning the trunk. Thus, when going from standing to walking, a few hundred milliseconds before the heel of the swing foot is raised, the calf muscles are silenced and the shin muscle contracts, which brings the CoP to the heels and accelerates the CoM forwards, even before the first step is taken (**Figure 3B**, Burleigh et al., 1994). This is in accordance with the mobility theory, since initially accelerating the CoM forwards allows one's own weight to be used to assist the movement. Indeed, this initial acceleration of the CoM is correlated with the speed reached at the end of the first step, and is larger if the person is asked to walk faster (Brenière et al., 1987).

#### The Ability to Use One's Weight for Movement Requires Practice

For walking, a movement which is learned very early on in life, the ability to displace the CoM at the initiation of the movement emerges over the course of development (Ledebt et al., 1998; Bril et al., 2015). The amplitude of the initial backwards shift of the CoP thus increases over the first several years of life as children learn to walk faster (Ledebt et al., 1998; Bril et al., 2015). It then decreases with age, and with certain neurological diseases such as Parkinson's disease (Halliday et al., 1998; Mancini et al., 2016).

For leaning the trunk, the sequential muscle contraction, which allows for the displacement of the CoM at the initiation of the movement, seems to be dependent on learning. Indeed, when control subjects are asked to lean backwards, a movement for which they presumably have less practice than leaning forwards, then the calf and dorsal trunk muscle contractions are simultaneous (**Figure 7D**), and the movement is performed twice as slowly as leaning forwards (Pedotti et al., 1989). This is presumably because the CoM was not displaced backwards (**Figure 7E**). However, when gymnasts are asked to lean backwards, then their calf muscles contract first, and they perform the movement faster than controls (**Figures 7F–H**, Pedotti et al., 1989). Moreover, the ability to displace one's CoM during movement seems to remain plastic throughout life, and to depend on the possibility to use one's weight to assist movement. Thus, when astronauts return from a several months journey in space (during which they could not use their weight to assist their movements), the forwards displacement of the CoM when leaning forwards is no longer observed (Baroni et al., 2001).

Finally, for movements requiring skill learning, the temporal coordination which enables using one's weight to provide impetus for movement seems to develop with skill learning. Thus, when learning a complex gymnastics skill, such as the swings under parallel bars, in bent inverted hang position (**Figure 8**), beginners swing their legs and arms in synchrony, whereas experts swing their legs out of phase with their arms, which allows them to use the work of their own weight to provide impetus to the swing (Delignières et al., 1998).

#### BALANCE REQUIRES MOBILITY RATHER THAN IMMOBILITY

According to the immobility theory, if postural control does not immobilize the CoM at a unique equilibrium position, then the person must fall (Nashner et al., 1989; Massion et al., 2004; Horak, 2006; Bouisset and Do, 2008). We have shown however that in quiet standing, people can keep their balance over a range of positions of the CoM (Schieppati et al., 1994), and actually displace their CoM when their balance is challenged in a predictable (Carpenter et al., 2001; Welch and Ting, 2014). Moreover, we have shown that in well-practiced movements, people accelerate their CoM at the initiation of the movement, without this leading to a loss of balance (Cordo and Nashner, 1982; Crenna et al., 1987; Pedotti et al., 1989; Lee et al., 1990). We will now show that the response to an external perturbation should be considered as a movement in its own right, and therefore also benefits from the ability to use one's weight for movement, rather than to immobilize it.

#### Responding to External Perturbations Straightening Up after a Platform Translation

When the platform on which someone stands is translated backwards, the CoM ends up in a forward position relative to the feet (**Figures 9A,B**), as seen in Section The standing posture is adjusted in anticipation of movement. A response which is commonly observed is to straighten up (Horak and Nashner, 1986). The backwards acceleration of the CoM is performed through a sequential contraction of the dorsal muscles, starting with the calf muscles (**Figure 9A**), then the dorsal thigh then dorsal trunk muscles (Horak and Nashner, 1986). This contraction pattern is usually not considered as an actual movement, since it moves the CoM closer to its initial position, in accordance with the immobility theory. However, we believe it should be considered as a movement in its own right. Indeed, straightening up after a platform translation requires producing the appropriate backwards torque. The sequential contraction pattern allows for the CoM to be initially accelerated backwards, which increases the net backwards torque for the movement. Further details are provided in Appendix—Horizontal acceleration of the CoM. Moreover, contrary to the immobility theory, returning the CoM to its initial position is not the only way of preventing a fall.

#### Stepping after a Platform Translation

Indeed, another response which is also commonly observed is to take a step forwards (Maki et al., 2003): the CoM is then not returned to its initial position, without this causing a loss of balance. This response takes advantage of the forwards position of the CoM, such that the CoM needs not be accelerated backwards, and indeed the initial calf muscle contraction and forwards CoP shift is much reduced (**Figure 9B**) compared to when the person straightens up (**Figure 9A**); nor does the CoM need to be accelerated forwards, and indeed the shin muscle contraction lasts much less long and the backwards shift of the CoP is much smaller (**Figure 9C**, Burleigh et al., 1994) than when the person takes a step without the platform translation (**Figure 9D**).

## Emergence Over Development and Impairment with Aging

The ability to mobilize one's weight emerges over development. Thus, when straightening up after a backwards platform translation, both the systematic recruitment of the dorsal muscles and their temporal sequencing emerge during development. They are not observed in pre-walking infants, but are seen in children with a few years' walking experience (Burtner et al., 1998).

This ability is then deteriorated with aging, and with Parkinson's disease. The elderly, and even more so Parkinsonian patients, are less capable of moving their CoM, either when asked to adjust their quiet standing posture by leaning forwards or backwards (Schieppati et al., 1994), or during voluntary movement, such as gait initiation (Halliday et al., 1998). They are however quite as capable as young healthy adults of remaining immobile in quiet standing (Schieppati et al., 1994), and adjust the position of their pelvis to compensate for trunk curvature

such that the CoM remains above the middle of the feet (Schwab et al., 2006). Nevertheless, they have a heightened risk of falling. Thus, although the elderly and Parkinsonian subjects are quite as capable as young adults of maintaining their CoM immobile during quiet standing, we suggest that their higher risk of falling is due to a limited capacity to move when this becomes necessary to prevent a fall. Therefore, not only is immobilizing the CoM unnecessary for balance, it moreover seems that balance benefits from the ability to move one's CoM. This suggests that efficient balance training for the elderly can be achieved by practicing mobility (Xu et al., 2005).

#### DISCUSSION

## Posture Is Adjusted in View of Mobility Rather than Immobility

Although the position of the CoM is adjusted by the nervous system, this postural control does not serve to immobilize the CoM. On the contrary, the position of the CoM is adjusted so as to use the torque of one's own weight both for selfinitiated movements and for responding to external perturbation forces.

Thus, in quiet standing, when the direction of the torque to be produced cannot be anticipated, the CoM is maintained above the middle of the foot (Schieppati et al., 1994), allowing for the torque of one's weight to be used both for forwards and backwards movements. This position is actively maintained despite short-term changes in slope (Sasagawa et al., 2009) or long-term changes in trunk curvature (Schwab et al., 2006). However, when the direction of the torque to be produced can be anticipated, then the CoM is shifted in that direction. There is thus a small backwards shift of the CoM when someone is placed in front of a drop (Carpenter et al., 2001), or on a platform which is repeatedly translated backwards (Welch and Ting, 2014). Skill learning leads to much larger shifts in the position of the CoM, with the CoM placed forwards of the feet in anticipation of sprinting (Slawinski et al., 2010).

Moreover, during movement, we have shown that the postural responses which were thought to immobilize the CoM despite movement are actually temporally organized so as to accelerate the CoM at the initiation of the movement, in the appropriate direction such that the torque of one's weight can be used for the movement (Cordo and Nashner, 1982; Crenna et al., 1987; Pedotti et al., 1989; Lee et al., 1990). These postural responses should therefore be understood as providing impetus to the movement.

Finally, we have shown that in order to respond effectively to external perturbation forces, the CoM need not be immobilized, since the person can take a step (Maki et al., 2003). When the person straightens up without taking a step (Horak and Nashner, 1986), this requires producing forces to counteract the external perturbation, and may benefit from the ability to mobilize one's CoM rather than immobilize it. Balance therefore requires mobility rather than immobility.

FIGURE 9 | Response to platform translation. When straightening up after a platform translation (A), the initial contraction of the calf muscle accelerates the CoM backwards which increases the potential net backwards torque. When stepping forwards in response to a platform translation (B,C), this initial calf muscle contraction is reduced (B). Then, the shin muscle contracts (C). This shin muscle contraction is smaller than the initial contraction of the shin muscle which accelerates the CoM forwards when the person steps forwards without a platform translation (D). The sequence of activation of the muscles is indicated by the numbers 1 to 2.

## Mobility Emerges Through Development and Skill Learning

The ability to use one's weight for movement emerges through development and skill learning, and remains plastic throughout life. The appropriate temporal organization of muscular contraction emerges during development both for walking and for balancing responses (Burtner et al., 1998; Ledebt et al., 1998). It is not observed for less practiced movements, such as when control subjects lean the trunk backwards (Crenna et al., 1987). The extent to which the CoM can be mobilized seems to depend on the level of skill: thus, both for sprinters at the initiation of a race (Slawinski et al., 2010) and acrobats performing handstands (Clément and Rézette, 1985), elite athletes place their CoM further forwards than well-trained athletes. Future work should address the following questions: how is this ability learned through development and practice? Does the impairment of this ability in aging result from a lack of practice, and could this

## REFERENCES

ability be maintained during aging through appropriate training regimes?

## AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

This work was supported by Agence Nationale de la Recherche (ANR-14-CE13-0003).

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fncom. 2017.00067/full#supplementary-material

extensors: quiet standing on inclined surfaces. Exp. Brain Res. 196, 537–544. doi: 10.1007/s00221-009-1876-4

**Conflict of Interest Statement:** 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.

Copyright © 2017 Le Mouel and Brette. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Propagation of Movement Variability in Time: A Methodological Approach for Discrete Movements with Multiple Degrees of Freedom

Melanie Krüger 1, 2, 3 \*, Andreas Straube2 and Thomas Eggert 2

1 Sensorimotor Neuroscience and Ageing Research Laboratory, School of Medicine, University of Tasmania, Hobart, TAS, Australia, 2 Department of Neurology, University Hospital Munich Großhadern, Munich, Germany, 3 Department of Sport and Health Sciences, Technical University of Munich, Munich, Germany

In recent years, theory-building in motor neuroscience and our understanding of the synergistic control of the redundant human motor system has significantly profited from the emergence of a range of different mathematical approaches to analyze the structure of movement variability. Approaches such as the Uncontrolled Manifold method or the Noise-Tolerance-Covariance decomposition method allow to detect and interpret changes in movement coordination due to e.g., learning, external task constraints or disease, by analyzing the structure of within-subject, inter-trial movement variability. Whereas, for cyclical movements (e.g., locomotion), mathematical approaches exist to investigate the propagation of movement variability in time (e.g., time series analysis), similar approaches are missing for discrete, goal-directed movements, such as reaching. Here, we propose canonical correlation analysis as a suitable method to analyze the propagation of within-subject variability across different time points during the execution of discrete movements. While similar analyses have already been applied for discrete movements with only one degree of freedom (DoF; e.g., Pearson's product-moment correlation), canonical correlation analysis allows to evaluate the coupling of inter-trial variability across different time points along the movement trajectory for multiple DoF-effector systems, such as the arm. The theoretical analysis is illustrated by empirical data from a study on reaching movements under normal and disturbed proprioception. The results show increased movement duration, decreased movement amplitude, as well as altered movement coordination under ischemia, which results in a reduced complexity of movement control. Movement endpoint variability is not increased under ischemia. This suggests that healthy adults are able to immediately and efficiently adjust the control of complex reaching movements to compensate for the loss of proprioceptive information. Further, it is shown that, by using canonical correlation analysis, alterations in movement coordination that indicate changes in the control strategy concerning the use of motor redundancy can be detected, which represents an important methodical advance in the context of neuromechanics.

Keywords: movement coordination, variability, canonical correlation, reaching, sensory loss

#### Edited by:

Matthew Tresch, Northwestern University, United States

#### Reviewed by:

J. Michael Herrmann, University of Edinburgh, United Kingdom Daya Shankar Gupta, Camden County College, United States

#### \*Correspondence:

Melanie Krüger<EMAIL_ADDRESS>
Received: 31 December 2016 Accepted: 26 September 2017 Published: 13 October 2017

#### Citation:

Krüger M, Straube A and Eggert T (2017) The Propagation of Movement Variability in Time: A Methodological Approach for Discrete Movements with Multiple Degrees of Freedom. Front. Comput. Neurosci. 11:93. doi: 10.3389/fncom.2017.00093

## INTRODUCTION

Analyzing movement variability to gain insights into movement planning and control processes has been in the focus of researchers in the field of (computational) human motor control ever since Bernstein's (1967) famous Blacksmith example, which describes that, even in highly skilled movements with high outcome stability, execution variability can be observed across repetitions. Originating from that observation, numerous studies have been conducted to investigate changes in the amount and structure of movement variability with changing external task constraints (Gera et al., 2010; van der Steen and Bongers, 2011; Krüger et al., 2012) or under manipulated availability of sensory information during movement planning and execution (Tseng et al., 2002; Krüger et al., 2011), showing that variability in movement execution is an inherent characteristic of human performance. In the last recent years, this research has significantly profited from the emergence of a range of different mathematical approaches to analyze the structure of movement variability. Approaches such as the Uncontrolled Manifold method (Scholz and Schöner, 1999) or the Noise-Tolerance-Covariance decomposition method (Müller and Sternad, 2004) allow to detect and interpret changes in movement coordination by analyzing the structure of withinsubject, inter-trial movement variability. Referring to these methods, it was shown that changes in the structure of movement variability can be related to changes in movement planning and control processes due to learning, aging, and pathology (Cirstea and Levin, 2000; Müller and Sternad, 2009; Stergiou and Decker, 2011; Krüger et al., 2013).

In that context, particularly the investigation of the time course of movement variability during movement execution has stimulated theory-building in motor neuroscience and advanced our understanding of the synergistic control of the redundant human motor system. Research on reaching and pointing movements (Domkin et al., 2002, 2005; Cohen and Sternad, 2009; Krüger et al., 2012) provided empirical support for theories of motor control which postulate that the human motor system exploits its inherent redundancy to cost-optimize movement execution, such that only variability in task-relevant dimensions is minimized, a principle referred to as "minimum intervention principle" (Todorov and Jordan, 2002; Todorov, 2004).

For time course analyses of movement variability, different time points between movement start and end, usually corresponding to either certain percent of the normalized time between movement start and end or to distinct time points in the movement, such as e.g., time point of maximum velocity, are examined. Following the definition of relevant time points between movement start and end, two principally different approaches can be followed to analyze the time course of intertrial movement variability: first, the variance structure of the effector position at a single time point during the movement can be analyzed across many movement repetitions, and conclusions can be drawn from changes of this variance structure between different time points. As an example for this approach, a range of studies using the Uncontrolled Manifold method were able to show significant differences in the control of task-relevant and—irrelevant variability in the effector space across the time course of movement execution for, e.g., sit-to-stand, shooting, or goal-directed reaching movements (Scholz and Schöner, 1999; Scholz et al., 2000; Tseng et al., 2002). However, this approach focuses on the relations between variability in the effector and task space at one or multiple points in time, but does not allow direct conclusions about statistical coupling of effector variables across time.

In contrast to this first approach, the second approach focuses on the coupling of movement variability between different time points during movement execution. Importantly, while for cyclical movements (e.g., locomotion), mathematical approaches exist to investigate the propagation of movement variability in time (e.g., time series analysis, or Lyapunov exponent, Stergiou and Decker, 2011), similar approaches are missing for discrete, goal-directed movements with a redundant effector system, such as reaching with the arm. Here, we propose canonical correlation analysis as a suitable method to analyze the propagation of within-subject variability across different time points during the execution of discrete movements. While similar analyses have already been applied for discrete movements with only one degree of freedom (DoF; Messier and Kalaska, 1999; Richardson et al., 2011; Kuang and Gail, 2015; Eggert et al., 2016), canonical correlation analysis allows to evaluate the coupling of intertrial variability across different time points along the movement trajectory for multiple DoF-effector systems, such as the arm.

This method will be illustrated by empirical data from a study on reaching movements under normal and disturbed proprioception. Proprioception about joint positions is an important source of information for the control of complex reaching movements (Ghez and Sainburg, 1995; Bagesteiro et al., 2006). Studies on chronically deafferented patients suffering from severe peripheral sensory neuropathy showed impaired motor control of arm movements, including slowed movement execution (Gentilucci et al., 1994; Hepp-Reymond et al., 2009), increased movement variability (Gentilucci et al., 1994; Medina et al., 2010) and deteriorated movement coordination (Sainburg et al., 1993, 1995; Ghez and Sainburg, 1995). Studies of temporary peripheral deafferentation of healthy humans showed immediate adjustment to the loss of proprioception on a behavioral (Moisello et al., 2008, applying limb immobilization) and cortical level (Björkman et al., 2004a,b, applying a local anesthetic cream; Ziemann et al., 1998b, applying an ischemic nerve block). However, these studies mainly requested the production of simple motor tasks with a limited range of kinematic DoF. Studies on the production of complex motor behavior, such as reaching movements, are rare. Here, we investigated the effect of temporary proprioceptive deafferentation, induced by an ischemic block at the upper arm level, on the time course of joint angle variability of complex arm movements.

## METHODS

## Participants

Fifteen healthy volunteers (mean age ± SD: 26 ± 5 years; 8 female) participated in the study. All participants were righthand dominant as determined by the Edinburgh Handedness

Inventory (Oldfield, 1971) and had normal or corrected-tonormal vision. None of the participants had any record of neurological disorder. All participants were paid for their participation and had given written informed consent prior to participation. The experimental procedure was in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of the Medical Faculty of the Ludwig-Maximilians University Munich.

 A linear table track was mounted on the table, with a spherical object (reaching target, diameter: 80 mm) attached to it. Due to its geometric properties, the reaching target constrained final hand position but not final hand orientation. The size of the target forced the participants to grasp it with the whole hand, and not just with two fingers, which is why single finger motion was not of interest in the current study. The reaching target could be freely moved horizontally (in the fronto-parallel plane) between the bounds of the table track. These bounds (distance: 39 cm) were the two positions at which the reaching target could be located. The sitting position of the participants was adjusted so that: (a) trunk movement was not necessary to reach the target, and (b) body midline was centered to the table track. To minimize within-subject inter-trial variability due to differences in the initial position, the starting position was defined by a wooden lever, attached to the right side of the chair, which had to be grasped with the dominant right hand before each trial (see **Figure 1A**).

Joint angle motion of the arm in its seven degrees of freedom was recorded by an ultrasonic sound-emitting system (Zebris Medical, Isny, Germany). Six sound-emitting markers were attached to the arm and hand of the participant; each marker recorded at a frequency of 33 Hz. The following marker positions were chosen and are depicted in **Figure 1B**: marker 1 and 2 were attached to the metacarpophalangeal joints of the index (1) and little finger (2). The third marker was at the center of the wrist. Marker 4 and 5 were attached to the medial (4) or lateral (5) end of a bracelet directly above the elbow. The sixth marker was attached at the acromion. From those marker positions, the individual length of the participant's upper arm, lower arm, and hand could be determined. Based on these lengths, a geometrical model of the arm was created, as described in more detail below (see section Data Analysis). Further, the signal of the first marker was used to trigger the opening and closing of shutter glasses (Translucent Technologies, Toronto, Canada) that were used to prevent visual online control of the movement. The first contact with the reaching target was detected by changes in the electrical resistance between the participant and the target (sampled at 1 kHz).

## Procedure and Design

Participants repeatedly had to reach toward and grasp the reaching target with their dominant right hand. At the beginning of each trial, participants had to adopt the starting position (see **Figure 1A**). Subsequently, participants were instructed to press a button with their non-dominant hand, after which the target changed its position. After an acoustic go-signal, participants had to perform the reaching movement in a natural manner. To provoke the most natural movement behavior, participants were informed before movement recording that movement speed and reaction time were not of interest in the study. Shutter glasses occluded as soon as the participants started their movement, thus preventing visual online control of the movement. After the participants had grasped the target, the shutter glasses opened again and the participants returned to the starting position. A new trial was initiated by pressing the button again.

All participants participated in two experimental conditions in separate sessions, the order of which was counterbalanced across participants. Experimental sessions were separated by 1 to 2 days. In the first condition ("ischemia"), a customary blood pressure cuff was applied to the upper arm of the participants and inflated up to 150–160 mmHg (i.e., slightly above systolic blood pressure) to induce a transient ischemic block. Ischemic nerve block is an established experimental technique to study sensory control of movements (e.g., Glencross and Oldfield, 1975; Ziemann et al., 1998a). It is known to first affect the large, fast conducting afferent fibers, especially Ia afferents arising from the muscle spindle afferents (Fellows et al., 1993). In contrast to acute limb ischemia, a sudden decrease in limb perfusion caused by e.g., thrombosis, transient ischemic block is an incomplete block of limb perfusion caused by externally applied pressure and is non-threatening to the limb. Glencross and Oldfield (1975) showed that 20–25 min of ischemic nerve block results in complete dropout of finger sensation, and significant sensory decrease in wrist and elbow. Decrements in the exerted force were observed only after complete sensory dropout. In the current study, the duration of inflation was in a range of 20–25 min. This timeframe included 10 min of preparation to guarantee impairment in the global sensory afference, and a subsequent 10–15 min of movement recordings. Consequently, other effects of the ischemic nerve block, such as changes in producible muscle force (Björkman et al., 2004a), can be disregarded in the currentset-up because of the brevity of the ischemic block. Before movement recording started, the proprioceptive impairment was tested indirectly by assessing participants' touch sensitivity with von-Frey filaments (Marstock, Schriesheim, Germany, Rolke et al., 2006). On the back of the participants' hands it was tested which of the 12 logarithmically scaled filaments (range: 0.25–512 mN) participants were at least able to perceive. On average, participants were able to perceive a minimum pressure of 0.5 mN before the application of the blood pressure cuff (i.e., mean filament number ± SD: 1.91 ± 0.53). Participants' touch sensitivity had to be reduced by at least one filament (i.e., increase by a factor 2) before the experiment was continued. On average, participants perceived a minimum of 1 mN at the start of the movement recordings (i.e., mean filament number ± SD: 3.09 ± 0.70). This procedure allowed us to be sure about the effectiveness of the ischemic block. At the same time, the duration of preparation was minimized, which was of importance to prevent unwanted side effects of the ischemic block, as for example ischemic pain. The second experimental condition ("control") served as a control condition, executed identically but without inflated blood pressure cuff.

Two blocks of 40 trials in each block were recorded in each session (i.e., 80 trials per session). Each experimental block consisted of 20 trials of each of the two target positions, arranged in a random order to avoid predictability of the target position. Between the blocks a break of maximally 5 min was offered to avoid fatigue. Before movement recording started, participants were allowed to perform five trials to familiarize themselves with the experimental task and apparatus.

#### Analysis Data Analysis

Data analysis was calculated using Matlab 7.9.0 (Mathworks, Natick, USA) and was in line with earlier studies by our group (Krüger et al., 2011, 2012). In a first step, the seven joint angles of the arm were computed from the marker positions using a three-segment rigid body model, and expressed as seven consecutive Cardan angles. The order of the angles was as follows: two angles for the wrist (vertical, and horizontal), two angles for the elbow (torsion, and flexion), and three angles for the shoulder (torsion, horizontal, and vertical). The zero position of the arm was defined as the arm pointing straight forward with the elbow extended and the palm facing up. Based on that, positive joint angle indicated the following directions: vertical upward, horizontal rightward, and torsion clockwise. The vector containing the seven joint angles is hereafter referred to as arm posture. The position of the hand in space (i.e., 3D) was defined by the center of the two hand markers in world fixed Cartesian coordinates. In addition, the orientation of the hand in space was defined in Helmholtz coordinates relative to the external world.

Temporal and spatial movement characteristics were analyzed separately for each condition, participant, target position and trial. Overall movement duration was defined as the time between movement initiation and movement end. To determine movement initiation, movement start was defined as the time point at which the hand velocity first exceeded 10% of its peak velocity. Movement initiation was then determined by subtracting 10% of the acceleration time (i.e., the time between movement start and reaching peak velocity) from movement start. Movement end was defined as the last sample recorded before the first contact with the reaching target, as determined by the change in electrical resistance (see section Experimental Set-Up). Subsequently, duration of acceleration and duration of deceleration were calculated. In addition, peak velocity was analyzed. Thus, temporal characteristics of the reaching movements will be described by four measures: (1) overall movement duration, (2) duration of acceleration, (3) duration of deceleration, and (4) peak velocity. To determine spatial characteristics of the reaching movement, movement amplitudes were determined by calculating the absolute value of the difference between the maximum and minimum joint angle separately for each of the seven joint angles. Subsequently, mean movement amplitude was calculated as the average movement amplitude across the seven joint angles. In addition, to evaluate the changes in the diversion from shortest trajectory between starting and end position, the total path length in the 7D-joint space of the arm was calculated.

Within-subject inter-trial movement variability during the time course of movement execution and at movement end was analyzed separately for each condition, participant, and target position. On that account, the full temporal resolution of the joint angle motion was reduced to 10 equidistant sampling points between movement initiation and movement end. To account for small inter-trial variations in the actual starting position of the arm and in movement duration, a correction of the joint angle trajectories was calculated as described in detail in Krüger et al. (2011). Briefly, the within-subject deviations of the joint angles from their mean were submitted to a linear regression analysis with the predictor initial arm position and movement duration (i.e., 7 + 1 = 8 continuous predictor variables). Separate regression analysis were conducted for each participant, experimental condition, target position, and for each of the 10 sampling points, thus, containing the data of 40 trials. Subsequently, the joint angle deviations from the mean that were predicted by this linear model were subtracted from the actual joint angles. After this correction, the covariance matrix of the starting position (first sample) reduced to zero and was not considered in further analytical steps. Thus, the covariance matrix of the joint angles was analyzed at nine equidistant sampling points during the movement.

Afterwards, movement variability was analyzed by means of the two approaches described in the Introduction: first, analyzing the amount of variability at the nine single sampling points during movement execution, and second, analyzing the coupling of movement variability across different sampling points. To achieve the first, the square-root of the mean within-subject variance, averaged across the seven joint angles of the arm (hereafter referred to as: standard deviation of arm posture), was calculated. Further, the square root of the mean within-subject variance, averaged across its three dimensions was calculated for the task variables (a) hand position (standard deviation of hand position) and (b) hand orientation (standard deviation of hand orientation). In addition, the coupling between joint angles within the arm posture at a given sampling point was analyzed by calculating a principal component analysis on the 7 × 7 covariance matrix of the arm posture. Subsequently, the variances for each of the seven eigenvalues of the covariance matrix were averaged across sampling points, and the percentage of total variance explained by the first two eigenvalues was calculated. A relative increase of this percentage is closely related to a relative decrease of the number of kinematic DoF with respect to the mechanical DoF.

To achieve the second, the coupling between the arm posture at a given sampling point and the final arm posture was examined. On that account, the coupling between the arm posture during the movement and the final arm posture was assessed by canonical correlation analysis evaluating the percentage of inter-trial variance of the final arm posture that could be explained by the variance of arm posture at a given sampling point. The canonical correlation analysis returns a coefficient of determination, which equals the mean R 2 across the multiple regressions explaining the final arm posture as linear functions of the arm posture at a given sampling point.

#### Statistical Analysis

Statistical analysis was calculated using SPSS 9.0. Pairwise comparisons were calculated for the temporal and spatial measures of the reaching movements, the percentage of total variance explained by the first two eigenvalues, as well as for the coefficient of determination. A repeated measurement ANOVA with condition (control vs. ischemia) as the between-group factor, and sampling point as the repeated factor was calculated for the following dependent variables: (1) standard deviation of arm posture, (2) standard deviation of hand position, and (3) standard deviation of hand orientation. Bonferroni corrected pairwise comparisons were calculated for post-hoc analysis of significant interactions. A Greenhouse-Geisser adjustment was made if the sphericity assumption was rejected by Mauchly's sphericity test. Standard deviations of arm posture, hand position and hand orientation were tested for normal distribution with the Lilliefors-test. Data was normally distributed for both conditions and for almost all sampling points. The critical value for significance was set at p < 0.05. Participants were excluded from single analyses in case of data corruption, i.e., if the data matrix for each participant contained <10 valid trials in each condition and for each target position.

## RESULTS

Since the influence of target position on complex reaching movements was not of interest in the current study, and was already discussed elsewhere (Krüger et al., 2011, 2012), only the results for reaching toward the left target position will be presented here. Similar results were found for reaching movements toward the right target position, though in general the observed differences were smaller for the right target position as compared to the left target position.

#### Temporal Movement Characteristics

Overall movement duration was 778 ± 167 ms (mean ± SD) for the ischemia condition and 713 ± 142 ms for the control condition (see **Figure 2**). This difference was significant [paired t-test: t(14) = −3.55, p < 0.01] and based on a significantly increased duration of the acceleration phase under ischemia [403 ± 83 vs. 352 ± 84 ms, t(14) = −3.08, p < 0.01]. Neither duration of the deceleration phase (375 ± 116 vs. 360 ± 101 ms), nor peak velocity (1,068 ± 198.57 vs. 1,095 ± 195 mm/s) differed between the ischemia and control condition.

## Spatial Movement Characteristics

When reaching toward the target, trajectories for five out of the seven joint angles of the arm (shoulder torsion, shoulder horizontal, shoulder vertical, elbow torsion, and wrist horizontal) showed a continuous increase or decrease between movement initiation and movement end, with the trajectories slightly curved. For elbow flexion and wrist vertical, joint angle trajectories showed a reversal in movement direction during the movement. Under ischemia, total path length in the 7D-joint space was decreased by 15% (control: 40.8 ± 6.8◦ vs. ischemia: 34.3 ± 5.4◦ ). Associated with that, the participants' mean movement amplitude decreased significantly under ischemia as compared to the control condition [26.3 ± 4.1◦ vs. 31.3 ± 4.7◦ , t(14) = 5.32, p < 0.01, see **Figure 3A**]. This difference was especially pronounced in four of the seven joint angles: shoulder torsion [t(14) = 2.46, p = 0.03], shoulder vertical [t(14) = 2.95, p = 0.01], elbow torsion [t(14) = 3.93, p < 0.01], and elbow flexion [t(14) = 5.50, p < 0.01, see **Figure 3B**].

#### Movement Variability

First, movement variability at the different sampling points during movement execution was analyzed with respect to three

measures: (a) standard deviation of arm posture, (b) standard deviation of hand position, and (c) standard deviation of hand orientation. The amount of movement variability did not differ between the two experimental conditions (i.e., no significant main effect of experimental condition) for any of the three measures neither across the nine sampling points nor at movement end. However, for each of the three measures, a significant main effect of sampling point became evident: (a) F(2.39, 23.92) = 21.21, p < 0.01, (b) F(2.36, 23.62) = 53.35, p < 0.01, and (c) F(2.48, 24.83) = 22.93, p < 0.01. In all cases, movement variability increased until the middle of the movement and decreased afterwards. Variability was smallest at movement initiation and on an intermediate level at movement end (see **Figure 4**).

The interaction of experimental condition × sampling point was significant for standard deviation of hand position. Qualitatively, this effect became evident as a weaker modulation of hand position variability across the nine sampling points in the ischemia condition (see **Figure 4B**). Post-hoc analysis revealed that, under ischemia, only the first two sampling points differed largely from the other sampling points, whereas in the control condition almost all sampling points differed significantly from each other (see **Table 1**). No other effects reached the level of significance.

The coupling of joint angles within the arm posture was analyzed using a principal component analysis applied to the inter-trial 7 × 7 covariance matrix of the arm posture. Under ischemia, the first two eigenvalues explained 88.90 ± 2.44% of total joint angle variance compared with 83.40 ± 2.27% in the control condition (see **Figure 5A** for group mean and **Figure 5B** for a representative participant). This difference was significant [t(8) = −18.43, p < 0.01].

The coupling of the arm posture was analyzed using the coefficient of determination between of final arm posture with respect to the arm posture during movement execution. As a matter of course, the coefficient of determination increased toward movement end and finally reached the level of 1 (see

**Figure 5D** for a representative participant). The coefficient of determination of final arm posture with respect to the variance of the arm posture at the first sampling point was smaller in the control condition (R 2 = 0.4) than under ischemia (R 2 = 0.7). Consequently, the subsequent increase in the coefficient of determination up to the value 1 at movement end was steeper in the control condition than under ischemia. For group comparison, only the coefficient of determination with respect to the fifth sampling point, when the standard deviation of arm posture was maximal, was analyzed. Under ischemia the coefficient of determination was significantly higher than in the control condition [R 2 : 0.82 ± 0.18 vs. 0.54 ± 0.06, t(7) = −3.89, p < 0.01, see **Figure 5C**].

#### DISCUSSION

In the current study, we introduced a method to investigate the coupling of joint angle variability across the time course of discrete, goal-directed, natural reaching movements. The method

FIGURE 4 | Movement variability. (A) Standard deviation of arm posture, representing the mean across participants, is shown. Shaded areas represent the respective confidence intervals. (B) Standard deviation of mean hand position (+ confidence interval) is depicted. Hand position variability was less modulated under ischemia. (C) Standard deviation of mean hand orientation and the respective confidence interval is shown.

TABLE 1 | Post-hoc analysis for the significant interaction Condition × Sampling point for the measure: Standard deviation of hand position.

Data in the upper right half of the table represents p-values of significant pairwise comparisons of single sampling points for the ischemia-condition. Data in the lower left half of the table represents p-values of significant pairwise comparisons of single sampling points for the control condition. Bonferroni correction for multiple comparisons was applied to all calculations. Asterisks indicate close to significant interactions.

was exemplified on a dataset that was collected to study the influence of temporary proprioceptive deafferentation on the control of a complex reaching movement. In the following, the main outcomes of the study will be discussed first, followed by a discussion on the canonical correlation method introduced to analyze the propagation of movement variability across time.

#### Adjustment of Movement Duration

Movement duration was increased by the ischemia as a result of increased acceleration duration. The influence of proprioception on the duration of acceleration was already recognized by Bagesteiro et al. (2006) and associated with sensory-based online-correction of the movement. Movement's peak velocity was not increased under ischemia. Increased duration of acceleration without increased peak velocity indicates decreased peak acceleration and, consequently, decreased peak force. A reduction in total force applied during movement execution is accompanied by a reduction in signal-dependent noise (Harris and Wolpert, 1998). This may be advantageous under ischemia, as the precision of movement planning is of greater importance when movement online-control based on proprioceptive feedback is impaired. Our results suggest that healthy participants are able to immediately and efficiently adjust the precision of movement planning to the lack of proprioceptive information in elbow and wrist.

## Adjustment of Movement Amplitude

Movement amplitude was decreased under ischemia. Importantly, this was not only true for joints distal to the applied blood pressure cuff (i.e., elbow torsion and elbow flexion), which were directly affected by the ischemic block, but also for two joint angles proximal to the cuff (i.e., shoulder

torsion and shoulder vertical), which were not directly affected by the ischemia. In combination with the finding of stronger inter-joint coupling under ischemia, this suggests a more global change in the strategy of joint angle coordination involving all joints of the arm to compensate for the ischemia. A reason for planning a reaching movement with decreased mean movement amplitude may be the associated decrease in the signal-to-noise ratio (Harris and Wolpert, 1998), facilitating the control of movement endpoint variability. This assumption is also supported by Fitts' Law (Fitts, 1954), which describes the relationship between movement amplitude, movement duration and movement accuracy. According to this law, in order to keep movement endpoint variability constant in a task with increased task difficulty, movement duration and/or movement amplitude must be adjusted. Assuming that the ischemia may have increased task difficulty, as an

arm posture variance for each sampling point is shown for a representative participant for both conditions.

important source of sensory information was disabled, planning a movement with decreased movement amplitude and increased movement duration may have allowed the participants to keep movement endpoint variability constant, as observed in our study.

## Adjustment of Movement Variability

Another important finding of our study was that the modulation of hand position variability during movement execution was altered under ischemia in such a way that the initial increase and subsequent decrease of hand position variability was less pronounced. The increase-decrease pattern of movement variability was already described in earlier studies by our group (Krüger et al., 2011, 2012) and is a sign of successful minimization of variance at movement end. It indicates that signal-dependent noise (Harris and Wolpert, 1998), introduced by forces during the acceleration period, is successfully compensated by feedback control acting primarily during the deceleration phase (Elliott et al., 2001, 2010; Eggert et al., 2016). The fact that this increase-decrease pattern of hand position variability is less pronounced under ischemia (see **Figure 4**) is probably related to both reduced acceleration forces, resulting in a reduced increase of variability, and impaired proprioceptive feedback, resulting in a reduced decrease of variability. Interestingly, both of these changes compensated for each other in such a way that endpoint variability was almost identical in the control condition and under ischemia. This is in contrast to findings of studies with chronically deafferented patients (Gentilucci et al., 1994; Gordon et al., 1995; Nougier et al., 1996; Medina et al., 2009) and reflects the ability of the motor control system of healthy participants to immediately and efficiently adjust to the loss of proprioceptive information in parts of the effector.

#### Adjustments of Movement Coordination

Movement coordination was altered under ischemia, a finding similar to that observed in studies on deafferented patients (Sainburg et al., 1993; Ghez and Sainburg, 1995; Sarlegna et al., 2006). In the current study, the alterations in movement coordination became evident for the coupling between single joints of a specific arm posture as well as for the coupling of arm posture during movement execution with that at movement end: for both parameters the coupling was stronger under ischemia. Increasing the strength of joint angle coupling under ischemia, i.e., increasing the synergistic coordination of the redundant DoF at the same time point and across time points, may reflect a change in the control strategy concerning the way motor redundancy is used under impaired proprioceptive feedback. Alternatively, it may reflect the limited capacity of the brain to plan and control coordinated movement with a naturally higher number of DoF with decreased proprioceptive feedback. Similar to that assumption, Gupta (2014) highlighted the relevance of feedback information for the precise temporal control of complex motor actions. Independent of which explanation holds true, the observed effect of stronger joint angle coupling under ischemia can be interpreted as a reduction of the number of kinematic DoF of the redundant effectorsystem arm and consequently as a facilitation of its onlinecontrol.

### Methodological Considerations

In this study, canonical correlation analysis was introduced as a method to investigate the propagation of movement variability for discrete movements involving a redundant effector system. It was shown that, by using canonical correlation analysis, alterations in movement coordination that indicate changes in the control strategy concerning the use of motor redundancy could be detected. To the best of our knowledge, this is the first time such an approach is suggested for a multiple DoF-effector system in the context of discrete movements. Previous approaches to capture the propagation of movement variability for discrete movements have studied either eye movements (West et al., 2009; Richardson et al., 2011; Eggert et al., 2016) or arm movements with a limited number of DoF (Messier and Kalaska, 1999; Heath, 2005; Kuang and Gail, 2015). However, to advance our understanding of the human motor control system in its complexity, natural movements as used in this study have to be analyzed. In this context, the introduction of canonical correlation analysis as a suitable approach represents an important methodical advance in the context of neuromechanics. One aspect contributing to the importance of the introduced method is that it does not substitute or extend other approaches that have been previously suggested to identify structure in movement variability, but represents a new approach with the potential to also broaden the application of already existing methods and to increase their significance: In this study, it was analyzed how strongly final arm posture variability was determined by overall variability of arm posture during movement execution. Previous studies identified different components of overall movement variability (Müller and Sternad, 2004, 2009) and the relevance of independent components of effector variability for task variability (Scholz and Schöner, 1999; Latash et al., 2002). Following this line of thinking, future research could use canonical correlation analysis to investigate the propagation of certain components of overall effector variability in time.

## CONCLUSIONS

In this study, we introduced canonical correlation analysis as a method to investigate the temporal propagation of movement variability in reaching movements under normal or impaired proprioceptive feedback. As general findings, we found increased movement duration due to increased acceleration duration, decreased movement amplitude, as well as changes in movement coordination under reduced proprioceptive afference due to ischemia. The changes in movement coordination became evident as an increased coupling between arm postures during movement execution with final arm posture, resulting in a decreased number of kinematic DoF of the effector-system. Movement endpoint variability was not influenced by the ischemia. Thus, the canonical correlation analysis revealed that healthy participants are able to immediately and efficiently adjust their control strategy to the impaired flow of proprioceptive information. In conclusion, canonical correlation analysis provides a valuable method to advance our understanding of human movement control by offering an approach to analyze the temporal propagation of movement variability during discrete movements executed by multiple DoF-effector systems, such as the arm.

## AUTHOR CONTRIBUTIONS

All authors of this study contributed substantially to its conception, the interpretation of data, as well as drafting and revising earlier versions of the manuscript. In addition, MK and TE were responsible for data acquisition and analysis. MK, AS, and TE all approved the version to be published and agreed to be accountable for all aspects of the work.

#### FUNDING

This work was financially supported by the Research Training Group (GRK) 1091 "Orientation and motion in space" of the

#### REFERENCES

German Research Foundation (DFG) as well as by the Robert Bosch Foundation (Grant No. 32.5.G412.0007.0).

### ACKNOWLEDGMENTS

A modified version of this manuscript has been included as part of the dissertation thesis of MK (Krüger, 2013).

**Conflict of Interest Statement:** 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.

Copyright © 2017 Krüger, Straube and Eggert. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Predictive Simulation of Reaching Moving Targets Using Nonlinear Model Predictive Control

Naser Mehrabi\*, Reza Sharif Razavian, Borna Ghannadi and John McPhee

Systems Design Engineering, University of Waterloo, Waterloo, ON, Canada

This article investigates the application of optimal feedback control to trajectory planning in voluntary human arm movements. A nonlinear model predictive controller (NMPC) with a finite prediction horizon was used as the optimal feedback controller to predict the hand trajectory planning and execution of planar reaching tasks. The NMPC is completely predictive, and motion tracking or electromyography data are not required to obtain the limb trajectories. To present this concept, a two degree of freedom musculoskeletal planar arm model actuated by three pairs of antagonist muscles was used to simulate the human arm dynamics. This study is based on the assumption that the nervous system minimizes the muscular effort during goal-directed movements. The effects of prediction horizon length on the trajectory, velocity profile, and muscle activities of a reaching task are presented. The NMPC predictions of the hand trajectory to reach fixed and moving targets are in good agreement with the trajectories found by dynamic optimization and those from experiments. However, the hand velocity and muscle activations predicted by NMPC did not agree as well with experiments or with those found from dynamic optimization.

#### Edited by:

Massimo Sartori, University of Göttingen, Germany

#### Reviewed by:

Leonardo Abdala Elias, University of Campinas, Brazil Ton Van Den Bogert, Cleveland State University, USA

#### \*Correspondence: Naser Mehrabi

<EMAIL_ADDRESS>
Received: 14 September 2016 Accepted: 20 December 2016 Published: 13 January 2017

#### Citation:

Mehrabi N, Sharif Razavian R, Ghannadi B and McPhee J (2017) Predictive Simulation of Reaching Moving Targets Using Nonlinear Model Predictive Control. Front. Comput. Neurosci. 10:143. doi: 10.3389/fncom.2016.00143 Keywords: reaching, NMPC, prediction horizon, motor control

## INTRODUCTION

The human central nervous system (CNS), consisting of brain, and spinal cord, is responsible for controlling and maintaining body motions. As first formulated by Bernstein (1967), the CNS simultaneously coordinates the kinematics and kinetics of body motions, despite uncertain/unknown (future) trajectories and the redundancy in muscle actuators. As an example, in a goal-directed planar reaching task, where only the final position of the hand is specified, an infinite solution set of hand trajectories and muscle activation patterns exist to reach the final position. The early observations of reaching and pointing tasks led to the well-known "Minimum-X" models [e.g., minimum-jerk model (Flash and Hogan, 1985; Wada et al., 2001), minimum-torque-change model (Uno et al., 1989), minimum-variance model (Harris and Wolpert, 1998), and minimumwork model (Soechting et al., 1995)] to predict the hand trajectory. These models hypothesize that the CNS coordinates the body movement such that an exertion (X) is minimized. Later, this hypothesis is extended to consider physiologically-motivated exertions such as muscle activation effort (Crowninshield and Brand, 1981; Happee and Van der Helm, 1995; Ackermann and van den Bogert, 2010), metabolic energy expenditure (Anderson and Pandy, 2001; Peasgood et al., 2006), and muscle fatigue (Sharif Razavian and McPhee, 2015).

In computer simulations, the Minimum-X model has been successfully implemented using dynamic optimization (DO) to predict the average human motion for a given task. A common DO approach parameterizes the muscle activation profiles for the period of motion and searches the feasible space to find the profiles that minimize X (Davy and Audu, 1987; Yamaguchi and Zajac, 1990; Neptune and Hull, 1998; Anderson and Pandy, 2001; Kaplan and Heegaard, 2001; Sha and Thomas, 2013; Kistemaker et al., 2014). This approach provides an open-loop (feedforward) command of muscle activations to control the given task. This command can represent the descending command of a wellrepeated/well-learned task [e.g., platform diving (Koschorreck and Mombaur, 2011)]. In this model, the CNS only recalls the learned information, and does not intelligently adjust the commands in real-time. However, during conscious voluntary movements, the CNS has to continuously update the motor commands to correct for errors (Todorov, 2004). For example, previous studies (Sarlegna and Pratik, 2015) on pointing and reaching have shown that the CNS constantly updates the hand trajectory based on sensory (feedback) information. This sensory information can be received from vision, proprioception, audition, the vestibular system, and internal models that can predict the motion (Desmurget and Grafton, 2000).

Dynamic optimization implementation of minimum-X models raises an interesting question: does the CNS predict the trajectory at the beginning of the motion? Or does it constantly readjust the trajectory? If the latter is true, how far in advance does the CNS predict the motion, and how does that affect the motion? This article focuses on these questions and provides a computational platform to study the effects of the prediction horizon using optimal feedback control theory. Optimal control methods have been previously used to find a unique solution for motor coordination (Meyer et al., 1988; Loeb et al., 1990; Sporns and Edelman, 1993; Kuo, 1995; Anderson and Pandy, 2001; Todorov and Jordan, 2002b; Liu and Todorov, 2007); however, there are few applications of optimal feedback control to a nonlinear redundantly-actuated musculoskeletal model. The LQR (linear quadratic regulator) and LQG (linear quadratic Gaussian) control methods have been applied to a linear arm model to describe the hand trajectory (Harris and Wolpert, 1998; Todorov and Jordan, 2002a; Liu and Todorov, 2007). Later, to control the nonlinear dynamics of the neuromuscular system, an iterative LQG (iLQG) controller has been developed, in which the nonlinear model is iteratively linearized (Todorov and Li, 2005).

In the present research, a nonlinear model predictive controller (NMPC) with a finite prediction horizon is employed. Predicting infinitely into the future is highly improbable in humans, and a finite prediction horizon allows more realistic simulation. The NMPC allows us to consider the complexity and nonlinearity of the musculoskeletal system without compromising the accuracy and optimality, as occurs using model linearization. It can be formulated to simulate trajectory tracking and goal-oriented tasks with both fixed and moving targets where it only corrects the deviations from the task goal (Todorov and Jordan, 2002a). The NMPC is a simultaneous control method because the optimal trajectory and its required muscular activities are calculated at the same time. To the best of the authors' knowledge, this work is the first use of NMPC for fully predictive simulation of human reaching tasks.

In this research, we are not focused on the source of the sensory information; we assume that the current biomechanical states (posture and velocities) are available to the CNS when necessary. This assumption seems to be valid for healthy individuals, as a wide range of sensory organs is available to sense and transmit information to the CNS. However, a pathological condition might limit the CNS access to this available sensory information. For instance, in a deafferented patient, the sense of position (and therefore the motor skills) is largely lost due to the loss of somatosensory inputs (Bringoux et al., 2016).

This paper is organized as follows. In the Method section, the experimental procedure, the planar arm model, the nonlinear model predictive controller, and forward dynamic simulation framework are provided. Next, in the Results and Discussion section, the use of NMPC as the motor control unit in human reaching tasks is presented and discussed. This study investigates the use of anticipatory planning with continuous error correction by the CNS during reaching tasks. The first goal of this study is to study the effects of varying the prediction horizon on the hand trajectory and muscle activities in a reaching task. Therefore, we ran a number of NMPC simulations with various prediction horizons, as well as a DO simulation to obtain a "gold standard" for comparison. Secondly, the capability of the NMPC as an optimal feedback controller for tracking predefined trajectories has been investigated. This ability is useful when an expected/desired trajectory is available. Lastly, the effectiveness of the proposed NMPC for the simulation of reaching to moving targets is studied. We hypothesize that the anticipatory behavior of the CNS can be modeled by NMPC and verified by comparing the hand trajectory predicted by NMPC to those collected in experiments. Finally, Conclusions and Future Work are presented.

## METHODS

## Experiments

To examine the accuracy of the NMPC predictions, a 27 year old male subject was selected to perform reaching tasks. An Optotrak Certus motion capture system (Northern Digital Inc., NDI) was used to measure the arm trajectory at 30 Hz. In these experiments, the subject was seated with the arm elevated at the shoulder level. An active marker attached to the back of the hand (as shown in **Figure 1A**) has been used to capture centerout hand motion trajectories. The subject was asked to move his hand from an initial central position to one of eight final targets spread evenly on a circle of 20 cm radius at a self-selected convenient speed. The experiment was repeated 10 times for each target with 2 min rest intervals between each set. The subject also performed reaching to moving targets. He was instructed to reach to a target, which was relocated to another position midway through the movement. The subject was instructed to adjust his motion to reach to the moving target. The subject was also given 5 min to rest before performing the reaching moving targets experiment.

During the reaching trails, the electromyography (EMG) activity of seven muscles (anterior/middle/posterior deltoid, long/lateral triceps brachii, biceps brachii, and brachioradialis) were collected at 2000 Hz using a Trigno portable EMG system (Delsys Inc.). The EMGs were band-pass filtered (5–800 Hz cutoff), rectified, low-pass filtered (4 Hz cut-off), and normalized to maximum voluntary contractions (MVCs). We performed a Pearson correlation analysis to investigate the correlation between the EMGs and muscle activity predicted by NMPC and DO simulations. We resampled both the captured EMGs and simulation results with a sampling rate of 100 Hz, and performed the Pearson correlation analysis using the "Corr" command in MATLAB. The experiments have been approved by the Office of Research Ethics at the University of Waterloo and carried out with written informed consent from the subject. The subject gave written informed consent in accordance with the Declaration of Helsinki.

#### Planar Arm Model

In this research, a planar arm model similar to the one developed by Ghannadi et al. (2015) was used. The model consisted of torso, upper arm, and forearm to simulate the hand motion. The torso was fixed and the shoulder and elbow were modeled using revolute joints. Six muscle groups including shoulder and elbow mono-articular flexors/extensors and two bi-articular flexors/extensors were used to actuate the arm as shown in **Figure 1B**. A modified Hill-type muscle model with muscle excitation-to-activation dynamics was used to simulate the skeletal muscle contraction dynamics (see Appendix A for details). The muscle parameters of the planar arm model (i.e., insertion and origin positions, maximum isometric force, fiber optimal length, slack length, and pennation angle) were tuned to represent the dynamics of the upper extremity in the experimental condition (reaching targets in a horizontal plane elevated at the shoulder level). These parameters were tuned through a series of optimizations so that the planar model provides the same joint torques as a high-fidelity three-dimensional upper extremity model (Ghannadi et al., 2015). Kinematic and dynamic parameters of the arm model and the Hill muscle model parameters used here can be found in **Tables 1**, **2**, respectively.

## Principles of Nonlinear Model Predictive Control

The arm motion is controlled by complex commands descending from the CNS, which are the combination of the motion prediction (feedforward control) from an internal representation

\*About the center of mass, the mechanical y axis is assumed to be perpendicular to the plane of movement.

#### TABLE 2 | Hill-type muscle model parameters.

The details of how the upper-extremity muscle groups are lumped into these representative muscles can be found in Ghannadi et al. (2015).

of the body and environment (or so-called internal model; Desmurget and Grafton, 2000), and the corrective command from the sensory organs to correct any errors due to uncertainty or unknown environment (feedback control). This complexity is captured here by a model-based NMPC with a receding horizon. The NMPC uses a control-oriented model (COM) representing the human's internal model to predict the optimal trajectory, and feedback information to correct the prediction errors. The NMPC predicts the optimal dynamics of the system (x¯, u¯) over a prediction horizon as shown in **Figure 2A** by minimizing the following cost function:

$$J = \Psi \left( t\_0 + t\_{ph} \right) + \int\_{t\_0}^{t\_0 + t\_{ph}} \psi \left( \varkappa \left( t \right), \mu \left( t \right) \right) dt \tag{1}$$

$$\text{subject to: } 0 < \mu(t) < 1 \tag{2}$$

where 9 is the cost evaluated at the end of prediction horizon, ψ is the cost evaluated during the prediction horizon, and tph is the length of prediction horizon. The state variables at the current time (to) are obtained from the current sensory information. The input (u¯) is an optimal open-loop solution over the prediction horizon. If there are no external disturbances and no model uncertainty in the system, with infinitely long prediction horizon, the open-loop solution can be applied to the system for all time t > to. However, for the finite horizon case and in the presence of noise and uncertainty, the open-loop solution should only be applied until the next sampling time (t0 + δ). At the new time step, the optimal solution is re-evaluated with the new initial conditions for the receding horizon and iteratively applied to the system. By incorporating the feedback information, the NMPC is converted from a completely open-loop controller to an optimal closed-loop controller. The NMPC can handle constraints on both the states and the inputs. In musculoskeletal models, the muscle activation command must be non-negative and less than one, and constraints on states can be added to avoid unphysiological movements.

#### Mathematical Formulation of NMPC

In this article, the optimal dynamics over the prediction horizon were calculated using the GPOPS-II optimal control package that utilizes an orthogonal collocation method (Patterson and Rao, 2014). This method is a direct (simultaneous) optimization method in which both states (x) and inputs (u) are parameterized using a series of connected Legendre polynomials and become part of a Nonlinear Programing (NLP) problem. Here, the arm model described in the previous section was used as the COM of the NMPC and the simulation model for planar reaching and pointing tasks. The dynamic equation of the arm model can be described by:

$$\dot{\boldsymbol{x}} = \boldsymbol{f}\left(\boldsymbol{\varkappa}\left(\boldsymbol{t}\right), \boldsymbol{\mu}\left(\boldsymbol{t}\right)\right), \quad \boldsymbol{\varkappa}\left(\boldsymbol{0}\right) = \boldsymbol{\varkappa}\_{0} \tag{3}$$

where x ǫ R 10×1 are the arm model state variables consisting of shoulder angle and angular velocity plus elbow angle and angular velocity, and muscle activation states, and x0 is the vector of the initial states. The muscle excitation inputs u ǫ R 6×1 represent the ratio of excited motor units to the maximum number of motor units in that muscle.

In this research, for the particular case of a goal-directed reaching task, the terminal cost of the NMPC cost function (9) was removed, and the integral part (ψ) is computed from the summation of two terms: (i) a choice of specific physiological cost function, and (ii) a trajectory tracking error. Therefore, the NMPC cost function shown in (1) is converted to:

$$J = \int\_{t\_0}^{t\_0 + t\_{ph}} \left( p(\zeta(t) - \zeta\_{des})^2 + q \, G^M(u \,\left(t\right)) \right) \, dt \tag{4}$$

where p and q are cost function weightings, and ζ and ζdes are the hand position and its desired final value (in the Cartesian coordinate system), respectively. The simulated hand position ζ varies on the prediction horizon, while the desired final value

position trajectory. The hand moves from its original position to the marked position on its left.

ζdes was kept constant. The hand position is calculated from.

$$\begin{aligned} \xi &= \left[ L\_1 \cos(\theta\_1) + L\_2 \cos(\theta\_1 + \theta\_2) \right. \left. L\_1 \sin(\theta\_1) \right. \\ &+ L\_2 \sin(\theta\_1 + \theta\_2) \right]^T \end{aligned} \tag{5}$$

where θ1 and θ2 are shoulder and elbow angles, and L1 and L2 are upper arm and forearm lengths, respectively. The physiological cost G M (u (t)) is defined as:

$$G^{\mathcal{M}} = \mathfrak{u}^2 \tag{6}$$

The term G M in the cost function represents the neural excitation effort to perform the reaching tasks.

GPOPS-II finds the optimal dynamics of each given horizon by minimizing Equation (4) while satisfying inequality constraints related to muscle excitation (2) and equations of motion (3) using the Sparse Nonlinear Optimization (SNOPT) solver (Gill et al., 2005). An hp-adaptive mesh refinement method (Liu et al., 2015) has been used within GPOPS-II to refine the individual interval widths and the polynomial degree to reach a final optimal solution. Then, the first five-percent (e.g., 50 for 1000 ms prediction horizon) of optimal activations are applied to the muscles, and the arm motion is simulated. The new position and orientation of the arm are measured and sent back to the NMPC as initial conditions of the next iteration. In this research, we have assumed that the sensory organs can measure the exact joint angles and angular velocities. Uncertainty can be added to the measurements to account for the noise within the sensory organs, and to simulate the variability in the movement repetition. However, this has not been included in the scope of this work. The optimal muscle activations are shifted and considered as the initial guess of the next iteration.

#### Dynamic Optimization

In addition to the NMPC simulations, a dynamic optimization (DO) using GPOPS-II was performed to simulate the same task. Unlike the NMPC simulations, which continue until the position tracking error passes a certain threshold, the final simulation time and the final position of the hand are explicitly specified in the DO simulations. The DO cost function is:

$$J = \int\_0^{t\_f} G^M(\mu\_\cdot(\mathbf{t})) \, dt \tag{7}$$

subject to : 0 < u(t) < 1, and x tf = xf (8)

where tf is the final simulation time, and xf is the state vector corresponding to the target position. The same physiological cost function as in the NMPC simulations (Equation 6) was used in DO to compute the optimal hand position trajectory and muscle activations.

#### RESULTS

#### Effects of the Prediction Horizon Length

In this section, a goal-directed reaching task is simulated and the effect of prediction horizon length variation on the hand trajectory and muscle activation is studied. Here, the hand is initially at rest and in a natural position (θ1 = 44o and θ2 = 58o ) and moves toward a target 20 cm to the left of its initial position as shown in **Figure 2B**. This task was simulated using NMPC with 0.2, 0.3, 0.4, 0.5, and 0.8 s prediction horizons, and DO with a fixed time duration. In the NMPC simulations, the cost function weightings (p = 20 and q = 1) were kept the same. The DO final simulation time was chosen to be 1.5 s, in accordance with the experimental reach duration (1.431 ± 0.176 s in 10 repetitions).

**Figure 2B** demonstrates the hand trajectories for the aforementioned prediction horizons (tph). Despite the slight differences between the trajectories, the solutions with 0.4, 0.5, and 0.8 s prediction horizons closely correlate with the experimental results. The Pearson correlation factors between NMPC simulations and the experimental trajectories are 0.9351, 0.9758, 0.9947, 0.9991, 0.9994 respectively for 0.2, 0.3, 0.4, 0.5, 0.8 s. The Pearson factor between the DO and experiment is 0.867.

As shown in **Figure 3**, the elbow and shoulder angle variations with 0.2 s prediction horizon are less than those with longer prediction horizons. This signifies the importance of prediction horizon; with short prediction horizons, the controller is more cautious and takes longer to reach a desired position. Simultaneously, this results in smaller muscle activations for shorter horizons since the transient time to reach the final position is longer (see **Figure 4**). The reaching error at 1.5 s of the simulation is about 23% for 0.2 s prediction horizon, and reduces to 0.36% for 0.8 s prediction horizon. As expected for this motion, the shoulder mono-articular flexor (muscle 1), elbow mono-articular flexor (muscle 3) and bi-articular flexor (muscle 5) are activated at the beginning to accelerate the body; then, the antagonistic muscles are activated to reach a full stop at the desired position.

In the DO, similar to NMPC, flexor muscles are active at the beginning of the motion to accelerate the hand toward the target, then the extensor muscles are activated to stop the hand movement (a bang-bang control strategy). Since DO has to stop at the specified final time (1.5 s) the extensor muscle activities are larger than NMPC predictions at the decelerating phase of motion. On the other hand, the flexor muscle activities at the accelerating phase of the simulation are larger for NMPC than DO because the trajectory error at the beginning of the motion is large and exponentially reducing when it gets closer to target.

As shown in **Figure 4**, the muscle activations predicted by NMPC and DO simulations can capture the general trends of the experimental measurements. It can be observed that as the prediction horizon increases, the NMPC predicts larger muscle activities at the beginning of the motion. A Pearson correlation analysis was performed between the flexor muscle activities predicted by the NMPC and DO simulations and the EMGs from experimental measurements; the correlation coefficients are presented in **Table 3**. The correlation coefficient for extensor muscles are not reported since the EMG activity of these muscles were minimal in the experiments. As shown in **Table 3**, the correlation coefficient of flexor muscles reduces when the prediction horizon increases in the NMPC simulations, while the DO predictions correlate better with the experiments.

FIGURE 4 | Optimal activations of flexor and extensor muscles: (A) Optimal activation of shoulder mono-articular flexor (muscle 1), (B) Optimal activation of elbow mono-articular flexor (muscle 3), (C) Optimal activation of bi-articular flexor (muscle 5), (D) Optimal activation of shoulder mono-articular extensor (muscle 2), (E) Optimal activation of elbow mono-articular extensor (muscle 4), (F) Optimal activation of bi-articular extensor (muscle 6)

TABLE 3 | Pearson correlation analysis of flexor muscle activations predicted by NMPC and DO vs. experimental measurements.

As suggested by Morasso (1981), subjects tend to move in straight lines with bell-shaped tangential-velocity profiles when reaching for a target. **Figure 2B** shows that the NMPC with long-enough prediction horizon can realistically predict the hand trajectory, while **Figure 5** shows the differences in the velocity profiles of the NMPC and DO predictions. In these simulations, the NMPC tends to accelerate the hand more quickly than the DO; in DO, the optimization knows the final time and can distribute the acceleration over a longer time. On the contrary, in the receding horizon NMPC, the controller can only predict the motion as far as the prediction horizon. This results in the fast acceleration at the beginning the motion due to large tracking errors and slow deceleration at the end of motion due to small tracking errors.

## Reaching Targets Using a Predefined Trajectory

Reaching to eight different directions was also simulated using NMPC with 0.8 s prediction horizon. The target positions of the hand are located on a circle centered at the initial position of the hand with a radius of 20 cm (see **Figure 1**). In the NMPC simulations, a smooth 5th-order polynomial with zero initial and final velocities and accelerations is used as the desired straight-line hand trajectory. These trajectories begin at the initial position of the hand and end at the target positions. As shown in **Figure 6**, the NMPC is able to follow the desired trajectories, which are qualitatively correlated with the measured trajectories in experiments.

## Reaching a Moving Target

We have assumed that the CNS plans a trajectory to reach a target and constantly monitors the deviations from this trajectory and the target position. In this section, we study the case where the target position is suddenly relocated. Here, the hand is initially at rest at point O (**Figure 7A**) and moves toward the target at point A. Then 1 s later, the target position suddenly moves to the point B. This protocol was achieved in the lab by manually moving the target from its initial point A to B when the subject reached half the way to A. In this simulation, the time delay related to the visual cognition of this change [about 150 ms, (Jeannerod, 2006)] has not been considered.

**Figure 7** depicts that the NMPC controller can track the location of the target and correct the hand trajectory to reach the new target. As shown in **Figure 7A**, it seems that when the target moves, the subject over-compensates by moving the hand to the right, while the NMPC finds a trajectory that minimizes both position error and control effort. It is not possible to simulate this scenario with DO; it is one of its disadvantages compared to NMPC, which is able to make online adjustments to the hand trajectory.

## DISCUSSION

In this research, we presented a NMPC to mimic the human motor control system. The results showed that it can successfully replicate certain features of human motor control such as path planning and target tracking. This controller is a fully predictive

left of its initial position.

optimal controller that simultaneously solves the kinematic redundancy and muscle-sharing problem.

In hierarchical models (Menegaldo et al., 2006; Guigon et al., 2007; Mehrabi et al., 2015), the computational burden is reduced by separating the computation into two steps. In the first step, an optimal trajectory is generated based on a kinematic criterion; then, in the second step, the muscle sharing problem is solved based on another criterion (kinetic criterion). The main drawback of hierarchical models is that they follow a preplanned trajectory (output of the first step); therefore the online movement corrections are in favor of trajectory tracking. However, the NMPC controller simultaneously takes into account both kinetic and kinematic exertions to determine an optimal path along with the optimal muscle activations to achieve it. Therefore, it can be argued that this controller is more similar to the human CNS, as it receives proprioceptive information to adjust the predicted trajectory that satisfies the new condition and in favor of the end goal.

Thelen et al. (2003) developed a feedback/feedforward controller (computed muscle control, CMC) that uses inverse dynamics and static optimization to find muscle activities that track a set of desired kinematics. However, such an approach is applicable only if the kinematics is known or if there is a desired kinematics. One advantage of NMPC is its ability to control the motion with and without a prescribed motion, or when the target position moves.

In our study and for the first time, the effect of varying the prediction horizon on the path planning ability in reaching tasks has been investigated. As expected, increasing the prediction horizon improves the tracking performance, but makes the solution computationally more expensive. The prediction horizon length can be adjusted to capture the characteristic of a desired motion. Here, simulation results showed that the resultant hand trajectory with long enough prediction horizons resembles those found from the experiments. However, NMPC accelerates the hand faster and decelerates it slower than the bell-shaped speed trajectories reported by Morasso (1981) and observed in our experiments. This can be due to the fact that the reach time is not specified in the NMPC, or due to the selection of minimum control effort as the physiological cost function. The proposed NMPC is not limited to the suggested cost function; various cost functions can be implemented. For example, Kistemaker et al. (2014) studied the effect of different cost functions on the trajectory of the hand while performing a reaching task using a DO approach.

The Pearson correlation coefficients reported in **Table 3** show that muscle activities predicted by an NMPC with a short horizon can better predict (in a temporal sense) the experimental measurements. In contrast, the hand trajectory predictions (in a spatial sense) of an NMPC with larger prediction horizons can more closely replicate those from the experiments as shown in **Figure 2B**. Finally yet importantly, the NMPC simulations can be used to reproduce the experimental hand trajectories with a moving target as shown in **Figure 7**. The differences between the experiments and simulations may be due to the subject anticipation of another movement of the target point, while the final target position is known to the controller immediately following the shift. This unique feature of NMPC simulations can advance our theoretical understanding of hand movements and enables the next generation of assistive.

## Online Implementation of NMPC