# NEUROMECHANICS AND CONTROL OF PHYSICAL BEHAVIOR: FROM EXPERIMENTAL AND COMPUTATIONAL FORMULATIONS TO BIO-INSPIRED TECHNOLOGIES

EDITED BY : Massimo Sartori, Francisco J. Valero-Cuevas, Alfred C. Schouten, Matthew Tresch, Yoshihiko Nakamura and Manish Sreenivasa PUBLISHED IN : Frontiers in Computational Neuroscience and Frontiers in Neurorobotics

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-946-9 DOI 10.3389/978-2-88945-946-9

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# NEUROMECHANICS AND CONTROL OF PHYSICAL BEHAVIOR: FROM EXPERIMENTAL AND COMPUTATIONAL FORMULATIONS TO BIO-INSPIRED TECHNOLOGIES

Topic Editors:

Massimo Sartori, University of Twente, Netherlands Francisco J. Valero-Cuevas, University of Southern California, United States Alfred C. Schouten, Delft University of Technology, Netherlands Matthew Tresch, Northwestern University, United States Yoshihiko Nakamura, The University of Tokyo, Japan Manish Sreenivasa, University of Wollongong, Australia

Cover image: Modified from Jane\_Mori/Shutterstock.com

The term "neuromechanics" defines an integrative approach that combines the neuromuscular control and the biomechanical aspects of physical behavior in humans and animals. Crucial to this approach is a detailed description and modeling of the interaction between the nervous system and the controlled biomechanical plant. Only then do we have the broader context within which to understand evolution, movement mechanics, neural control, energetics, disability and rehabilitation.

In addition to enabling new basic science directions, understanding the interrelations between movement neural and mechanical function should also be leveraged for the development of personalized wearable technologies to augment or restore the motor capabilities of healthy or impaired individuals. Similarly, this understanding will empower us to revisit current approaches to the design and control of robotic and humanoid systems to produce truly versatile human-like physical behavior and adaptation in real-world environments. This Research Topic is therefore poised at an opportune moment to promote understanding of apparently disparate topics into a coherent focus.

Citation: Sartori, M., Valero-Cuevas, F. J., Schouten, A. C., Tresch, M., Nakamura, Y., Sreenivasa, M., eds. (2019). Neuromechanics and Control of Physical Behavior: From Experimental and Computational Formulations to Bio-inspired Technologies. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-946-9

# Table of Contents

## CHAPTER 1

#### INTRODUCTION

*06 Editorial: Neuromechanics and Control of Physical Behavior: From Experimental and Computational Formulations to Bio-inspired Technologies* Manish Sreenivasa, Francisco J. Valero-Cuevas, Matthew Tresch, Yoshihiko Nakamura, Alfred C. Schouten and Massimo Sartori

## CHAPTER 2

#### EXPERIMENTAL AND COMPUTATIONAL TECHNIQUES

*10 The Dynamics of Voluntary Force Production in Afferented Muscle Influence Involuntary Tremor*

Christopher M. Laine, Akira Nagamori and Francisco J. Valero-Cuevas

Keivan Mojtahedi, Bryan Whitsell, Panagiotis Artemiadis and Marco Santello

*46 Forearm Flexor Muscles in Children With Cerebral Palsy are Weak, Thin and Stiff*

Ferdinand von Walden, Kian Jalaleddini, Björn Evertsson, Johanna Friberg, Francisco J. Valero-Cuevas and Eva Pontén

## CHAPTER 3

#### MODELLING FORMULATIONS

*97 Predictive Simulation of Reaching Moving Targets Using Nonlinear Model Predictive Control*

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

*109 Evaluation of a Neuromechanical Walking Control Model Using Disturbance Experiments*

Seungmoon Song and Hartmut Geyer

## CHAPTER 4

#### BIO-INSPIRED TECHNOLOGIES

Alexei V. Alexandrov, Vittorio Lippi, Thomas Mergner, Alexander A. Frolov, Georg Hettich and Dusan Husek

Andrés Úbeda, José M. Azorín, Dario Farina and Massimo Sartori

# Editorial: Neuromechanics and Control of Physical Behavior: From Experimental and Computational Formulations to Bio-inspired Technologies

Manish Sreenivasa1 \*, Francisco J. Valero-Cuevas 2 , Matthew Tresch3,4 , Yoshihiko Nakamura5 , Alfred C. Schouten6,7 and Massimo Sartori 7 \*

*1 Department of Mechanical, Material, Mechatronics and Biomedical Engineering, University of Wollongong, Wollongong, NSW, Australia, 2 Division of Biokinesiology and Physical Therapy, Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, United States, 3 Shirley Ryan AbilityLab, Chicago, IL, United States, 4 Department of Biomedical Engineering, Northwestern University, Evanston, IL, United States, 5 Department of Mechano-Informatics, School of Information Science and Technology, University of Tokyo, Tokyo, Japan, 6 Department of Biomechanical Engineering, Delft University of Technology, Delft, Netherlands, 7 Department of Biomechanical Engineering, University of Twente, Enschede, Netherlands*

Keywords: sensory-motor system, biorobotics, neuromechanical modeling, motor control, neuro-technologies, neural control of movement

**Editorial on the Research Topic**

**Neuromechanics and Control of Physical Behavior: From Experimental and Computational Formulations to Bio-inspired Technologies**

## INTRODUCTION

#### The motivation behind this research topic was to cut across conventional boundaries that separate movement neuroscience, biomechanics, and robotics. The aim was to underscore that brain and body collaborate to produce behavior in biological organisms. While this is a simple idea, compartmentalization in education and science has often artificially separated brain from body. We also bring forward research paradigms to investigate physical behavior at the interface between humans and interacting robots. Understanding human-robot physical interaction requires the understanding of the complex interplay between brain, body, and the external environment. This could be achieved by employing a neuro-mechanical approach to the study of human and robot movement.

Within the context of a neuromechanical approach, we aimed to collate a research corpus that included work on experimental and neurophysiological analysis, computational modeling, and applications in rehabilitation and bio-inspired robotics. The 22 contributions to this research topic provide a wide range of perspectives and methodologies. The high-caliber contributions to this research topic also highlight the existence of a significant community of researchers interested in an interdisciplinary view toward the study of brain-body and human-robot interactions.

## THE CONTRIBUTIONS

Several of the studies enclosed in this research topic combined experimental work with state-of-the-art technologies and novel methodologies to measure and analyze neural control in humans (Banks et al.; Mojtahedi et al.; Logan et al.; Krüger et al.; Reyes et al.; Úbeda et al.).

Edited and reviewed by: *Si Wu,*

#### *Peking University, China* \*Correspondence:

*Manish Sreenivasa<EMAIL_ADDRESS>Massimo Sartori<EMAIL_ADDRESS>
Received: *23 January 2019* Accepted: *15 February 2019* Published: *19 March 2019*

#### Citation:

*Sreenivasa M, Valero-Cuevas FJ, Tresch M, Nakamura Y, Schouten AC and Sartori M (2019) Editorial: Neuromechanics and Control of Physical Behavior: From Experimental and Computational Formulations to Bio-inspired Technologies. Front. Comput. Neurosci. 13:13. doi: 10.3389/fncom.2019.00013* A common theme among many studies was a cohesive and complimentary use of numerical modeling, signal processing and experimental approaches, with one supporting and enhancing the findings of the other to understand the basic mechanisms of movement (Guarín and Kearney; Golkar et al.; Song and Geyer; Petric et al. ˇ ; Laine et al.; Mehrabi et al.; Sreenivasa et al.; Alexandrov et al.; Lippi and Mergner; Pasma et al.; Von Walden et al.). Yet another inspiring aspect of this research topic was the crossover toward efforts in biomarkers of pathology and rehabilitation (Laine et al.; Reyes et al.; Sreenivasa et al.; Von Walden et al.; Shuman et al.; Banks et al.; Úbeda et al.) and robot control (Pasma et al.; Lippi and Mergner; Alexandrov et al.; Szczecinski et al.). Taken together, the contributions to this research topic illustrate the wide scope of research being conducted in neuromechanics and the vital role this will play in the future of rehabilitation and robotics.

The remainder of this editorial provides an overview of the scientific contributions. These are presented within three macro categories including:

## NEUROMECHANICS

Although there is detailed understanding of the mechanisms taking place in the central nervous system or in the musculoskeletal system during movement, there is far less knowledge of the interplay between these two systems during complex motor tasks (Tresch and Jindrich, 2014; Sartori et al., 2016, 2017; Sreenivasa et al., 2016; Valero-Cuevas and Santello, 2017; Cohn et al., 2018; Schouten and Mugge, 2019). This motivates a neuromechanical approach to the study of human behavior when producing movements and forces to interact with the environment.

One way to better understand the interplay between the neural and biomechanical apparatus is by identifying the system dynamics using a combination of models and experimental observations. Guarín and Kearney developed a methodology to dissociate the time-varying intrinsic and reflex components during movements at the human ankle joint. Their results indicate that the joint stiffness is modulated dynamically during non-stationary movements. The relevance of their work lies in proposing methodologies for extracting joint stiffness profiles during dynamic conditions, a step forward with respect to current methods that predominantly investigate quasi-stationary scenarios. Golkar et al. also focused on dynamic stiffness about the ankle joint, using human experiments to study the modulation of intrinsic and reflex stiffness with muscle activation. Logan et al. applied harmonic transfer function analysis to investigate how humans control upright posture and speed during locomotion. Their results support the existence of a temporal hierarchy of subtasks during locomotion, with the control upper-body posture taking precedence over others. Schumacher and Seyfarth used a feedback-based neuromuscular modeling approach to study reflex gains during hopping. Interestingly, their results indicate that different combinations of gain values enable optimizing specific hopping characteristics (e.g., performance and efficiency), and that this was invariant with respect to changes in the model's mechanical properties. With another feedback-based approach, Song and Geyer investigated the response of a walking model to simulated disturbances such as electrical stimulation and perturbations applied to the legs and to the whole-body. Their results show that the model behavior compares favorably to experimental recordings of similar disturbances during human locomotion. In a conclusion that resonates with this research topic, the authors comment on how such model approaches can complement the experimental study of human motor control.

Another branch of contributions investigated whole-body human movements (either free or when interacting with the environment), which require the coordination of multiple degrees of freedom, and are often directed toward completing specific tasks accurately and efficiently. The study of such movements provides the means to answer fundamental questions about motor control and brain-body interaction. Krüger et al. applied canonical correlation analysis to study movement variability in discrete goal-directed reaching tasks. They tested the effects of ischemia, (temporary) artificially reduced blood flow to the arms, and found that while this reduced the complexity of movement control, the endpoint variability did not increase. Petric et al. ˇ comment on the lack of movement studies on real-world tasks such as hammering a nail into wood. In their work, they investigated the kinematics and dynamics of hammering both from an experimental and modeling perspective. They found that for periodic impact tasks, a model that relates the distance moved as a logarithmic function of time provided better predictions than Fitts' Law. Le Mouel and Brette proposed the hypothesis that postural adjustments follow spatial and temporal patterns that provide impetus for future movement. In a thought-provoking article, they relate this ability to learned proficiency in movements and the impairment of these skills with aging. Mojtahedi et al. studied how two humans collaborate to complete a physical task and the strategies used to infer movement direction. Interestingly, their findings suggest that the modulation of arm stiffness during cooperative tasks could be an effective means to communicate intended direction of movement. The authors relate this to future applications in human-robot interaction.

Even the simplest of movements often require the coordinated and smooth control of several muscles. Investigations at the neuromuscular level can provide an interesting way to reverseengineer how the brain controls movements. Reyes et al. studied corticomuscular and intermuscular coherence between the muscles of the finger and the thumb during pinch grips. They found significantly reduced coherence when individuated control of the thumb and index finger was required, and interpret their findings as supporting the notion that the cortex bonds task-related motor neurons into task-dependent functional units. Laine et al. used complimentary experimental and modeling approaches to study how the neuromechanics of voluntary force production in the index finger were sufficient to produce involuntary tremor. The authors suggest that their findings represent a new category of tremor, one related to the viscoelastic response of a closed-loop musculotendon system to dynamic inputs. Moreover, they underscore the importance of coherence analysis, a well-established and neurophysiologically informative phenomenon (Farmer, 1998), to the study of muscle coordination in neuromechanics. Mehrabi et al. modeled a musculoskeletal planar arm and compared predictions using non-linear model predictive control (NMPC), to those from dynamic optimization and to experimental recordings. The authors motivate the NMPC approach by highlighting its effectiveness in simulating motion where the kinematics may or may not be prescribed, and where the target position may move during execution. Their results show that NMPC predicted the hand trajectory quite well, but not the hand velocities or muscle activations. Sreenivasa et al. developed optimal control-based predictive simulations of a child's gait with and without an anklefoot orthosis. The motivation behind their work was to provide modeling-based tools to identify the patient-specific optimal stiffness of the orthosis, while reducing the effort required for the clinical procedure of fitting an orthosis to a patient. These studies inform the debate on the extent to which optimality (as understood by engineers) occurs in the neural control of our bodies (Loeb, 2012).

## BIOTECHNOLOGY AND REHABILITATION

The knowledge gained from the study of human motor control, and the methods developed to analyze and synthesize human biomechanics are invaluable for the development of biomarkers of pathology, biotechnologies and rehabilitation processes that directly interact with the human body during complex movements. In a review article, Pizzolato et al. set an optimistic yet cautious note on the application of computational musculoskeletal models for the design and control of wearable devices. An interesting focus of the review was on bioinspired technologies that not only support or augment movement, but also foster tissue strengthening and repair through optimal tissue loading. Von Walden et al. studied the forearm flexor muscles in children with cerebral palsy (CP). They experimentally record the muscle strength and cross-sectional area, and use a modeling approach to estimate passive muscle stiffness independent of reflex activity. With a study title that conveys a concise takehome message, the authors found that the forearm flexor muscles of children with CP are weak, thin, and stiff, something that may inform personalized treatments. This underscores the plastic relationship between brain and body: damage to the neural controller can lead to physical changes in the anatomical actuators and plant.

The conclusions from Shuman et al. and Banks et al. are of significant importance to the research on muscle synergies during movements. Shuman et al. comment on the impact that filtering and scaling of recorded electromyography (EMG) data can have on the computed synergies for CP and typicallydeveloping children, and the implications this can have when comparing results across studies and laboratories. Banks et al. investigated 30 variations in the methodology used for muscle synergy analysis, and evaluated the impact of these variations in identifying responders in post-stroke population. Their results highlight the sensitivity of MSA to methodological choices, and a need to standardize and/or provide exacting detail about the methodology used in future work. Úbeda et al. proposed a method to decode spinal primitives of multi-muscle control from electroenchephalography (EEG) recordings. The authors apply their method to healthy individuals as well as to patients with incomplete spinal cord injury, and comment on the future usage toward a new class of brain-exoskeleton interfaces. This work, together with that of Reyes et al. and Laine et al. helps clarify the ongoing debate on whether and how one can detect prescriptive synergies of neural origin—as opposed to descriptive synergies of numerical origin (Brock and Valero-Cuevas, 2016).

## BIO-INSPIRED ROBOTICS

Robot sensors and actuators remain vastly different from those on humans and animals, yet bio-inspired control methods find rich application in robotics. Szczecinski et al. developed methodology to design subnetworks that perform specific mathematical operations, and can be assembled into larger networks to mimic to some extent an animal's nervous system. Humanoid robots, with their anthropomorphic form and capabilities for bipedalism, are an excellent example for the translation of insights from human movement research to robotics. One interesting focus is on maintaining balance and posture, with robots possibly learning from humans, but also providing a real-world platform to test models and hypotheses for human balance. Alexandrov et al. proposed the use of the eigen-movement control concept derived from human hip-ankle coordination, as a simpler and more stable way to control humanoid robot balance. Lippi and Mergner implemented a modular human-derived control architecture for maintaining balance and posture in a humanoid robot. Interestingly, their results show that the mechanical coupling from the robot's body was sufficient to stabilize movements in the frontal and sagittal planes, without a need for an explicit link between the respective control modules. Pasma et al. quantified human balance control using the independent channel (IC) model and implemented the IC model in a humanoid robot. Their results show that the IC model, a descriptive model in the frequency domain, can imitate human balance behavior in real world situations with a humanoid robot. This provides further evidence that the IC model is a valid description of human balance control.

## CONCLUSIONS

This research topic posed direct focus on the neuromechanics of movements and the forces generated in interaction with the environment. This is described as an integrative approach that combines the neuromuscular control and the biomechanical aspects of body and the physics of the tasks in humans. We assert that this approach is overdue and necessary to obtain the theoretical and experimental frameworks for understanding the evolution, mechanics, neural control, energetics, disability, and rehabilitation of physical behavior. Moreover, understanding the neuromechanics of physical behavior should also leverage the development of personalized wearable robotic technologies that can interact with biological tissues within the composite neuromuscular system. This is central to ultimately mimic, restore or augment motor capabilities in healthy or impaired individuals.

## AUTHOR CONTRIBUTIONS

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

## REFERENCES

## FUNDING

Research reported in this publication was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the National Institutes of Health under Awards Number R01 AR- 050520 and R01 AR-052345 to FV-C. This work was also supported by the Department of Defense CDMRP Grant MR150091 and Award W911NF1820264 from the DARPA-L2M program to FV-C. MCT acknowledges funding from the NIH under the grant number NS086973. This work was partly supported by the European Research Council (ERC) Starting Grant INTERACT (No. 803035) awarded to M. Sartori, under the European Union's Horizon 2020 research and innovation programme.

**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 © 2019 Sreenivasa, Valero-Cuevas, Tresch, Nakamura, Schouten and Sartori. 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) and the copyright owner(s) 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 Dynamics of Voluntary Force Production in Afferented Muscle Influence Involuntary Tremor

#### Christopher M. Laine1, 2, Akira Nagamori 2 and Francisco J. Valero-Cuevas 1, 2 \*

1 Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, USA, 2 Division of Biokinesiology and Physical Therapy, University of Southern California, Los Angeles, CA, USA

Voluntary control of force is always marked by some degree of error and unsteadiness. Both neural and mechanical factors contribute to these fluctuations, but how they interact to produce them is poorly understood. In this study, we identify and characterize a previously undescribed neuromechanical interaction where the dynamics of voluntary force production suffice to generate involuntary tremor. Specifically, participants were asked to produce isometric force with the index finger and use visual feedback to track a sinusoidal target spanning 5–9% of each individual's maximal voluntary force level. Force fluctuations and EMG activity over the flexor digitorum superficialis (FDS) muscle were recorded and their frequency content was analyzed as a function of target phase. Force variability in either the 1–5 or 6–15 Hz frequency ranges tended to be largest at the peaks and valleys of the target sinusoid. In those same periods, FDS EMG activity was synchronized with force fluctuations. We then constructed a physiologically-realistic computer simulation in which a muscle-tendon complex was set inside of a feedback-driven control loop. Surprisingly, the model sufficed to produce phase-dependent modulation of tremor similar to that observed in humans. Further, the gain of afferent feedback from muscle spindles was critical for appropriately amplifying and shaping this tremor. We suggest that the experimentally-induced tremor may represent the response of a viscoelastic muscle-tendon system to dynamic drive, and therefore does not fall into known categories of tremor generation, such as tremorogenic descending drive, stretch-reflex loop oscillations, motor unit behavior, or mechanical resonance. Our findings motivate future efforts to understand tremor from a perspective that considers neuromechanical coupling within the context of closed-loop control. The strategy of combining experimental recordings with physiologically-sound simulations will enable thorough exploration of neural and mechanical contributions to force control in health and disease.

Keywords: tremor, force control, muscle models, closed-loop system, dynamics

## INTRODUCTION

It is well known that humans cannot produce a perfectly stable force. Within the context of precise, goal-directed actions, involuntary force fluctuations can reveal clinically relevant information about neuromuscular control in disorders such as dystonia (Xia and Bush, 2007; Chu and Sanger, 2009), Parkinson's disease (Vaillancourt et al., 2001; Ko et al., 2015), bruxism (Laine et al., 2015b),

#### Edited by:

Sliman J. Bensmaia, University of Chicago, USA

#### Reviewed by:

Andrew Joseph Fuglevand, University of Arizona, USA Juan Álvaro Gallego, Northwestern University, USA Helen Huang, North Carolina State University and University of North Carolina at Chapel Hill, USA

\*Correspondence:

Francisco J. Valero-Cuevas<EMAIL_ADDRESS>
Received: 29 April 2016 Accepted: 02 August 2016 Published: 19 August 2016

#### Citation:

Laine CM, Nagamori A and Valero-Cuevas FJ (2016) The Dynamics of Voluntary Force Production in Afferented Muscle Influence Involuntary Tremor. Front. Comput. Neurosci. 10:86. doi: 10.3389/fncom.2016.00086 and essential tremor (Héroux et al., 2010), among others. In such tasks, the nature of force variability may be influenced by both central and peripheral components of sensorimotor integration.

Unfortunately, the utility of measuring involuntary force fluctuations (i.e., tremor) within scientific or clinical settings has been limited due to the large and often ambiguous set of factors which can influence such measures. In some cases, tremor may reflect a mechanical resonance whose frequency depends on the physical characteristics of the muscle/limb in question (Lakie et al., 2012; Vernooij et al., 2013). At the same time, tremor may stem from cycles of excitation around the stretch-reflex loop (Lippold, 1970; Young and Hagbarth, 1980; Christakos et al., 2006; Erimaki and Christakos, 2008). The two mechanisms likely interact, since reflex activity is itself influenced by muscle/tendon compliance (Rack et al., 1983), limb loading (Joyce and Rack, 1974), contraction history (Gregory et al., 1998), and the temporal dynamics of force production (Xia et al., 2005).

It is clear that the specific type and extent of neuromechanical coupling influencing performance of a given task are of key importance for understanding the generation of force variability. Understanding the factors which influence dynamic force control is especially important given that this is the basis of manual dexterity during activities of daily living. However, the neural and/or mechanical origins of unintended force variability are not always clear, particularly within the context of dynamic force control.

In this study, we investigated the relationship between voluntary force production and involuntary force variability in a group of healthy adults engaged in a dynamic, isometric force tracking task. Given the various links between reflex activity, contraction dynamics, and tremor, our hypothesis was that involuntary force variability would depend upon voluntary contraction dynamics. In order to better understand the potential sources of force variability within our experimental task, we used a physiologically-realistic computer simulation to determine the sufficiency of muscle-tendon mechanics and reflex pathways to reproduce our experimental results. The simulation also allowed us to characterize the sensitivity of force variability to parameters such as reflex gain.

The significance of our study is two-fold. First, we describe a novel source of tremor along with a method for its experimental induction, and strong evidence for its origin in musculotendon dynamics. Second, the sensitivity of this tremor to both neural and mechanical factors within our simulation implies that simple force tracking tasks, such as described here, may represent a novel approach to investigating peripheral components of sensorimotor integration in health and disease.

#### METHODS

All procedures were approved by the institutional review board at the University of Southern California and all participants gave informed written consent prior to participation. Ten healthy participants were recruited (4 female, 6 male, aged 23–31 years) to carry out force tracking experiments.

## Physiological Data

#### Task

Participants were seated ∼1 m from a 17-inch computer monitor which displayed a sinusoidal target with a vertical range representing forces from 5 to 9% of the maximum force that each individual could exert with the index finger of their self-reported dominant hand (see **Figure 1A**). Visual feedback of exerted force was provided in the form of a cursor which moved left to right across the computer screen for 40 s before looping back to the left. Prior to recordings, participants practiced tracking several target cycles to become familiar with the task. Each participant then tracked the 0.25 Hz sinusoidal target for two 80 s trials separated by several minutes of rest. A slow sinusoidal target is a rich behavior that is ideal for probing dynamic dependencies. For example, if tremor depended on force velocity, then tremor amplitudes would appear to follow the derivative of the target sinusoid (i.e., a cosine). If one direction of force (increasing vs. decreasing) were tied to tremor amplitudes, then tremor amplitudes would be largest along either the rising or falling phase of the target sinusoid. If the magnitude of force were most relevant, tremor amplitudes would essentially follow the target trajectory, being largest at the peaks and smallest at the valleys.

As depicted in **Figure 1A**, a miniature single-axis force transducer was fixed to the top of a plastic cylinder and located under the tip of the finger. Participants were asked to produce a downward force perpendicular to the force sensor, an action requiring contraction of the index finger slip of the flexor digitorum superficialis (FDS) muscle. This particular muscle and joint action were chosen because flexion at the proximal interphalangeal (PIP) joint is necessary for manipulation activities of daily living, and because this straightforward mechanical action is well suited for simulation.

#### Force and EMG Measurements

Surface EMG recordings were made over the distal portion of the index finger slip of the FDS muscle using an active bipolar electrode (Biometrics Ltd, Newport, UK) grounded at the wrist. Confirmation of correct electrode positioning was accomplished via palpation of the distal muscle belly (∼7 cm proximal to the crease of the wrist, on the ulnar side) during index finger flexion, as well as observation of ongoing EMG signals during PIP joint flexion/extension and during our isometric task. The EMG signals were acquired at 1000 samples per second using a Biometrics DataLog system and associated software. The measurement of force, and the display of visual feedback to participants, was accomplished using custom MATLAB (The MathWorks, Natick, MA, USA) scripts to acquire force signals from a miniature load cell (ELB4-10, Measurement Specialties, Hampton, VA, USA) using a USB-DAQ (National Instruments, Austin TX, USA). The data acquisition unit sent a synchronization pulse to the biometrics system at the start of each recording. The data were analyzed offline using custom MATLAB scripts.

## Data Analysis

#### Conversion of Sinusoidal Force to Instantaneous Phase

To uncover the slow, voluntary force associated with the intended target trajectory, the force produced by each participant was low-pass filtered at 0.5 Hz. Using the Hilbert transform, the instantaneous phase of this tracking force was calculated and expressed in degrees (0–360◦ ) over the course of each target cycle. This conversion was useful since instantaneous phase holds information about the actual dynamics of force production at a given time, regardless of tracking error. Although tracking error was not a focus of this investigation, it was still important to eliminate poorly tracked target cycles. Target cycles in which the absolute tracking error exceeded 4% of a participant's MVC level at any time point were excluded from all further analysis.

#### Calculation of Instantaneous Tremor Amplitude

To quantify the presence of involuntary force fluctuations, the sinusoidal force trajectories produced by each participant were filtered into two different frequency bands.

First, we investigated the presence and magnitude of force fluctuations at high frequencies (>6 Hz), which cover the frequency range of physiological tremor (Lippold, 1970; Elble and Randall, 1976; Burne et al., 1984; Christakos et al., 2006). These force fluctuations were extracted by band-pass-filtering the force produced by each individual between 6 and 15 Hz (zero-phase, 4th order Butterworth filter). Beyond about 15 Hz, the amplitude of force fluctuations is essentially negligible due to the low-pass filtering effects of tissue (finger pad, tendon) and muscle.

Second, we quantified slow (1–5 Hz) force fluctuations. Generally, slow fluctuations in force stem from changes in the overall drive to motor neurons (Allum et al., 1978; De Luca et al., 1982; Miall et al., 1993; Slifkin et al., 2000; Squeri et al., 2010). These slow fluctuations do include voluntary correction of tracking errors, but our analysis focused on fluctuations that were consistently present at particular phases of the target cycle, and therefore reflect an involuntary process. To extract these fluctuations, we used a 1–5 Hz band-pass-filter (zero phase, 4th order Butterworth filter).

An example of band-pass filtered force traces in relation to the target sinusoid is depicted in **Figure 1A** (right).

#### Calculation of Tremor Modulation as a Function of Tracking Phase

As described above, the tracking force produced over time by each participant was converted to a trace of instantaneous phase angles where each complete target cycle was represented as a progression from 0 to 360◦ . Each cycle was then divided into 36 phase-bins (each representing 10◦ ). Again, it should be noted that all analyses are based on the temporal dynamics of the force produced by the participants and not on the displayed target. This renders any positive or negative tracking lags irrelevant [although they would be minimal given the highly feed-forward nature of this type of task (Erimaki et al., 2013)]. To examine the relationship between tracking phase and force variability, we first converted each band-pass filtered force signal into an instantaneous amplitude signal by rectification

FIGURE 1 | Experimental task and simulation. Participants were asked to produce isometric force against a small load cell with the tip of the index finger (A). Using feedback of their applied force, participants tracked a 0.25 Hz sinusoid during 80 s trials. The range of the target sinusoid was 5–9% of each individual's maximum voluntary contraction (MVC) force. An example of the recorded force from one participant is depicted to the right (top trace). Below, the same force trace has been filtered into 2 different bands, a 1–5 Hz band which captures slow tremor and tracking error (middle), and a 6–15 Hz band which reveals physiological tremor and fast twitches (bottom trace). In this example, both high and low frequency force fluctuations appear to depend upon the phase of the tracked sinusoid. (B) Depicts the implementation of a control loop used to simulate the tracking experiment above. At the "spinal cord," three sources of input are summated. The first two sources are proprioceptive signals from the muscle spindle (via group Ia and II afferents), and Golgi tendon organ (GTO) (via Ib afferents) models, which provide positive and negative feedback to the "spinal cord," respectively. The third input is a tracking signal sent from a supraspinal controller. This "tracking controller" sends an output (C) which is continually updated according to the difference between the target and the force (F) produced by the muscle (see text for details). From the "spinal cord," a neural drive (ND) signal is sent to a Hill-type muscle model via an "activation filter" which further shapes the neural drive to account for calcium dynamics within physiological muscle. The filtered muscle activation (MA) signal sent from the "activation filter" generates contraction (Fm), accounting for physical properties of the muscle-tendon complex, such as a series and parallel elastic element (SE and PE), mass (M), viscosity (B), and pennation angle (alpha). The delays and sign associated with each feedback loop are also depicted.

and smoothing with a 200 ms Gaussian window. The magnitude of the resulting smoothed/rectified signal also serves as a simple estimation of instantaneous variance within the specified frequency band, given the equivalence between total signal power (in frequency domain) and total signal variance (in time domain) (i.e., Parseval's theorem). We then found the sum of the filtered/rectified force values associated with each phase interval, and divided each by the integral of the filtered/rectified force trace. This procedure gives the relative proportion of total force variability (within the specified frequency band) associated with each 10◦ phase interval of the target cycle.

Under the null hypothesis, the proportion of force variability in each phase bin does not depend on the phase progression of the target cycle. Thus each 10◦ phase bin would be expected to show about 2.8% (100%/36 bins) of the total force variance. To test the null hypothesis, we compared our recorded proportions to a phase-randomized distribution generated directly from the recorded data (i.e., shuffled versions of our own data). We constructed these shuffled distributions of proportion values by randomly selecting a different 10◦ phase bin in each tracked target cycle, and then calculating the proportion of total force variability, as previously described. The process was repeated 5000 times, creating a distribution of shuffled proportion values which allowed us to determine a 95% confidence interval. Proportions falling outside of this interval would then represent statistically significant deviations from chance level. Our use of a Monte–Carlo method provides a direct, conservative, and assumption-free statistical analysis. Similar methods are often used in neuroscience, where analysis of real vs. shuffled/randomized neural activity is common (Perkel et al., 1967; Tam et al., 1988; Türker et al., 1996; Rivlin-Etzion et al., 2006; Laine et al., 2012). In our case, alternative methods such as testing for differences between individual phase bins, would be ill-suited for identifying the timing of tremor modulation with respect to the target phase, and would also not account for the fact that the proportion within each phase bin is not strictly an independent measurement.

The above methodology was applied to individual participants. To evaluate the population as a whole, the proportions for each phase bin were averaged across individuals. As a statistical evaluation, we calculated, for each phase bin, the number of individuals whose tremor proportion fell above or below the 95% confidence interval. For any given phase bin, a 5% error rate might be expected. Since our analysis included 10 individuals, it could be expected that at least one may have exceeded the confidence level purely by chance. However, the binomial probability that 2 of 10 individuals should show a (false positive) significance at the 95% confidence level is 0.015. For this reason, our population significance level was set to 0.015, or 2 out of 10 individuals, for our consistency analysis.

In addition to analyzing the proportion of force variance in each phase interval, we also calculated the cross-cycle average tremor amplitude in each phase interval. This analysis yielded an amplitude profile for each individual (and frequency band), similar to the proportion profiles described above. We then recorded the maximum and minimum values observed in the amplitude profile of each individual, regardless of the particular phase at which these values were found. This allowed us to evaluate the actual extent of tremor amplitude modulation, uncoupled from any particular pattern of tremor modulation across target phases.

#### Force to EMG Coherence across Target Phases

Coherence is a frequency-domain measure of synchronization (primarily phase-locking) between signals, and is bounded between 0 (no correlation between signals) and 1 (perfect linear correlation). Coherence between rectified EMG activity and force is useful for identifying the frequency content of force-relevant neural drive to muscles, since action potential shapes/sizes and other recording artifacts only influence the EMG spectrum, but would not be synchronized with force. In addition, coherence between FDS activity and force tremor provides validation that our simulation of a dynamically activated FDS muscle is appropriate for exploring the potential origins of recorded force fluctuations.

To calculate EMG to force coherence, the force and EMG signals were concatenated across all trials from all subjects to form two long signals. These signals were then converted to timefrequency-representations (TFRs) via wavelet analysis. We chose a wavelet approach so that we could precisely determine which frequencies of force were synchronized with EMG, and at what times. The technique is common where temporal variation of spectral power or synchronization is of interest (e.g., Siemionow et al., 2010; Tscharner et al., 2011). This was accomplished through convolution of each original signal x(t) with a Gaussianwindowed complex sinusoid (a Morlet wavelet), the duration of which was set to span 3 cycles of each frequency (f) from 1 to 20 Hz. The process can be expressed by the following formula:

$$TFR\left(t,f\right) = \int \propto \left(t\right) \frac{1}{\sigma\sqrt{2\pi}}e^{-\left(t-\tau\right)^2}e^{-j2\pi f\left(t-\tau\right)}dt$$

where the standard deviation (σ) of the Gaussian window is set to 3/(2πf). The force trace (band-pass filtered between 1 and 20 Hz) as well as the EMG activity (rectified, normalized per subject to have unit variance) were thus converted to complex-valued TFRs (herein defined as TFR\_Force and TFR\_EMG, respectively). The spectral power of each signal can be calculated as follows:

Power\_Force (t, f) = TFR\_Force (t, f) • conj (TFR\_Force (t, f)) Power\_EMG (t, f) = TFR\_EMG (t, f) • conj (TFR\_EMG (t, f))

Where conj refers to the complex conjugate.

Likewise, the time-frequency cross-spectrum can then be defined as:

$$\text{TFR\\_cspec}\left(\mathbf{t}, \mathbf{f}\right) = \text{TFR\\_Force}\left(\mathbf{t}, \mathbf{f}\right) \bullet \text{ conj}\left(\text{TFR\\_EMG}\left(\mathbf{t}, \mathbf{f}\right)\right)$$

The time course of coherence can then be calculated per frequency as:

$$\frac{\text{TFR\\_Coherence (f)} = }{\text{ } \begin{aligned} \text{TFR\\_cspec (f)} \* W & \text{/} 2 \\ \text{(Power\\_Force (f)} \* W) & \text{\(Power\\_EMG (f)} \* W\) \end{aligned}}$$

Where the term ∗W represents convolution of the indicated time series with a rectangular window (W), the duration of which was set per frequency to be 7/f. The multiplication and division in the above equation are simply element-by-element operations on the time series data.

Prior to further analysis, coherence values were normalized using Fisher's r-to-z transform Fz = atanh(√ C) where C is the coherence at a given time-frequency point (Benignus, 1969).

Because of the short time scales involved in the calculation of wavelet coherence, it is best to recast coherence values as a statistical deviation from chance level. Here, the chance level was derived empirically by recalculating the time-frequency coherence after reversing the concatenated EMG signal in time. This causes EMG signals from one participant to be tested for coherence with force traces produced by a different participant, and completely misaligns the signals with respect to the phase progression of the sinusoidal target. The actual coherence values for each frequency were then converted to standard Z-scores with respect to the distribution of coherence values obtained from the "fake" time series. This method helps to emphasize any synchronization which varies significantly across target phases. Values > 1.65 (the one-sided 95% confidence level for a Z-test) indicate that the time-localized coherence between EMG and force was greater than expected by chance at a given phase and frequency.

#### Simulations

#### Closed-Loop Control Overview

We used a computational model of an afferented muscle to study the dependence of tremor on the dynamics of force production. A schematic diagram of the feedback-driven control loop is shown in **Figure 1B**. Briefly, a Hill-type muscle-tendon model was driven by a neural activation signal to produce force under isometric conditions. The simulation was intended to approximate the action of the FDS muscle in our experimental data. The muscle-tendon model describes changes in force as well as the magnitude and rate of associated changes in the length of the muscle fascicle and tendon, accounting for their viscoelastic properties. Our simulation includes two spinal proprioceptive systems; the muscle spindle and the GTO. Upon muscle fiber lengthening, the muscle spindle sends excitatory feedback through primary (Ia) and secondary (II) afferent fibers proportional to eccentric changes in muscle fiber length and velocity; while GTOs send inhibitory feedback (Ib) proportional to the force in the tendon. A tracking controller, whose operation includes conduction and synaptic delays appropriate for a transcortical loop (Lourenço et al., 2006; Pruszynski et al., 2011; Sohn et al., 2015), sends a command signal (C) to the "spinal cord" which is corrected at each time step according to the difference between the target force level and the actual force output from the muscle. This tracking control signal simply ensured that the afferented muscle-tendon model could follow the target force trajectory, and is not intended to model a specific neural pathway, or to recreate human visuomotor or voluntary tracking behaviors. Signals from the tracking controller, muscle spindle, and GTO, are integrated at the "spinal cord" to generate the αmotoneuron drive to the lumped-parameter muscle model. This neural drive (ND) at each ms (t) can be expressed in the following form:

$$\text{ND(t)} = \text{la(t-15)} + \text{II(t-25)} - \text{lb(t-17)} + \text{C(t)}$$

The output of the tracking controller (C) is calculated as:

$$\text{C(t) = C(t-1) + k \bullet (Target(t) - F(t-35))}$$

where F is the force on the tendon and k is a constant. Note that the above represents a simple "iterative learning control" (ILC) scheme (Wang et al., 2009).

To translate the neural drive into force, the signal was delayed by an additional 22 ms before reaching the muscle fibers to account for conduction time along efferent fibers. At the muscle, the signal was passed through an "activation filter" which shapes the signal to account for calcium dynamics in physiological muscle. Finally, the muscle-tendon model converts this muscle activation signal to the force output of the tendon. In this simulation, the delays for each pathway have been matched to physiological recordings from humans and reflex latencies from the FDS muscle in particular (Lourenço et al., 2006).

The muscle model, muscle spindle, and GTO elements of this control loop have been published previously by various groups and will be described and referenced individually below.

## Control Loop Elements

#### Muscle Model

Our Hill-type muscle-tendon model and its mathematical derivation were adopted from previous literature (He et al., 1991; Brown et al., 1996). The schematic diagram of this muscle-tendon model is presented in **Figure 1B**. The muscle fascicle consists of a mass (M), two passive elastic elements (PE in **Figure 1B**), a viscous element (B), and a contractile element (Fm), which is connected with a pennation angle (α) to a series elastic element (SE) representing tendon and aponeurousis.

The contractile element generates muscle force as a fraction of the maximal force that the muscle is capable of producing. This is defined as the product of its physiological cross-sectional area and a constant factor (45 N/cm2 ) (Holzbaur et al., 2005). Two parallel elastic elements characterize passive behaviors of muscle fascicles. The first (non-linear) spring acts against stretch of muscle fascicle, while the second (linear) spring resists compression (Brown et al., 1996). The series elastic element (SE) shown in **Figure 1B** is a lumped non-linear spring model of tendon and aponeurosis. The force produced by this element in relation to the length of the tendon has been implemented as in Brown et al. (1996). The contraction dynamics within the muscletendon unit are modeled as a second-order differential equation (He et al., 1991).

Taking the above factors into account, the output force function (F) can be summarized as follows

$$\begin{aligned} \mathbf{F(t)} &= \mathbf{MA(t)} \bullet \mathbf{FL(t)} \bullet \mathbf{FV(t)} + \mathbf{F\_PE1(t)} + \mathbf{F\_PE2(t)} \\ &+ (\mathbf{B} \bullet \mathbf{v(t)}) + (\mathbf{a(t)} \bullet \mathbf{M}) \end{aligned}$$

Where MA is the muscle activation (the output of the activation filter), FL is the force-length function, FV is the force-velocity function, F\_PE1 and F\_PE2 are the forces produced by the two elastic passive elements in the model, v(t) and a(t) are the velocity and acceleration of muscle fiber contraction, and the muscle mass (M) and viscosity (B) are constants.

Because the force produced by human participants ranged from 5 to 9% of maximal effort, we applied the same forces to the simulated FDS muscle. Given our focus on understanding the general nature of tremor modulation by dynamic force production, it was not necessary to calculate the precise, isolated contribution of the FDS muscle to the generation of index finger force in our experimental task. To simulate the FDS muscle, architectural parameters were set (**Table 1**) according to published anatomical data (Lieber et al., 1991; Holzbaur et al., 2005). For our purposes, the muscle fibers of the FDS muscle associated with the tendon acting on PIP joint of the index finger were combined into a single belly for simplicity.

#### Muscle Spindle Model

The muscle spindle model employed in this study is adapted from Mileusnic and Loeb (2006). This computational model was chosen because it is both physiologically realistic and, at the same time, is immediately compatible with the inputs and outputs of the other elements within our control loop. The model comprises three types of intrafusal fibers, namely, the bag1, bag2, and chain fibers, all of which are modeled as a second-order mechanical system (a mass, a viscous element, and parallel and series elastic elements), similar to a Hill-type muscle-tendon model. Each of the intrafusal fibers receives input describing the muscle fascicle length, velocity, and acceleration, as well as a fiber-type-specific fusimotor activation signal (dynamic or static). In this study, fusimotor activation was set to be constant during each simulation run. The fusimotor gains tested were 75, 150, and 350. Functionally, these are arbitrary units, but can be expressed conceptually as pulses per second. We chose to define our baseline value as 75, since this is near the previously published value of 70 (Mileusnic and Loeb, 2006), and we varied that parameter because fusimotor drive is known to depend upon task and individual psychology (Ribot et al., 1986; Ribot-Ciscar et al., 2000, 2009; Hospod et al., 2007). Because fusimotor drive is modified by the nature of the task independently of (and even without) α-motoneuron firing (neural drive, in our model), we chose not to assume obligatory α-γ coactivation. It is true that mechanisms other than fusimotor drive may change the effective gain of afferent activity (e.g., presynaptic inhibition). Here, variation in fusimotor drive is not only a likely physiological occurrence, but also serves to more generally

TABLE 1 | Architectural parameters of the slip to the index finger of the flexor digitorum superficialis (FDS) muscle.

represent the overall gain of spindle feedback to motor neurons. For integration with the feedback control loop, the final outputs of the spindle model were normalized to fall between 0 and 1.

#### Golgi Tendon Organ Model

The GTO model was adopted from Elias et al. (2014). This GTO model presents the overall behavior of a population of Ib fibers. It was placed in series with tendon, so that it receives tendon force as an input. The force was then converted into Ib fiber output. The transfer function described in Elias et al. (2014) was implemented using the c2d function in MATLAB. The Ib fiber output was scaled between 0 and 1, as was carried out for the spindle outputs.

#### Activation Filter

The activation filter adjusts the neural drive signal to account for the effects of calcium dynamics (release and reuptake) on crossbridge formation, as described in Song et al. (2008). The resulting muscle activation signal (MA in **Figure 1B**) is the "effective" drive delivered to the muscle model.

#### Simulation Analysis

Force tracking was simulated for 128 s with new values for each output parameter derived every ms. To be certain that only consistent, steady-state behavior was analyzed, only the last 30 cycles were used for analysis. The muscle forces produced by the simulation were analyzed in the same way as the experimentally recorded force tremor, providing a direct comparison.

## RESULTS

#### Tremor during Force Tracking

The phase of voluntary force modulation influenced both low (1–5 Hz) and high (6–15 Hz) frequency bands of involuntary tremor. For reference, the top panels of **Figures 2A,B** show the target isometric force sinusoid, which spanned from 5 to 9% of each individual's maximum voluntary force level. The panels immediately below the target in **Figures 2A,B** show the proportion of total tremor variance associated with each phase of the target sinusoid. The proportion of total force variance accumulated in each 10◦ phase bin is shown for low (**Figure 2A**, upper trace) and high (**Figure 2B**, upper trace) frequency bands. Each plot represents a grand average over all 10 participants, who together tracked a total of 357 target cycles. The largest tremor in either frequency band was observed at the beginning of the rising phase of the target. A second time period of increased tremor amplitude appeared at or slightly after the peak of the target cycle, mainly for high frequency tremor. The dashed line indicates the proportion of variance to be expected in each phase bin if force variability were evenly distributed over all phase bins.

Below these average tremor profiles are the cross-participant standard deviations associated with the mean proportions in the traces above. Variance across participants was highest at the base of the target sinusoid, where bursts of tremor (e.g., **Figure 1**) were often observed. The bottom panels of **Figures 2A,B** depict the number of participants (out of 10 total) whose tremor profiles showed statistically larger (above 0 line) or smaller (below 0

FIGURE 2 | Phase-dependent modulation of tremor amplitude. The amplitude of force variability within the 1–5 and 6–15 Hz frequency bands depended upon the phase of the target cycle (shown at the top of each column, for reference). The first data panel of (A) (left column) shows the mean proportion of total force variability (within the 1–5 Hz band) observed in each 10◦ phase bin, calculated over all participants. The same panel of the right column (B) shows this analysis for 6–15 Hz force fluctuations. The horizontal dashed line in each figure represents the proportion expected per bin if tremor amplitude were constant across target phases. The 1–5 Hz frequency band shows a clear peak at the initiation of the rising phase of the target sinusoid. A similar profile was observed for 6–15 Hz force fluctuations, but with a more pronounced tremor amplification at the peak of the target sinusoid. Below these grand averages are depicted the cross-participant standard deviation for each phase bin. The cross-participant variability is generally highest at the valley of the target sinusoid, when bursts of high amplitude tremor were more likely to occur (as can be seen in Figure 1A). The bottom panels are histograms which depict the number of individuals (out of 10), whose tremor profiles deviated significantly from chance level at each phase-bin, as determined by a Monte-Carlo test (see text for details). In these histograms, counts above the 0 line record the number of individuals who displayed greater than chance-level proportions. Counts below 0 indicate the number of individuals showing lower tremor than expected by chance at each phase. Note that these "less than chance" histogram counts are an effect of other phases showing high proportions of tremor, and do not indicate suppression of ongoing tremor, which was never observed. Histogram counts exceeding the shaded region indicate that significant effects at a given phase were more consistent across individuals than could have occurred by chance (according to a binomial test). The histograms show that the modulation of tremor by target phase was fairly consistent across individuals. The arrows on each histogram (A,B bottom), emphasize that tremor was most often larger than expected at the peaks and valleys of the target sinusoid, for both frequency bands.

line) proportions in each phase bin than expected by chance (as determined by a Monte–Carlo test, as previously described). The shaded region marks the number of participants that may have been expected to show significant effects by chance. That is, histogram counts exceeding the upper limits of the shaded range represent a consistent amplification of tremor occurred across the population of participants. Histogram counts below the 0 line are caused by the high proportions observed in other phase bins and do not represent suppression of tremor, which was never observed. Again, column A shows results for 1–5 Hz force fluctuations and column B shows results for 6– 15 Hz force fluctuations. For both frequency bands, a significant population effect was observed at the beginning of the rising phase of the target sinusoid, and at the beginning of the falling phase.

## EMG to Force Coherence

To confirm that the cross-cycle modulation of force variability was also reflected in the activation of the FDS muscle, we calculated EMG-to-force coherence. Using wavelet coherence, we were able to examine the coupling between signals at each frequency, and at each phase of the target cycle. The statistical magnitude of coherence (z-score with respect chance-level) shown in **Figure 3** for each time-frequency pixel was calculated from the full data set (all 357 tracked cycles). Pixels with values greater than 1.65 can be considered as significant at the 95% confidence level.

Importantly, the coupling between EMG and force signals closely resembles the phase progression of force tremor amplitudes, and reflects the same frequency profile. The phase-related modulation of coherence demonstrates temporal variation in synchronization between signals, which would be expected if the frequency content of neural drive depended on the phase progression of the tracking action.

## Actual Tremor Amplitudes during Force Tracking

Although our study is primarily focused on the modulation pattern of normalized tremor amplitudes as a function of voluntary force dynamics, it is also important to address actual tremor amplitudes, and the extent of amplitude modulation across target phases. Since the pattern of tremor modulation could vary across individuals (described below), we chose to record the maximum, minimum, and 1 amplitude (max-min). The latter was important for better comparability with our simulation results, since our simulation does not contain noise or ongoing physiological tremor, both of which are typically present in human participants. In general, we found tremor amplitudes fluctuated by a factor of about 2 over the course of a target cycle. **Figure 4** shows the mean and cross-participant SD for each measure.

## Tremor Modulation across Different Individuals

Although the phase-dependent modulation of force variability shown in **Figure 2** was representative of the population overall, tremor profiles did vary across individuals. **Figures 5A–D** depicts the tremor proportion profiles for high and low frequency bands in 4 individuals whose profiles differed from each other. Overall, most participants showed some degree of increased tremor

(in either frequency band) at the peaks and/or valleys of the target sinusoid. The modulation of tremor amplitudes in these individuals, shown in the bar graphs at the bottom of each column, indicate that 1–5 Hz tremor amplitude was consistently higher than 6–15 Hz tremor amplitude, but the ratio could vary across individuals.

#### Simulation Results

By embedding a modeled FDS muscle within a feedback-driven control loop, we were able to simulate the force tracking experiment carried out by human participants. Surprisingly, this simple model of an afferented FDS muscle sufficed to produce much the same pattern of tremor modulation in relation to the 0.25 Hz sinusoid as seen in **Figures 2**, **5**. **Figures 6A–C** shows the modulation of tremor obtained when the simulation was run using low, medium, and high fusimotor drive. Adjusting the fusimotor drive, and thus, the gain of afferent feedback from the muscle spindle, could produce variation in simulation results similar to the type of variation observed across different subjects (e.g., compare **Figure 5A** with **Figure 6A**, or **Figure 5C** with **Figure 6C**). Force fluctuations near the valley of the sinusoidal

deviations for these features. For both frequencies, there was a nearly two-fold

modulation of tremor amplitude across a target cycle, on average.

target were present in all cases, although the fluctuations occurring at the peak of the sinusoid was reduced as the afferent gain was increased. As with the experimental data, the rising and falling phases of the target sinusoid did not appear to be associated with consistent changes in tremor activity. At the bottom of each column in **Figure 6** are bar graphs showing the average extent of amplitude modulation. Since the minimum amplitude was nearly 0 in all cases, these bars also represent the average maximum amplitude across phases as well. Of particular importance is the fact that increasing fusimotor gain resulted in a doubling (B) and tripling (C) of high frequency tremor amplitudes, as compared with the low fusimotor drive condition (A). Low frequency tremor did not appear to be consistently influenced, but if anything, was actually reduced in amplitude as fusimotor drive was increased. It should be noted that the amplitudes measured from our simulation should not be expected to precisely match those recorded experimentally. Of greater importance is the relative relationship between high and low frequency tremor amplitudes, and how they vary across individuals or simulation parameters. That said, our simulated amplitudes appear to be smaller than those recorded experimentally by a factor of about 10, which is reasonable considering the highly reduced/simplified nature of the model and the absence of any noise.

We also ran the simulation after eliminating various elements of the control loop. **Figure 7** shows the resulting

tremor modulation pattern when the simulation was run completely feedforward (**Figure 7A**), using only feedback from the controller (**Figure 7B**), using only the controller and GTO feedback (**Figure 7C**), and using only the controller and spindle feedback (**Figure 7D**). Where spindle feedback is present, the fusimotor drive was set to 75 (as in **Figure 6A**). An increase in 6–15 Hz fluctuations occurred roughly at the peak and valley of the target sinusoid regardless of the feedback utilized in the control loop. However, the precise shape, timing, and magnitude of these fluctuations were altered by the type of feedback utilized. Inclusion of spindle feedback (**Figure 7D**) was necessary to produce realistic tremor variance patterns (compared with **Figures 2**, **5**) at the initiation of the rising phase of the target sinusoid. Also it is worth noting that tremor amplitudes (bar graphs at bottom of **Figure 7**) were drastically reduced in the absence of spindle feedback roughly by a factor of 10 for 1–5 Hz tremor and by a factor of about 50 for 6–15 Hz tremor. These observations are well aligned with previous findings where reduction of afferent feedback was associated with reduced/eliminated physiological tremor (Halliday and Redfearn, 1958; Sanes, 1985; Erimaki and Christakos, 2008). Tremor modulation was in general particularly sensitive to spindle feedback, since increasing its effective gain through fusimotor drive (**Figures 6B,C**) or removing spindle feedback

Also, the relationship between high and low frequency amplitude modulation is similar to that observed in the experimental data (Figures 4, 5, bottom).

entirely (**Figures 7A–C**) both influence the overall shape, timing, and amplitude of tremor fluctuations. This occurred even though the largest change in muscle fascicle length was only about 1.1 mm. Interestingly, removal of GTO feedback had minimal effect (**Figure 7D** compared with **Figure 6A**).

## DISCUSSION

In this study, we show that involuntary tremor can arise simply from the dynamic viscoelastic response of afferented muscles during voluntary production of isometric force. We characterized this tremor in two frequency bands as healthy adults performed a sinusoidal force tracking task. Furthermore, we simulated the spontaneous emergence of tremor from purely peripheral mechanisms using a computational closed-loop model comprised of well accepted musculotendon, spindle, and GTO computational modules. Our results extend the current understanding of how force variability arises independently of central mechanisms during production of isometric force. Importantly, our results suggest that simple force tracking tasks may provide a clinically and scientifically relevant

FIGURE 7 | Simulation results after removing sources of feedback control. Each panels (A–D) depicts the results of the simulation after removing elements of the control loop. (A) Depicts the tremor modulation profile that emerged when the simulation was run completely feedforward, with no feedback. (B) Depicts the results when using the controller output, but no neural feedback. Panel (C) shows the tremor modulation occurring when Gogli tendon organ feedback was added to the controller, while panel (D) shows the tremor profile that emerged when only controller and spindle feedback were used. The baseline condition (before removal of feedback) was the same as in Figure 6A. Removal of the GTO feedback had minimal influence on the overall shape and timing of force fluctuations (compare Figure 6A with panel D). In contrast, elimination of spindle feedback not only reduced the magnitude of force tremor, but also influenced the general shape and timing of tremor modulation (panels A–C), particularly at the initiation of the rising phase of the target sinusoid. Shown at the bottom of each column are bar graphs depicting the average tremor amplitude modulation for each frequency band (as in Figures 5, 6). The addition of spindle feedback (D) to the model was the only condition which greatly modified tremor amplitudes. Note the change of scale for the bar graph in (D) as compared with (A–C). While low frequency tremor modulation was increased by roughly a factor of 10 with the addition of spindle feedback, the high frequency tremor increased by about a factor of 50, with respect to any other condition.

window into the neural and mechanical factors which generate involuntary tremor.

Although tremor is not often attributed to the specific dynamics of voluntary force production, several investigations have suggested the existence, and potential importance, of such an interaction. For example, muscle stretch has been suggested to play a role in tremor modulation during dynamic force production. Specifically, declining isometric force is associated with at least some small degree of muscle fiber lengthening, toward resting state (Ito et al., 1998). Compared with shortening contractions, lengthening contractions are associated with increased force variability (Christou and Carlton, 2002) and increased motor unit coherence within the physiological tremor range (Semmler et al., 2002). However, it does not appear that muscle fiber lengthening, or associated spindle activity, could explain our results, since we would have expected a systematic and consistent increase in tremor during the descending phases of our sinusoidal target trajectory.

Similarly, in the production of bite force, 7–10 Hz jaw tremor is reduced in slowly increasing force ramps compared with constant or slowly decreasing force (Sowman et al., 2008). In the present study, a simple relation between force direction (increasing or decreasing) and finger tremor was not observed, likely due to differences in the physiology of bite vs. grip force control. Jaw tremor in the 7–10 Hz range depends upon the activity of periodontal mechanoreceptors (Sowman et al., 2006), and different bite-force dynamics may have led to different levels of dental intrusion (Schoo et al., 1983), and presumably adaptation of the periodontal mechanoreceptors (Sowman et al., 2008). Accordingly, both the dynamics of bite-force production and afferent feedback are important considerations when comparing healthy adults to those who suffer from bruxism (Laine et al., 2015b). While jaw tremor may depend on specific mechanical properties of the gums and their interaction with periodontal mechanoreceptors, these studies do demonstrate that tremor modulation can stem from afferent responses to dynamic force.

Stretch-reflex amplitudes have, in fact, been reported to change during the production of sinusoidal forces. For example, rhythmic (sinusoidal) pen-squeezing has been shown to produce stretch-reflex modulation in the FDS muscle (Xia et al., 2005). This is important in the context of the present study because oscillations of excitation around the stretch reflex loop are considered to be one of the major contributors to physiological tremor (Lippold, 1970; Young and Hagbarth, 1980; Christakos et al., 2006; Erimaki and Christakos, 2008). In the study of Xia et al. (2005), it was observed that stretch-reflex amplitudes were roughly modulated in a sinusoidal fashion such that increased reflex amplitudes were associated with higher background FDS EMG levels. A similar conclusion was reached by Stanislaus and Burne (2009), who reported a consistent relationship between stretch-reflex gain and overall contraction level regardless of force dynamics. If a similar sinusoidal modulation of reflex gain were responsible for tremor modulation in our study, we should have observed a sinusoidal modulation of tremor amplitude, which was not the case.

Few tasks involve only one muscle, as thus, it is possible that some tremors stem from an interaction among co-activated muscles. Due to the simplicity of our task and the posture of the hand, it is likely that any co-activated muscles were also co-modulated synergistically during tracking. It has been shown that synergistic muscles may share neural drive over a wide range of frequencies (Laine et al., 2015a). Moreover, we have shown that changes in the magnitude of an isometric fingertip force are likely produced by a simple scaling of a same muscle coordination pattern (Valero-Cuevas, 2000), while others have shown that isometric force magnitude does not influence the frequency content of shared neural drive among muscles of the hand (Poston et al., 2010). In addition to the potential for shared descending drive, neighboring co-activated muscles would likely show temporally-coordinated afferent activity and reflex responses as well. At the very least, it seems that the effects of slow sinusoidal contractions on musculotendon dynamics would be similar across all co-activated muscles, leading them to tremor at the same time relative to the slow voluntary action.

During voluntary force production, tremor may also stem from the recruitment, de-recruitment, and firing rates motor units. For example, the force level at which motor units are recruited can be lower than the force level at which they are de-recruited (De Luca et al., 1982), and due to the activation of persistent inward currents, the magnitude of neural drive needed to recruit a motor unit is often higher than the drive at which the same unit is de-recruited (Gorassini et al., 2002). Therefore, the population of motor units which generate a given force is partly determined by the recent contraction history of the muscle. Motor unit activity, especially the twitches of motor units near threshold, may contribute to isometric force tremor (McAuley and Marsden, 2000). The relevance of such mechanisms to the present study is not clear, but we can speculate that the contribution would be minimal, given the sufficiency of our simulation (which does not include firing motor units, their intrinsic properties, or any source of signal-dependent noise) to replicate the experimentally-observed tremor modulation. We would, however, assume that simulated tremor amplitudes would more closely match those observed experimentally if signaldependent noise and/or intrinsic motor unit properties were included in our closed-loop system. This is a topic which certainly merits future investigation.

It is of course possible, even likely, that many of the above mentioned sources of force variability were still present to some degree in our study, but were not very consistent across participants or across target cycles. In that sense, they may explain some of the variation in tremor profiles observed across participants. Similarly, the degree to which any differences in muscle/tendon strength, size, and compliance across individuals would have influenced our results remains a topic for future investigation. We can speculate, however, that tendon strain magnitudes would likely be important, given that the dynamics of muscle stretch influence spindle output, which had the largest influence on tremor in the present study. The fact that our model could at least partly replicate cross-subject variation in tremor modulation through manipulation of fusimotor gains, which would be expected to vary across individuals (Ribot et al., 1986; Ribot-Ciscar et al., 2000, 2009; Hospod et al., 2007), adds validity to our simulation results and helps to mitigate concerns about its simplifications/assumptions. While it is beyond the scope of this study to precisely match the tremor profiles of every individual, or to exhaustively test the influence of all possible parameters, our results should serve as an important proof of principle upon which to base future investigation.

Despite the many possible sources of tremor within our task, it is clear that the dominant phase-dependent source of tremor was peripheral neuromechanical coupling, rooted in the viscoelastic properties of muscle and tendon. While smooth tracking is initially disrupted as a mechanical consequence of musculotendon dynamics, the spindle reflex system plays an important role in determining the overall magnitude and timing of the resulting tremor.

Our results, therefore, motivate and justify the development of similar experimental paradigms for scientific and clinical applications. For example, we propose that the tremor induced by slow voluntary force modulation may provide a simple measure of reflex integrity, providing an alternative to direct, yet time-consuming and often uncomfortable, perturbations of nerves or tendons. Further, it may be that characteristic patterns of tremor modulation would emerge within the context of spasticity, dystonia, or within conditions such as Parkinson's disease or essential tremor. As a means of probing peripheral neuromechanical coupling, the type of tremor described in this study may hold potential as a tool for understanding and assessing dysfunctional sensorimotor control in those with congenital or developmental disorders, or in those with acquired dysfunction due to trauma or disease. Finally, neuromechanical coupling may contribute mechanistically to the maintenance or amplification of pathological tremor. While we did not test the effects of inserting a descending tremor-frequency input into our simulation, such investigation may be informative, and perhaps even suggest novel avenues for clinical intervention.

#### REFERENCES

## AUTHOR CONTRIBUTIONS

CL, AN, and FV contributed to the design, execution and drafting of this work, and approved the final manuscript. Experimental data was collected by CL and simulations were implemented by AN.

#### FUNDING

Research reported in this publication was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the National Institutes of Health under Awards Number R01 AR-050520 and R01 AR-052345. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

#### ACKNOWLEDGMENTS

The authors thank Matteo Bertucco for assistance with data collection and participant recruitment.

**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 © 2016 Laine, Nagamori and Valero-Cuevas. 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.

# Beta Band Corticomuscular Drive Reflects Muscle Coordination Strategies

Alexander Reyes 1 , Christopher M. Laine1 , Jason J. Kutch2 † and Francisco J. Valero-Cuevas 1 \* †

1 Brain-Body Dynamics Lab, Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, USA, 2 Applied Mathematical Physiology Lab, Division of Biokinesiology and Physical Therapy, University of Southern California, Los Angeles, CA, USA

During force production, hand muscle activity is known to be coherent with activity in primary motor cortex, specifically in the beta-band (15–30 Hz) frequency range. It is not clear, however, if this coherence reflects the control strategy selected by the nervous system for a given task, or if it instead reflects an intrinsic property of cortico-spinal communication. Here, we measured corticomuscular and intermuscular coherence between muscles of index finger and thumb while a two-finger pinch grip of identical net force was applied to objects which were either stable (allowing synergistic activation of finger muscles) or unstable (requiring individuated finger control). We found that beta-band corticomuscular coherence with the first dorsal interosseous (FDI) and abductor pollicis brevis (APB) muscles, as well as their beta-band coherence with each other, was significantly reduced when individuated control of the thumb and index finger was required. We interpret these findings to show that beta-band coherence is reflective of a synergistic control strategy in which the cortex binds task-related motor neurons into functional units.

#### Edited by:

Sliman J. Bensmaia, University of Chicago, USA

#### Reviewed by:

Gordon M. G. Shepherd, Northwestern University, USA Aaron J. Suminski, University of Wisconsin-Madison, USA

#### \*Correspondence:

Francisco J. Valero-Cuevas<EMAIL_ADDRESS>
† These authors have contributed equally to this work.

Received: 29 December 2016 Accepted: 10 March 2017 Published: 04 April 2017

#### Citation:

Reyes A, Laine CM, Kutch JJ and Valero-Cuevas FJ (2017) Beta Band Corticomuscular Drive Reflects Muscle Coordination Strategies. Front. Comput. Neurosci. 11:17. doi: 10.3389/fncom.2017.00017 Keywords: coherence, dexterity, beta-band, neuromuscular, synergy, EMG, EEG

## INTRODUCTION

Both corticomuscular and intermuscular synchronization, as quantified by coherence analysis (Rosenberg et al., 1989; Farmer et al., 1993; Conway et al., 1995), provide an important means of understanding the cortical drive to muscles (Conway et al., 1995; Baker et al., 1997; Brown, 2000; Boonstra et al., 2009a). Corticomotor drive contains a 15–30 Hz (beta-band) oscillatory component (Murthy and Fetz, 1992, 1996a,b; Sanes and Donoghue, 1993; Stancák and Pfurtscheller, 1996; Donoghue et al., 1998; Mima and Hallett, 1999; Lebedev and Wise, 2000; Witham et al., 2010), which entrains targeted motor neurons (Farmer et al., 1993, 1997; Mima and Hallett, 1999) and leads to synchronization between cortical and muscular activities in that frequency range (Murthy and Fetz, 1992; Conway et al., 1995; Baker et al., 1997; Salenius et al., 1997). Functionally-related muscles also share a common intermuscular beta-band input (Kilner et al., 1999; Boonstra and Breakspear, 2012; Boonstra, 2013), which is widely accepted as cortical in origin (Brown et al., 1999), as the motor cortex is the only well-established source for such beta-band drive.

Although beta-band cortical drive has received a great deal of attention, its functional significance for motor control remains unclear. Currently it is suggested that oscillations in this range functions to support a constant motor state (Kilner et al., 2000; Pogosyan et al., 2009; Engel and Fries, 2010). This is consistent with the observation that, during the production of a constant force, beta-band coherence is strengthened with continued sensory feedback and minimal voluntary movement (Gilbertson et al., 2005; Androulidakis et al., 2006, 2007; Lalo et al., 2007; Engel and Fries, 2010; Aumann and Prut, 2015). While the magnitude of beta-band corticomuscular coherence does correlate with force(Conway et al., 1995; Baker et al., 1997; Kilner et al., 1999, 2000, 2004; Baker, 2007; Kristeva et al., 2007; Witte et al., 2007), it disappears during movement (Baker et al., 1997; Kilner et al., 1999, 2000, 2004; Brown, 2000; Feige et al., 2000) and imagined movements (De Lange et al., 2008), and there is even evidence that the signal is not entirely feed-forward (Fisher et al., 2002; Baker et al., 2006; Witham et al., 2011). Such findings raise important questions as to the functional role of beta-band drive to muscles.

Numerous studies have investigated low force (< 5 N) precision pinch paradigms to characterize cortico-spinal interactions through the use of corticomuscular coherence (CMC) (Muir and Lemon, 1983; Lemon and Mantel, 1989; Lemon et al., 1995, 1998; Baker et al., 2003). Findings have revealed that beta-band CMC is modulated by digit displacement (Riddle and Baker, 2006), object compliance (Kilner et al., 2000), and similar studies suggest a dependence upon the history and time course of muscle contraction (Chakarov et al., 2009; Omlor et al., 2011; Nazarpour et al., 2012). In nearly all cases, findings of decreased beta-band CMC can be interpreted as reflecting either (1) a departure from steady-state control of a particular muscle, or (2) the cortical "unbinding" of muscles when individuated, rather than synergistic, activation is called for. Given that betaband cortical activity has been suggested as a "binding" signal for many years (Gray, 1994; Santello, 2014), and that such binding would naturally favor synergistic rather than individuated control of the fingers (Boonstra et al., 2009b; Danna-Dos Santos et al., 2010; Kattla and Lowery, 2010; Aumann and Prut, 2015), our overall hypothesis was that beta-band corticomotor drive should be reduced or eliminated when the degree of individuated muscle control is increased.

To address these issues, we studied beta-band corticomuscular and inter-muscular coherence (CMC and IMC, respectively) while participants applied low magnitude precision pinch forces to one of two different objects. The first object was a solid wooden dowel. Production of a constant pinch force against a solid object represents a relatively simple task for the nervous system. The second object was a custom-designed spring which buckles when compressed unless prevented from doing so through precise dynamic adjustment of thumb and index fingertip forces. The spring task described herein has been modified from Valero-Cuevas et al. (2003) and described in Dayanidhi et al. (2013b). This task requires the dynamic regulation of thumb and index fingertip force vectors in 3-D to stabilize the spring, which can be modeled as undergoing an instability similar to a subcritical pitchfork bifurcation (Venkadesan et al., 2007). The physical movement of the fingers remains negligible because large movements tend to increase the likelihood of buckling.

If beta-band coherence depends on the generation of a relatively stable pinch force, then coherence should change relatively little across these two tasks. If beta-band coherence is inherently an intermuscular binding signal, then we would expect to see little coherence during compression of the unstable spring, either in terms of CMC or IMC. Our findings support this hypothesis and suggest that the dynamic, mechanical relationships among muscles are likely critical factors shaping the frequency content of corticomotor drive.

## METHODS

## Participants

We recruited 15 healthy participants (30.3 ± 4.6 years, 6 females) who were self-reported as right-handed. There were no known prior or current neurological conditions in any of the participants, nor did they report any previous hand injuries or surgeries. This study was carried out in accordance with the recommendations of Institutional Review Board (IRB) at the University of Southern California (USC) with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Office for the Protection of Research Subjects at USC.

## Task 1: Establishing Maximal Unstable Spring Compression Force

The Strength-Dexterity test provides a quantitative measure of hand dexterity by requiring dynamic regulation of endpoint force direction and magnitude to stabilize a slender and compliant spring prone to bucking (Valero-Cuevas et al., 2003). The instability in the spring increases with compression force, and thus the maximal compression force reached is indicative of the greatest instability the neuromuscular system can control. The force required to bring the spring to solid length was approximately 3.7–3.8 N, however the maximal compression force healthy adults can reach is less than 3.0 N (Dayanidhi, 2012). We specifically chose a spring with a low strength requirement (<15% maximal precision pinch force) to focus primarily on cortical drive involved with dexterity demand, rather than strength. The resting length of the spring measured 4.2 cm in length, weighed approximately two grams and had a spring constant of 0.86 N/cm. Custom designed 3-D printed acrylonitrile butadiene styrene (ABS) plastic end caps were glued to both ends of the spring to create flat surfaces on which to attach a force transducer. Additional ABS end caps were attached on top of the sensor for two purposes: (1) to provide a place for subjects to grasp the object and (2) to serve as a thermal barrier to prevent body heat from adding a bias to the temperature-sensitive force transducer. With the addition of the end caps and sensor, the effective resting length of the spring was 5.7 cm (**Figure 1A**).

In this task, subjects rested their right arm on a table and used their index finger and thumb to compress the spring with their hand resting on the table. During the task, we asked participants to ensure that their 3rd–5th fingers did not assist in the task by tucking them into their palm. They were given four attempts (90 s

tolerance limits as with the spring. Gray areas represent data within tolerance and dowel-low and dowel-high refer to dowel-low and dowel-high, respectively.

each) to try to compress the spring as much as possible (Valero-Cuevas et al., 2003; Dayanidhi et al., 2013a,b; Lawrence et al., 2014). The average maximal compression force reached prior to spring buckling was taken as a normalized measure of their dexterous performance. We rounded this maximal value to the nearest tenth of a Newton and defined this as the subject-specific Fmax. We then calculated 40 and 80% of the subject-specific Fmax for use in the second phase of the study.

#### Task 2: Visuomotor Force Tracking

Subjects were seated comfortably in front of a computer providing visual feedback of their precision pinch compression force. They compressed either the same slender spring or a wooden dowel (length = 5.2 cm, diameter = 0.12 cm, **Figure 1C**) to visually track a series of randomly presented step targets set to 40 and 80% of their Fmax. For three trials for a given object, target force levels were presented at 30 s intervals with 5 s of rest in between (**Figures 1B,D**). The trial lasted until five repetitions of each force level were presented. During the resting periods, the subjects would hold the object with just enough force to prevent dropping it. Breaks were given in between trials when necessary to prevent fatigue effects. Subjects then repeated the procedure with the other object for the same number of trials and force levels. This two-by-two factorial design yielded four force conditions: spring-low, spring-high, dowel-low, and dowel-high.

## Recordings

#### Force

Normal compression forces were measured by affixing a uni-axial load cell (ELFF-B4-10L, Measurement Specialties, Hampton, VA) with double-sided tape to the index finger side of either the compliant spring or the wooden dowel. The circular load cell measured 0.41 cm in height, 1.27 cm in diameter and aligned perfectly with the diameter of the objects. Signals from the sensor were differentially amplified with a custom designed circuit operating in the 0–5 N range. Data were captured using a USB Data Acquisition (DAQ) system (National Instruments, Austin, TX) sampling at a rate of 2,048 Hz. Prior to data collection, the sensor voltage was converted to Newtons by removing the DC offset and calibrating the load-cell using a four-point linear regression with fixed weights. The offset and gain of the load cell were corrected periodically to ensure accurate force recordings.

#### Electromyography (EMG)

Bipolar surface EMG were collected using a Delsys Bagnoli system (Delsys, Natick, MA) from the first dorsal interosseous (FDI) and abductor pollicis brevis (APB). Data were filtered between 20 and 450 Hz, amplified by 1000, and then sampled at a 2048 Hz. The reference electrode for the recordings was placed on the olecranon of the right arm. Recording locations were identified by palpating the muscle during force production in the direction of mechanical action for each muscle.

#### Electroencephalography (EEG)

Sixty four channels of EEG were recorded at a sampling rate of 2048 Hz (ANT Neuro, Enschede, The Netherlands). The fixed recording sites were arranged according to the international 10– 20 system for scalp electrode placement. We ensured repeatable recordings of cortical areas across subject by taking skull measurements and placing electrode Cz at the cross section of the midway point between the nasion and inion and the midway point between the left and right tragus of the ear. Electrode impedances were kept below 10 k with respect to the reference electrode CPz.

Following digitization, EEG signals were bandpass filtered between 10 and 500 Hz and both EEG and EMG signals were notch filtered at 60 Hz and its harmonics up to 500 Hz using a 4th order Butterworth filter implemented in MATLAB (Mathworks, Natick, MA) and FieldTrip, a software package for EEG and EMG analysis (Delorme and Makeig, 2004; Oostenveld et al., 2011). Subsequently, EMG were rectified to extract group activity of motor units (Halliday et al., 1995; Mima and Hallett, 1999). Data formats collected from two separate systems were synchronized by configuring the NI-DAQ to send a trigger pulse to the EEG and EMG systems via a split BNC cable. Custom scripts were created to read in trigger events and synchronize all data.

#### Trial Selection

The start of a trial was defined as the time when the on-screen target transitioned from a resting value to either 40 or 80% of Fmax, and its end was defined as the time when the target value returned back to rest. Trial windows were 30 s in duration with a 5 s inter-stimulus interval. Sample traces of the matching paradigm are shown in in **Figure 1B** for the spring (blue trace) and in **Figure 1D** for the dowel (orange trace); each showing three randomized presentations of the low and high force levels (note that the force levels are the same across objects). Force data during these steady hold phases were visually examined to determine if the task was performed correctly. Our requirement was that the hold phase must be within a ±15% tolerance of the target force value and be held within this range for at least 5 s. The force profiles meeting these criteria are shown as the gray shaded areas in **Figures 1B,D**. Force profiles not meeting these criteria were excluded from analysis. The accuracy of the force matching task was analyzed by taking the root mean square error of the steady hold force compared to the target level for each of the four conditions.

Epochs of synchronized EEG and EMG from each condition that satisfied the steady-state criteria were normalized and pooled across conditions and subjects (Amjad et al., 1997). Normalizing each epoch of data gives equal weight to each section and effectively eliminates the possibility of coherence bias which favors sections with high EMG amplitude (Amjad et al., 1997; James et al., 2008; Schoffelen et al., 2011).

## Coherence Analysis

Synchronous oscillations between cortical activity and EMG indicate functional connectivity which can be assessed through coherence analysis (Nunez et al., 1997; Mima and Hallett, 1999). Briefly, coherence measures the temporal correlation between two signals through the strength of the consistency of their phase lag as a function of frequency. The result is a coherence spectrum bounded between 0 and 1 for each frequency of interest. A value of 1 indicates a perfect temporal correlation, while 0 indicates no correlation.

For a given time series, x(t), let the auto spectrum be represented as

$$P\_{\text{xx}}(f) = \frac{1}{L} \sum\_{i=1}^{L} X\_i(f) \cdot X\_i^\*(f)$$

where Xi(f) represents the Fourier transform of the signal at segment i of L, and ∗ indicates the complex conjugate. A similar spectrum exists for the signal y(t), represented as Pyy f . The cross spectrum between the signals x (t) and y(t) is represented as

$$P\_{\mathbf{x}\mathbf{y}^\*}\left(f\right) = \frac{1}{L} \sum\_{i=1}^L X\_i\left(f\right) \cdot Y\_i^\*\left(f\right).$$

Coherence is calculated by normalizing the square of the crossspectral density between signals x (t) and y(t) by the product of their individual auto spectral densities (Baker et al., 1997; Nunez et al., 1997) as represented by

$$C\_{\rm xy} \left( f \right) = \frac{\left| P\_{\rm xy} \left( f \right) \right|^2}{P\_{\rm xx} \left( f \right) \cdot P\_{\rm yy} \left( f \right)}$$

Corticomuscular coherence (CMC) refers to the specific case where cortical and muscular signals represent x (t) and y(t). CMC was computed for each EEG-EMG electrode pair as well as FDI-APB coherence using FieldTrip, an open-source toolbox in MATLAB for the analysis of EEG and MEG data (Oostenveld et al., 2011). We used discrete prolate spheroidal sequences (DPSS) or Slepian tapers (Slepian, 1978) for the calculation of the auto and cross spectra. The multi-taper method provides several measures of the spectral estimation by multiplying the data series by a series of orthogonal tapers prior to calculating the Fourier transform (Pesaran, 2008). Three tapers were used in our analysis, providing a spectral bandwidth of ±5 Hz (Maris et al., 2007; Schoffelen et al., 2011).

#### Selection of EEG Electrodes

Although CMC can be calculated between each muscle and every scalp electrode, we limited our selection to EEG electrodes that showed high FDI-EEG coherence during the dowel-low task. Similar low-force isometric precision task have been previously used in literature for CMC analysis (Baker et al., 1997; Kilner et al., 2000; Fisher et al., 2002; Riddle and Baker, 2006; Chen et al., 2013). Raw coherence values were first normalized by conversion to standard Z-scores using the following formula

$$Z = \frac{\operatorname{arctanh}\left(\sqrt{\text{C}}\right)}{\sqrt{1/2N}}$$

Where N is the total number of tapers used in the calculation of coherence (C) (Baker et al., 2003; Laine et al., 2013, 2014). The electrode locations selected for further analysis were those at which the average Z-score exceeded a Bonferroni-corrected 99% confidence level (Z = 3.6). The correction accounts for the total number of EEG channels.

## Linear Mixed-Effects Model

A linear mixed-effect model provides a method of describing a relationship for a measurable quantity as a function of the sum of weighted independent variables (Winter, 2013). We investigated the effects of task condition on beta coherence using a linear mixed-effect model with the following format:

$$\text{CMC}\_{\beta} = \beta\_0 + \beta\_1 \cdot \text{Condition} + \beta\_2 \cdot \left(1|\text{participant}\right) + \epsilon\_1$$

where CMCβ is the average beta-range coherence of an epoch, Condition (high and low force for both the spring and dowel) is the fixed-effect term, Participant is the random-effect term, βn terms are the coefficients for the independent variables, and ǫ is the error. The random-effects term was inserted to account for subject variability since several measurements were taken for each condition. Models were generated to estimate CMCβ for the FDI and APB well as beta intermuscular FDI-APB coherence.

## RESULTS

## Baseline Corticomuscular Coherence

To confirm the spatial sensitivity and validity of our analysis procedures, we calculated FDI-EEG coherence for the dowellow condition. **Figure 2A** shows a head map of the average Ztransformed corticomuscular coherence between the FDI and all EEG channels. The locations of the four electrodes that exceeded our significance threshold (described previously) are shown in white over the left primary motor cortex and labeled C3, C1, CP3, and CP1. The spatial localization and magnitude shown over the left sensorimotor cortex coincides with previous literature findings during low force production (Witte et al., 2007; Chakarov et al., 2009; Piitulainen et al., 2013).

Cortical drive to the FDI during the dowel-high and springhigh conditions was directly compared. **Figure 2B** shows the average Z-transformed coherence spectra (5–100 Hz) of the four significant electrodes for the dowel-high (orange trace) and spring-high (blue trace) conditions. The peak coherence for both conditions appeared in the beta frequency range (gray shaded area in **Figure 2B**). The peak Z-transformed coherence during the dowel-high condition was 6.84 at 20.1 Hz. Despite matched force levels, however, the peak coherence during the spring-high condition decreased to 4.59 at 21.1 Hz.

#### Linear Mixed-Effect Model

**Figure 3** shows the results of the linear mixed-effects model which tests the effects of condition on beta FDI-EEG (**Figure 3A**) and APB-EEG (**Figure 3B**) beta-range coherence. Using an Ftest, we compared the difference in effect of spring-low and dowel-low model coefficients on beta CMC. We found no significant difference in effect of low force coefficients on either FDI-EEG or APB-EEG beta coherence. It is important to note that, in these low force conditions, both objects remain in the stable domain as the spring has not been compressed enough to exhibit instability.

A similar comparison of the difference in the linear mixedeffect coefficients for the spring-high and dowel-high conditions revealed a significance difference in effect on both FDI-EEG (p = 6.8459e-08) and APB-EEG (p = 1.6889e-05) beta-range coherence. Overall, CMC was reduced in the unstable task (spring-high condition) compared to the stable task (dowel-high condition) with matched force levels. During the high force compression conditions, the dowel remains in the stable domain, but the spring has been compressed to the point of instability.

## EMG-EMG Coherence

We tested for changes in FDI-APB beta coherence across the higher force conditions. **Figure 4A** shows the FDI-APB intermuscular coherence for the spring and dowel objects. The dowel-high condition showed a peak beta range Z-transformed coherence of 6.59 at 24.6 Hz. This peak value was significantly higher than for that of the spring-high condition, which had a peak value of 2.33 at 24.4 Hz. **Figure 4B** depicts the results of the linear mixed-effect model showing that the effects of the dowel was significantly higher than for the spring at the high force level (p = 4.9415e-10).

## DISCUSSION

In this study, we have examined the functional meaning of betaband corticomotor drive. Specifically, we tested the hypothesis that beta-band cortical drive reflects the "binding" of motor

FIGURE 3 | Results of the linear mixed-effects model for the four conditions: dowel-low (DL), spring-low (SL), dowel-high (DH), and spring-high (SH). The model was constructed to predict mean beta range coherence using condition as the fixed-effect with corrections for the random effects of participant. In each bar graph, the mean beta range CMC is shown on the vertical axis and condition is on the horizontal axis. Bars representing one standard error are shown for each condition and significance above the matched force level column pairs represent the statistical difference in the linear mixed-effect coefficients as determined using an F-test. (A) Linear mixed-effect found no differences in FDI-EEG beta coherence for the low forces, but found a significant difference between the higher force compression conditions. (B) Linear mixed-effect model results for APB-EEG coherence show the same statistical results for low and high conditions as in FDI-EEG.

neurons into functional units for synergistic cortical control, and thus should depend upon the type/degree of muscle coordination required by a task. Our findings provide evidence that beta-band corticomotor drive is inherently a reflection of intermuscular "binding" rather than steady isometric force production.

It is known that synergistically-activated muscles share beta-band cortical drive, as measured by intermuscular EMG coherence (Kilner et al., 1999; Boonstra and Breakspear, 2012; Boonstra, 2013). It has been suggested that this phenomenon reflects the "binding" of these muscles into a functional unit (or synergy) at the cortical level (Gray, 1994; Santello, 2014). As such, we predicted that beta-band CMC should emerge more prominently when the task of the motor cortex can be reduced to scaling the activation of a functional group of motor neurons (Valero-Cuevas, 2000). The idea that cortical output should reflect the dimensionality of a task is in line with the work of Rathelot and Strick (2006, 2009) who proposed that the motor cortex in humans and some higher primates may have evolved specific pathways (an "old" and "new" M1) to allow optional and flexible utilization of muscle synergies or "motor primitives."

In this study we have used compression of a slender spring to obligate time-varying fingertip force adjustment while minimizing physical movement and changes in net force (Valero-Cuevas et al., 2003; Dayanidhi et al., 2013a,b; Lawrence et al., 2014). This paradigm has been shown in an fMRI study to increase the engagement diverse brain networks (Mosier et al., 2011) beyond what is required for less-demanding actions. Succeeding at the unstable task (spring-high) requires a level of independent dynamic control of the muscles of the thumb and index finger to dynamically regulate fingertip force vectors to stabilize the spring (Johanson et al., 2001; Venkadesan et al., 2007), which is not the case for isometric tasks like when a wooden dowel is pinched (Valero-Cuevas, 2000). Accordingly, our data showed a general reduction of intermuscular coherence between the EMG signals recorded from the APB and FDI muscles when the spring was compressed at a level that induced instability. The reduction of beta-band (15–30 Hz) coherence between muscles can be interpreted as stemming from an "unbinding" of these two muscles at the level of the cortex. Changes in force alone would not explain our results, since little change in IMC was observed in the dowel-low vs. dowel high condition. Further, increased spring compression force is associated with increased instability, but not changes in spring compliance. Our results almost certainly reflect a neural reaction to the demands of controlling instabilities.

Similar to our findings for IMC, CMC in the beta-band was markedly reduced for the FDI and ABP muscles when the spring became unstable in the spring-high condition. If this phenomenon were simply due to increased noise in the EEG signal (e.g., due to increased cortico-cortical communication during the more difficult task) we would expect that IMC would have been preserved (since the muscles would still share betaband cortical drive), but this was not the case. Further, a change from static to dynamic force production is known to increase high-frequency (30–50 Hz) CMC as beta-band CMC shifts to higher frequencies (Omlor et al., 2007). We did not observe this type of shift. The simplest interpretation would be that the betaband cortical drive was fundamentally an intermuscular signal in this task, and that a shift in cortical control strategy occurred when the thumb and index finger muscles required independent control. This interpretation is also in line with evidence from a pilot study in which within-muscle motor unit coherence was reduced in the APB muscle during the same high-force spring compression task (Laine et al., 2015). A relationship between beta-band cortical drive and control strategy also agrees with its dependence on psychosensory aspects of a task (Laine et al., 2014), although typically it is the higher frequencies of cortical drive (30–50 Hz) which show the greatest sensitivity to such features (Brown et al., 1998; Omlor et al., 2007, 2011; Patino et al., 2008; Mehrkanoon et al., 2014).

Recently it has been hypothesized that beta oscillations arise from closed loops from M1 to muscle synergies back to M1 (Aumann and Prut, 2015). As such, during sustained contraction of a muscle, groups of sensorimotor neurons oscillate in synchrony to maintain the current state (Engel and Fries, 2010). Conversely during movement, de-synchrony of the local group disrupts beta rhythms. The results obtained in this study support this closed-loop hypothesis given in the fact that during the stable conditions (i.e., spring-low, dowel-low and dowel-high), synergistic muscle activations were necessary for producing the target forces. However, in the spring-high condition, it was necessary to disrupt these synergies despite the fact that a relatively constant net force was maintained.

## Confounds and Alternative Explanations

Physical differences in the objects compressed are of key importance to this study. Our interpretation is that the primary factor was the decreased stability of the spring when compressed at higher forces, however, small differences finger position, movement, object compliance, etc. could influence coherence measures and must be considered.

In a previous study, Kilner et al. found positive correlations between the magnitude of beta-band CMC and object compliance during a ramp and hold precision pinch task (Kilner et al., 2000). In our study, we did not observe a difference in beta-band CMC between the dowel-low and spring-low conditions, and at high forces, we observed a reduction rather than an increase in beta-band CMC. While the study of Kilner et al. shares many similarities with the work presented here, the compliant object used in our study was free at both ends and had the propensity to buckle and slip out of the hand when compressed at applied forces >2.2 N. In contrast, two fixed levers with programmable compliance were compressed in the study of Kilner et al. Such differences make direct comparison between studies inexact, but it is clear that the higher compliance of the spring relative to the dowel should, if anything, favor beta-band coherence (Kilner et al., 2000), as should higher force relative to lower force (Witte et al., 2007). Accordingly, our observation of decreased beta-band coherence in the spring-high condition appears to be related to a change in neuromuscular control strategy in response to the reduced stability of the spring. If very small movements can influence CMC (Kilner et al., 1999), then it is possible that the underlying reason for this may relate to the fact that movement may not allow muscles to be "bound" in the same way as is possible during static isometric force production. Although the physical positioning of the fingers can also have an influence on coherence (Riddle and Baker, 2006) the final position of the fingers during spring vs. dowel compression did not differ by more than about 1 cm, and we have no reason to believe that this would have been a major factor in the present study. Moreover, prior fMRI work in this same task has shown that the presence of instability has an effect distinct from that of compliance (Mosier et al., 2011).

It is also important to consider potential drawbacks of using EEG and surface EMG signals. For example, surface EMG will not be as sensitive as single motor unit recordings (Keenan et al., 2012) and may be influenced by signal processing techniques such as filtering or rectification (Boonstra and Breakspear, 2012; Farina et al., 2013). Even so, a systematic difference in the sensitivity of EMG signals to beta-band drive across tasks seems an unlikely explanation for our results, especially given the low levels of force required and our use of the most appropriate signal processing methods (i.e., EMG rectification) under these conditions (Farina et al., 2013).

It is relevant that compression of the spring to the point that it becomes unstable is an inherently difficult and demanding task. While there is increasing evidence that M1 activity may be involved with the perception of task goals (Shen and Alexander, 1997; Cisek et al., 2003; Scott, 2003), our suggestion is that the perception of difficulty in motor tasks is secondary to more tangible factors, such as the need for time-critical sensorimotor corrections, and individuated vs. synergistic control of muscles. In one recent study, for example, beta-band CMC was reduced when the degree of bimanual muscle coordination was increased through a visual feedback manipulation (de Vries et al., 2016). The interpretation was that beta-band CMC relates to the control of individual muscles, rather than coordination. The opposite argument could also be made, however, because in that study, errors of force may have been perceived and corrected individually for each hand in order to achieve better overall bimanual coordination. This would be in line with our coordination-based interpretation of beta-band drive to muscles. The same study also found that low-frequency IMC (∼10 Hz) across hands was highest when bimanual coordination was highest. Because ∼10 Hz CMC is not present in that study nor our own, low-frequency drive to individual muscles cannot be measured or compared concurrently with inter-muscular drive (i.e., IMC). This makes it difficult to attribute changes in IMC to inherently coordination-related aspects of our tasks, rather than other factors which might influence the production of lowfrequency neural drive in general. That said, our **Figure 4A** does visibly show low-frequency IMC reduced along with the beta-band in the spring-high condition, which certainly justifies further investigation of this issue.

Overall, our study tests the notion that beta-band cortical drive essentially reflects the dimensionality of cortical commands, and our results strongly suggest that beta-band CMC should be interpreted carefully, with special attention to

#### REFERENCES

muscle coordination and time-critical sensorimotor demands. Our current speculation, supported by the findings of this investigation, is that beta-band CMC reflects the use of a lowdimensional mode of cortical control over groups of muscles, rather than the maintenance of a steady-state force output. While fully understanding the neurophysiology and task-dependence of cortico-motor oscillations requires further study, our results justify such future work and provide a springboard for more focused investigations into the relationship between the physical requirements of a task and neural control strategies necessary to satisfy them. Finally, such insights into the origin and modulation of cortico-motor oscillations would not only clarify fundamental mechanism for sensorimotor control, but perhaps provide well-founded tasks and analyses directly translatable to clinical measures and diagnostic tests (Norton and Gorassini, 2006; Hammond et al., 2007; Pogosyan et al., 2009; Fisher et al., 2012; Ko et al., 2016).

## AUTHOR CONTRIBUTIONS

AR, JK, and FV designed research; AR performed research; AR, CL, and JK analyzed data; AR, CL, JK, and FV wrote the paper.

#### ACKNOWLEDGMENTS

This material is in part based upon work supported by NSF Grant EFRI-COPN 0836042, NIH Grant R01-052345, NIH Grant R01-050520 to FV, and Supplement R01-050520-W1 to AR.

coupling between neuronal spike trains. Prog. Biophys. Mol. Biol. 53, 1–31. doi: 10.1016/0079-6107(89)90004-7

**Conflict of Interest Statement:** FV holds US Patent No. 6,537,075 on some of the technology used in this study that is commercialized by Neuromuscular Dynamics, LLC.

The other 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 Reyes, Laine, Kutch and Valero-Cuevas. 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.

*Keivan Mojtahedi 1 , Bryan Whitsell 2 , Panagiotis Artemiadis 2 and Marco Santello1 \**

*1Neural Control of Movement Laboratory, School of Biological and Health Systems Engineering, Arizona State University, Tempe, AZ, USA, 2Human Oriented Robotics and Control Laboratory, Mechanical and Aerospace Engineering, School for Engineering of Matter, Transport and Energy, Arizona State University, Tempe, AZ, USA*

Of particular interest to the neuroscience and robotics communities is the understanding of how two humans could physically collaborate to perform motor tasks such as holding a tool or moving it across locations. When two humans physically interact with each other, sensory consequences and motor outcomes are not entirely predictable as they also depend on the other agent's actions. The sensory mechanisms involved in physical interactions are not well understood. The present study was designed (1) to quantify human–human physical interactions where one agent ("follower") has to infer the intended or imagined—but not executed—direction of motion of another agent ("leader") and (2) to reveal the underlying strategies used by the dyad. This study also aimed at verifying the extent to which visual feedback (VF) is necessary for communicating intended movement direction. We found that the control of leader on the relationship between force and motion was a critical factor in conveying his/her intended movement direction to the follower regardless of VF of the grasped handle or the arms. Interestingly, the dyad's ability to communicate and infer movement direction with significant accuracy improved (>83%) after a relatively short amount of practice. These results indicate that the relationship between force and motion (interpreting as arm impedance modulation) may represent an important means for communicating intended movement direction between biological agents, as indicated by the modulation of this relationship to intended direction. Ongoing work is investigating the application of the present findings to optimize communication of high-level movement goals during physical interactions between biological and non-biological agents.

Keywords: human–human interaction, impedance, leader and follower, physical interaction, stiffness

## INTRODUCTION

Collaboration, defined as the act of cooperation among multiple agents toward the attainment of a common goal, is one of the most sophisticated behaviors exhibited by biological organisms. Although cooperation is ubiquitous among a wide range of species ranging from ants to primates, the level of sophistication reached by humans in their ability to cooperate is unparalleled in the animal kingdom. Of particular interest to the neuroscience and robotics communities is the understanding of how humans collaborate to perform motor tasks.

Physical collaboration between two homologous biological agents, such as two humans holding a tool or moving it across locations, entails complex sensorimotor processes. Specifically, the

#### *Edited by:*

*Matthew Tresch, Northwestern University, USA*

#### *Reviewed by:*

*Robert Gregg, The University of Texas at Dallas, USA Marc Daniel Killpack, Brigham Young University, USA*

#### *\*Correspondence:*

*Marco Santello<EMAIL_ADDRESS>
*Received: 30 September 2016 Accepted: 29 March 2017 Published: 13 April 2017*

#### *Citation:*

*Mojtahedi K, Whitsell B, Artemiadis P and Santello M (2017) Communication and Inference of Intended Movement Direction during Human–Human Physical Interaction. Front. Neurorobot. 11:21. doi: 10.3389/fnbot.2017.00021*

problem of physically collaborating with another agent to perform a given motor task introduces control problems that go well beyond those encountered when controlling one's own limb. For example, planning and execution of reaching or grasping movement are thought to occur through an internal model of the agent's limb that allows prediction of the sensory consequences of the motor action (Johansson and Flanagan, 2009; Wolpert et al., 2011). Examples of such phenomena are the temporal coupling of grip and load forces associated with moving an object denoting anticipation of movement-related inertial forces (Flanagan and Wing, 1997), or the anticipatory control of torque prior to lifting an object with an asymmetrical center of mass (Salimi et al., 2003; Bursztyn and Flanagan, 2008; Fu et al., 2010, 2011; Fu and Santello, 2012; Mojtahedi et al., 2015). However, when two humans physically interact with each other, sensory consequences and motor outcomes are not entirely predictable as they also depend on the other agent's actions. Therefore, the question arises as to how the central nervous system of each agent factors in the other agent's actions when physically interacting with each other to perform a collaborative task. A better understanding of this problem can help developing biologically inspired controllers supporting human–robot physical interactions, e.g., exoskeletons used for neurorehabilitation or physical augmentation, and optimizing the way these interactions can be performed.

Physical interaction between humans and robots has been mainly investigated using the notion of mechanical impedance. Hogan (1985) first proposed robot impedance controllers as a way to guarantee stable and robust behavior of a robot that interacts with a human. Since then, a plethora of robot applications involving physical human–robot interaction use control of impedance, and in most cases, this is done to purposefully impose a specific dynamic behavior to the human agent. For example, the MIT-MANUS—used extensively for upper limb rehabilitation—uses the concept of impedance control in a back-drivable system to restrict the motion of the patient's arm along a specific path, while the patient tries to reach a target *via* a manipulandum attached to his/her paretic arm (Krebs et al., 1998). For this scenario, impedance control is used to assist the human subject to reach a pre-defined target and imposes high resistive forces to motion that is not congruent with the desired trajectory.