University of Southern California Dissertations and Theses

Coordination in multi-muscle systems: from control theory to cortex

(USC Thesis Other)

Coordination in multi-muscle systems: from control theory to cortex

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content COORDINATION IN MULTI-MUSCLE SYSTEMS: FROM CONTROL
THEORY TO CORTEX
by
Sarine Babikian
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
August 2017
Copyright 2017 Sarine Babikian
Dedication
To my mother, father, and brother.
ii
Acknowledgements
There are many individuals who have helped, encouraged, and supported me throughout
my time as a PhD student. The day-to-day stress of graduate school makes it easy to
forget what a dream come true this is; but it would not have been possible without the
support of those people and for that I am eternally grateful.
First and foremost, I would like to thank my parents, Maral and Harout, for their
endless support in following my dreams. I want to thank my mother for always supporting
me in my education and career (but also for providing me with her delicious recipes). I
want to thank my father for being so generous and helpful (and for frequently asking me
if I've been eating). I can not forget to thank my brother Sevag; I am lucky and proud
to have someone so smart and dedicated whom I can always learn something from. I am
also extremely appreciative of my extended family in Los Angeles who welcomed me with
open arms. My aunt Vergine, grandparents, and cousins who helped me out so much
especially when I rst moved here. My Lebanese cousins, Talar, Shaghig and Patil for
all the lighthearted group chats they contribute to that always make my day (thanks for
keeping the inside jokes running).
Second, I want to thank my advisor and mentor, Eva Kanso. I met Eva in 2008 when
I was an undergraduate summer researcher in her lab. Since then, she has helped open
iii
up many opportunities over the years. I would not be part of this program, and I would
not be doing the research I enjoy doing, had it not been for Eva. She has never missed
the chance to give useful advice, while also delivering constructive criticism with a good
sense of humor. I am incredibly grateful to have met and worked with such an amazing,
caring and smart woman who has become a role model.
I'd like to also thank my co-advisor, Jason Kutch. I am extremely lucky to have had
the opportunity to work with someone who is the true denition of a scientist. Over the
past two years where I've been advised by Jason, I have come to realize what a supportive
and caring mentor he is, who not only teaches students how to be good scientists (and
extremely organized scientists), but also helps us develop new skills that we never thought
we would have the opportunity to learn. I am thankful for being given the chance to work
with him and for all the help and encouragement he has provided.
I'd like to thank my committee members, Profs. Paul Newton and Niema Pahlevan,
for their valuable input on my work; and previous research collaborators, Francisco Valero-
Cuevas and members of the Brain-Body Dynamics Lab. I want to thank Alex Reyes for
being a mentor and supportive friend.
It's so easy to take for granted the people I see every day; in reality they are the
ones who keep me going and motivated. I want to thank my lab mates who have become
friends; Hanliang, Yangyang, Alan, and all other previous and current members of the
Biodynamics Lab. Our postdoc and my good friend/Pasadena neighbor Lionel and his
wife Kristina, for all the dinners and lively conversations we've had. Members of the
Applied Mathematical Physiology Lab, Moheb, Sonja and Skulpan. All those who helped
by participating in my studies and allowing me to record their brain signals. My friends
iv
from biomedical engineering and previous roommates, Amber and Cassie: thank you for
being some of the coolest women in the engineering PhD program. All other unnamed
people I met at USC who were there one way or another since the start of graduate
school. My Lebanese friends Roubina, Ramela and Nanar for still engaging in meaningful
conversations even from miles away, but most importantly, for keeping me up-to-date with
all the gossip from back home.
The past few months have been some of the busiest and most stressful months of
my life. I couldn't have done it without the support of Aram, who even while being
thousands of miles away, has been incredibly encouraging and patient when I constantly
felt overwhelmed. Thank you for listening to my endless rants and for always believing
in me.
v
Table of Contents
Dedication ii
Acknowledgements iii
List Of Tables viii
List Of Figures ix
Abstract xii
Chapter 1: Introduction 1
Chapter 2: Energy optimization and non-smoothness in bio-inspired limb
movements 6
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Chapter 3: Experimental evidence of non-smoothness in nger motion 32
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
Chapter 4: Goal-equivalent nger force variation in EEG activity 39
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Chapter 5: Voluntary nger force variation in EEG activity 60
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
vi
Chapter 6: Goal-equivalent nger force variation in fMRI brain activity 70
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
Chapter 7: Goal-equivalent variability explained by interhemispheric in-
teractions 81
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Chapter 8: Conclusions 104
Appendix Chapter A: Common Spatial Patterns and Linear Discriminant
Analysis for the EEG prediction of movement intention 108
A.1 Common Spatial Patterns (CSP) . . . . . . . . . . . . . . . . . . . . . . . 108
A.2 Linear Discriminant Analysis (LDA) . . . . . . . . . . . . . . . . . . . . . 110
A.3 Classication of simulated EEG example . . . . . . . . . . . . . . . . . . . 118
Bibliography 123
vii
List Of Tables
4.1 AUC values and number of repetitions in each class for each participant. . 53
A.1 LDA example: product dimensions and quality control result . . . . . . . 115
viii
List Of Figures
2.1 Finger and leg models showing kinematics and muscle routing. . . . . . . 9
2.2 A sketch of the static equilibrium and energy optimization problem. . . . 14
2.3 Variable pre-tensioning update rule. . . . . . . . . . . . . . . . . . . . . . 18
2.4 Finger model conguration and endpoint space. . . . . . . . . . . . . . . 20
2.5 Endpoint trajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Finger model: optimal stiness values for the quasi-static trajectories. . . 22
2.7 Leg model: optimal stiness values for the quasi-static curved trajectory. 23
2.8 Finger model with constant moment arms: optimal stiness parameters
needed to achieve each posture in conguration and endpoint space. . . . 24
2.9 Finger model with posture-dependent moment arms: optimal stiness pa-
rameters needed to achieve each posture in conguration and endpoint
space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.10 Leg model with constant and posture-dependent moment arms: optimal
stiness parameters needed to achieve each posture in conguration space. 26
2.11 Finger model: achievable versus desired trajectories overlaid on top of the
reachable postures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1 Accelerometer setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2 Averaged nger accelerations. . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Power spectral densities of motion trials for all four speeds and two types
of motion (normalized across all speeds). . . . . . . . . . . . . . . . . . . . 37
ix
4.1 EEG Experiment setup and force output. . . . . . . . . . . . . . . . . . . 44
4.2 Covariance ellipses and force variance ratio of 18 EEG participants. . . . . 49
4.3 EEG-based prediction of force output and averaged common spatial pat-
terns (CSP). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.4 Individual ROC curves of all EEG participants and overall ROC curve
across all participants. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.1 EEG Experiment setup and force output. . . . . . . . . . . . . . . . . . . 64
5.2 Covariance ellipses and force variance ratio of 18 EEG participants. . . . . 66
5.3 EEG-based prediction of force output and averaged common spatial pat-
terns (CSP). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.1 fMRI Experiment setup and force output. . . . . . . . . . . . . . . . . . . 73
6.2 Covariance ellipses and force variance ratio of 9 fMRI participants. . . . . 76
6.3 Example participant BOLD ratio as a function of force ratio. . . . . . . . 77
6.4 All participants' BOLD ratio as a function of force ratio. . . . . . . . . . . 78
7.1 Interhemispheric inhibition measured by TMS. . . . . . . . . . . . . . . . 83
7.2 Interhemispheric inhibition in the hand motor area. . . . . . . . . . . . . . 85
7.3 A sketch of the interhemispheric connections in the M1-M1 model. . . . . 87
7.4 Example of the sigmoid feedback function in the M1-M1 model. . . . . . . 89
7.5 Bifurcation analysis of mutually-inhibiting neuronal network. . . . . . . . 90
7.6 Example phase-plane of the bistable system. . . . . . . . . . . . . . . . . . 91
7.7 Bistable system repetitions without noise. . . . . . . . . . . . . . . . . . . 92
7.8 Bistable system repetitions with the addition of noise. . . . . . . . . . . . 94
7.9 Example simulation and covariance ellipse. . . . . . . . . . . . . . . . . . . 96
7.10 Covariance ellipses corresponding to 1000 simulations. . . . . . . . . . . . 97
x
7.11 Inhibition functions of valid and invalid parameter sets . . . . . . . . . . . 99
7.12 Relationship between inhibition scaling and input amplitude. . . . . . . . 100
7.13 Boxplot showing the relationship between inhibition scaling and input am-
plitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
A.1 LDA example: features before projection . . . . . . . . . . . . . . . . . . . 116
A.2 Projected features with LDA . . . . . . . . . . . . . . . . . . . . . . . . . 117
A.3 Spatially ltered EEG signals . . . . . . . . . . . . . . . . . . . . . . . . . 120
A.4 Projected time series of simulated EEG signals . . . . . . . . . . . . . . . 121
A.5 Class prediction of simulated EEG . . . . . . . . . . . . . . . . . . . . . . 122
xi
Abstract
It remains unknown how the human brain coordinates the many biomechanical degrees
of freedom (e.g. muscle activations) to generate movements. Many theories explaining
motor coordination have been proposed, ranging from optimizing for a unique muscle
activation pattern, to exploring the subspace of plausible activation patterns that equally
accomplish the same task. This dissertation focuses on investigating the feasible solution
space of muscle activation patterns that produce the same movement or force production.
We use both computational and experimental approaches to understand the mechanisms
and neural origins of varied patterns that achieve the same task.
We rst present a mathematical model of a tendon-driven limb controlled using energy
minimization in muscles; and investigate the conditions in muscle parameters that allow
for the stability of limb postures and transitions between postures to generate slow move-
ment. We nd that the precision with which individual muscle activations are altered
changes substantially during the trajectory to produce smooth limb motions, suggesting
a possible explanation to the non-smoothness observed experimentally in slow movements
of human ngers.
We then turn our attention to the observation that human movements are performed
reliably and repeatedly while the details of the movement may vary at each repetition.
xii
This goal-equivalent variability has been extensively studied at the level of motor output
and has been shown to be crucial in avoiding injury risk. However, the neural mecha-
nism producing this variability remains unknown. To this end, we measure brain signals
using two modalities, electroencephalography (EEG) and functional Magnetic Resonance
Imaging (fMRI), while healthy subjects perform a bimanual task involving repeatedly
applying left and right index nger forces to keep the sum of forces constant. Using ma-
chine learning with EEG, we build predictive models of goal-equivalent variability based
on cortical signals, and nd motor cortical representation in the signals contributing to
the variability. Using fMRI with much higher spatial resolution, we verify our EEG re-
sults and conclude that motor cortical activity is correlated with automatic variation in
nger force output.
Finally, we hypothesize that goal-equivalent variation in left and right nger forces
is driven by interhemispheric interactions between brain signals in contralateral motor
cortices. We build a mathematical model of goal-equivalent variability based on the
interhemispheric inhibition between left and right motor cortices of the human brain,
and tune parameters to match our experimental results of nger force variation. Our
results suggest that in addition to the non-linearity requirement of the inhibition function
observed in previous research, a linear relationship between the level of the target force
and the strength of interhemispheric inhibition must also be met.
xiii
Chapter 1
Introduction
Human movements are complex in the sense that there exist many more degrees of free-
dom (i.e., muscle activations, forces) than needed to control a movement or a force pro-
duction task. A long standing problem in studying human motor control is understanding
how the many degrees of freedom are coordinated to perform movement (Bernstein 1967);
in particular, how the human central nervous system (CNS) chooses one out of many pos-
sible solutions to control movements. This redundancy in the human body can be found,
for instance, at the level of individual joint rotations (Bernstein 1967, Berkinblit, Feld-
man & Fukson 1986) and at the level of individual muscle activation patterns (Buchanan,
Rovai & Rymer 1989, Nichols 1994). Many dierent theories of motor coordination have
been developed to address motor redundancy; some attempting to eliminate redundant
degrees of freedom using computational principles that optimize for a unique solution,
while others suggesting that the CNS uses dierent solutions each time the movement is
repeated, thus allowing for a manifold of plausible solutions that may correctly perform
the same task.
1
The optimization approaches usually involve nding a unique solution to control move-
ment by minimizing or maximizing a cost function, often based on mechanical, physiolog-
ical or more complex aspects (Seif-Naraghi & Winters 1990, Rosenbaum, Loukopoulos,
Meulenbroek, Vaughan & Engelbrecht 1995, Prilutsky & Zatsiorsky 2002, Todorov 2004).
For example, optimal control modeling in biomechanics most often focuses on minimizing
the energy used by the muscles (Anderson & Pandy 2001, Chow & Jacobson 1971, Collins
1995, Davy & Audu 1987, Hatze & Buys 1977, Popovic, Stein, Oguztoreli, Lebiedowska
& Jonic 1999, Rasmussen, Damsgaard & Voigt 2001). Others utilize a smoothness cost,
e.g. in modeling control of arm movements, and minimize the time-derivative of accel-
eration (Flash & Hogan 1985, Hogan 1984b, Pandy, Garner & Anderson 1995, Smeets
& Brenner 1999, Todorov & Jordan 1998), joint torque (Nakano, Imamizu, Osu, Uno,
Gomi, Yoshioka & Kawato 1999, Uno, Kawato & Suzuki 1989), or muscle force (Pandy
et al. 1995).
Another theory developed to understand motor redundancy has derived from the ob-
servation that repeated movements, even when performed exactly the same way each time,
are characterized by varying movement details at each repetition (Bernstein 1967, Latash
2008). In other words, dierent muscle activation patterns may be used to repeat the same
arm-reaching movement. The nervous system appears to allow dierences in the ne de-
tails of the movement while performing the same task correctly at each repetition (Latash
& Anson 2006). An optimal feedback model where the only feedback is to correctly reach
the goal instead of instantaneously controlling for error during movement, has shown
agreement with experimental data where variable trajectories were observed across rep-
etitions, while the task (or goal) was achieved without error (Todorov & Jordan 2002).
2
Furthermore, there is increasing evidence that this \goal-equivalent" variability, rather
than being noise, is associated with lower injury risk or pain because it decreases muscle
fatigue and tissue loading (Hamill, Palmer & Van Emmerik 2012, Lipsitz 2002). As such,
one study demonstrated evidence that healthy runners used more variability in trunk
to pelvis rotations than runners with a history of low back pain and those with current
low back pain (Seay, Van Emmerik & Hamill 2011). Others have found low variability
in patients with knee pain (Hamill, van Emmerik, Heiderscheit & Li 1999), and elderly
populations vulnerable to falls (Lipsitz 2002).
This dissertation presents both an optimization approach to motor coordination em-
ploying an energy minimization technique, and a goal-equivalent approach where vari-
ability of movement patterns across repetitions is explored at the level of brain signals.
In Chapter 2, we present a mathematical model focusing on the energy minimization of
a tendon-driven limb, exploring the role of precision in the control of individual param-
eters in producing smooth and slow movement. Our main result her is that producing
smooth endpoint trajectories often requires dierences in the precision of muscle param-
eter control. We suggest this as a possible neural explanation to experimentally-observed
non-smoothness in slow nger movements in humans. Based on these experimental re-
sults, we then further explore this non-smoothness in Chapter 3 where we investigate
index nger accelerations during dierent speeds of movement performance in humans.
Chapters 4 to 7 focus on the goal-equivalent aspect of human movements. Goal-
equivalent variation has been well-studied at the level of motor output (kinematics or end-
point forces) within dierent tasks such as multi-nger force production, multi-joint reach-
ing, standing and stepping (Latash 2008); and at the level of muscle commands (Krishnamoorthy,
3
Latash, Scholz & Zatsiorsky 2003). However, the mechanism controlling this variation in
the nervous system has not been established. Thus, in Chapters 4, 5 and 6, we measure
human brain signals using two modalities, electroencephalography (EEG) and functional
Magnetic Resonance Imaging (fMRI) while participants perform a two-nger force match-
ing task; with the aim of investigating the neural mechanism(s) involved in the production
of goal-equivalent variability. The bimanual task requires applying forces using both in-
dex ngers to repeatedly match a single target of the sum of both forces, which results
in a clear variability in output forces. Our EEG signals indicate that goal-equivalent
nger force variation is correlated with cortical signals, in particular, localized hand mo-
tor cortex activity in generating nger force variation. We use fMRI for better spatial
resolution of brain signals and analyze Blood Oxygen Level Dependent (BOLD) signals
from subjects performing the same task in an MRI scanner. We nd that left and right
motor cortical activity is correlated with right and left automatic nger force variations,
thus validating our EEG results.
Lastly, in Chapter 7, we develop a mathematical model describing neural connections
in the brain that may produce goal-equivalent variability. We hypothesize that the vari-
ability in the brain is driven by interhemispheric interactions between brain signals in
contralateral motor cortices. We base our model on the experimental observations in
numerous studies which concluded that there exists mutual inhibition between the two
motor cortices, and that this inhibition follows a non-linear shape with respect to con-
tralateral brain activity. Our mathematical model of mutually-inhibiting motor cortices
with the addition of noise produces goal-equivalent variability similar to the nger force
output we observed experimentally. We tune the parameters of the model to match our
4
experimental results, and nd condition(s) that the parameters must meet to make the
output experimentally valid. Our results suggest that while non-linearity of the inhibition
function is a necessary condition, it is not enough for the output to be experimentally
valid. Rather, the strength of the mutual inhibition must be proportional to the experi-
mental task demand, e.g., the target force to be matched.
5
Chapter 2
Energy optimization and non-smoothness in bio-inspired
limb movements
2.1 Introduction
Bio-inspired tendon-driven limbs are widely utilized in the development of several robotics
devices, from robot hands and arms, to end eectors for surgical robots (Lee, Choi, Chung
& Youm 1994, Jung, Kang, Lee & Moon 2015, Frecker & Snyder 2005, Simaan, Xu,
Wei, Kapoor, Kazanzides, Taylor & Flint 2009). Many of these mechanisms employ an
actuation based on the length change of elastic { spring-like { tendons, e.g., by using
motors that behave like linear springs to generate movement. Changing the stretched
length of elastic tendons generates a force that in turn creates a torque at the joints.
The same force may be produced by changing the spring constant, or stiness, of the
tendon. Here we propose a computational framework for slow movements of tendon-driven
limbs where the controlled parameters are the tendon stiness values with a particular
resolution.
6
Accurate and slow limb movements are quasi-static in the sense that the limb reaches
a new state of equilibrium instantaneously in response to the control input. Here, the
limb is neither in a state of xed equilibrium under constant tendon forces, nor in a
dynamic state with time-dependent tendon forces resisted by inertial accelerations of the
limb. Inertial accelerations are negligible for slowly-actuated, quasi-static, limb motions.
These slow motions belong to a class of mechanical systems, known as driftless, where
motion stops when actuation stops. Driftless systems arise in many applications, including
satellite dynamics (Teel, Murray & Walsh 1995, Boscain & Chitour 2002), robotic vehicles
(De Luca, Oriolo & Samson 1998), and terrestrial and underwater biolocomotion (Shapere
& Wilczek 1989, Kanso, Marsden, Rowley & Melli-Huber 2005, Kanso 2009, Jing &
Alben 2013). We apply these principles to examine the neuromechanical properties of
quasi-static limb movements. A multi-muscle, multi-articular tendon-driven system has
an actuation space (tendon stiness parameters) with a dimension that is larger than
the number of kinematic degrees-of-freedom. This is in contrast to most studies on
driftless systems that consider underactuated motions, i.e., fewer actuation parameters
than kinematic degrees-of-freedom (Kanso et al. 2005, Jing & Alben 2013). It is also
unlike most robotics-type systems that have strictly 2 opposing tendons per kinematic
degree-of-freedom, or analyses that study limb function using rigid body dynamics driven
by torques applied directly at the joints (Hogan 1985). In tendon-driven limb movement,
the forward problem of controlling the stiness parameters that produce a desired slow
limb movement is non-trivial because the solution is not unique. The rotation of each
limb joint does not uniquely determine the lengths of all tendons crossing it.
7
We formulate the problem of accurate and slow movements in tendon-driven limbs
as an overactuated (i.e., underdetermined) driftless system, controlled by changing the
tendons' stiness parameters. We follow in the tradition of simple models that address
the fundamental physics of biological function. Our mathematical approach is reminiscent
to the rate control method proposed by Whitney (Whitney 1969). It is also related to
the control of limb impedance, another mathematical formulation well known in robotics
(Hogan 1985). Here, we focus, from a physics and mechanics perspective, on the range
of possible functions given the minimal assumption that muscle stiness can be changed.
We demonstrate that the control of stiness suces to produce stable and accurate limb
postures and quasi-static transitions among them. This is in agreement with the results
from impedance control. However, we go beyond these results to oer a proof that stable
postures are achievable only when tendons are pre-tensioned.
Altering tendon stiness to produce equilibrium limb postures does not admit a unique
solution due to the over-actuated nature of the tendon-driven system, thus, the same
posture can be achieved with multiple combinations of stiness parameter values. Each
combination places the system at a distinct strain energy level. Thus a same posture may
be realized with dierent levels of strain energy. Therefore, we formulate the problem
of stiness control as an optimization problem (Ivaldi, Morasso & Zaccaria 1988, Mussa-
Ivaldi & Hogan 1991) that minimizes the strain energy of the limb. Further, we nd
that, for constant joint moment arms and tendon pre-tensioning values, optimal stiness
values depend on the change in the joint angle but not on the reference posture. This
symmetry with respect to reference postures is broken in the case of posture-dependent
moment arms.
8
m
17
m
8
m
18
m
4 m
23
m
20
m
22
d
m
4
m
1
m
3
m
2
a
θ
1
θ
2
l
1
l
2
y
x
y
x
θ
1
θ
2
θ
3
l
1
l
2
l
3
b
c
Figure 2.1: Finger (a,b) and leg (c,d) models: kinematics (a,c) and muscle routing (b,d)
overlaid on top of the physiological system.
Finally, we investigate the smoothness of quasi-static transitions among limb postures
when the precision with which the stiness is controlled is not innitely smooth, that is to
say, when the stiness parameters are constrained to vary by small but nite increments.
Under such constraints on precision, we observe discontinuities in the limb's quasi-static
trajectories, indicating the existence of nearby unreachable postures because of insu-
cient stiness precision. We also nd that the precision with which stiness needs to be
controlled depends non-uniformly on the location in the limb's workspace.
2.2 Methods
We examine the slow movements of tendon-driven limbs in the context of two planar
model systems: an idealized nger model and a more anatomically-realistic leg model.
Figure 2.1 shows a simplied set of multi-articular muscles in both contexts of the nger
and leg. For notational convenience, we consider the general case of a planar limb of n
joints and let l
1
;:::;l
n
be the lengths of the individual limb segments. The limb posture
9
is dened by the column vector of joint angles = (
1
;:::;
n
)
T
where the superscript ()
T
denotes the transpose operator. The joints are driven by m muscles (m>n) of variable
stiness parametersk
i
,i = 1;:::;m, whose tendon paths produce annm moment arm
matrix R() (Valero-Cuevas 2009),
R =
0
B
B
B
B
B
B
@
r
11
r
12
r
1m
.
.
.
.
.
.
.
.
.
.
.
.
r
n1
r
n2
r
nm
1
C
C
C
C
C
C
A
: (2.1)
Each entry r
ij
(here, i = 1;:::;n, j = 1;:::;m) denotes the \rotation" the j
th
tendon
produces at the i
th
joint. It is positive when pulling the j
th
tendon induces a positive
rotation (i.e., counterclockwise rotation per the right-hand rule). Depictions of the nger
(n = 2) and leg (n = 3) models and their tendon paths are shown in Figure 2.1. Note
that the extensor mechanism of the ngers is rather complex (Valero-Cuevas, Yi, Brown,
McNamara, Paul & Lipson 2007), and Figure 2.1(b) corresponds to an idealized case
of multi-articular muscles whose moment arms change signs. The tendon paths of the
leg model in Figure 2.1(d) are more anatomically-realistic (Kuo & Zajac 1993, Kutch &
Valero-Cuevas 2012, Arnold, Ward, Lieber & Delp 2010). For clarity, Figure 2.1(d) shows
only seven of the fourteen muscles/muscle groups, amounting to a total of twenty-three
individual muscles that we use in this leg model.
In both the nger and leg models, we examine two types of joints: simple hinge joints
where the moment arm matrix R is constant for all values of , and non-circular joints
10
where the moment arm matrix is posture-dependent. In the nger model, the constant
moment arm matrix is chosen to be
R =
0
B
B
@
0:8 0:8 0:7 0:7
0:5 0:2 0:5 0:2
1
C
C
A
; (2.2)
in length units of cm. These values lie in the range of average moment arm values for
the index nger reported in (An, Ueba, Chao, Cooney & Linscheid 1983, Valero-Cuevas,
Zajac & Burgar 1998). For the posture-dependent moment arms, we dene joint ellipses
with semi-major and semi-minor axes in the range of measured index moment arms found
in (An et al. 1983). In the leg model, we utilize the posture-dependent moment arm values
of the model developed by (Arnold et al. 2010).
A change in limb posture corresponds to a rotation =
1

o
of the joints from a
reference limb posture
o
to a new limb posture
1
. A rotation fully determines the
length changes s = (s
1
;:::; s
m
)
T
of all tendons from their reference length. The
length changes of the tendons are dened { without loss of generality { at the reference
posture of the limb (An et al. 1983),
s =
Z

1
o
(R())
T
d =R
T
; (2.3)
where the negative sign indicates that a positive rotation of the joint will shorten the
tendons that induce it and vice versa.
For constant moment arms, the second equality in (2.3) is evident. In the case
of posture-dependent moment arms, because we only consider small posture changes
11
, we use the trapezoid integration rule such that equation (2.3) holds with R =
(R(
o
) + R(
1
))=2, where R(
o
) and R(
1
) are the moment arm matrices evaluated at
postures
o
and
1
, respectively.
We assume that each tendon is elastic and acts as a linear spring whose stiness
parameter is controlled to produce movement. We also allow each tendon to be pre-
stretched from its resting length by an amount l
o
= (l
1
;:::; l
m
)
T
. Basically, each
tendon has a baseline tension even at its reference length, and does not go slack at any
posture.
The total strain energy of the system, E, in the presence of tendon length changes
s is then given by
E =
1
2
(s + l
o
)
T
K(s + l
o
) :
(2.4)
As per the impedance (stiness) control formulation (Hogan 1984a), the stiness matrix
K is anmm diagonal matrix of entriesk
i
corresponding to the stiness of each muscle
i.
K =
0
B
B
B
B
B
B
@
k
1
::: 0
.
.
.
.
.
.
.
.
.
0 ::: k
m
1
C
C
C
C
C
C
A
: (2.5)
The associated muscle/tendon forces f = (f
1
;:::;f
m
)
T
are given by
f =K(s + l
o
) : (2.6)
12
These tendon forces, in turn, produce torques = (
1
;:::;
n
)
T
at the joints dened by
= Rf.
A stable posture is achieved by satisfying the static equilibrium condition at the joints,
namely, the total torque at each joint must be zero,
+() = 0 ; (2.7)
where () = (
1
;:::;
n
)
T
denotes externally-applied torques, such as gravity-induced
torques at the joints of larger limbs. Slow, i.e., quasi-static, limb movements are achieved
by sequentially satisfying this static equilibrium at each posture. It should be noted
that the equilibrium condition means that the strain energy E is at a local minimum.
That is to say, equation (2.7) is equivalent to being at a local minimum of the strain
energy paraboloid described by equation (2.4). However, a posture can be at equilibrium
at dierent energy levels. One has, at each posture, an admissible family of strain en-
ergy paraboloids, achievable at dierent combinations of spring stinessK and excursion
lengths s. Figure 2.2 illustrates this concept schematically.
In the absence of external loading (for = 0), and if the tendons are not pre-tensioned
(l
o
= 0), substituting (2.3) and (2.6) into Rf = 0 gives
RKR
T
= 0 ; (2.8)
which only admits the trivial solution = 0 for all K because the nn square matrix
RKR
T
is full rank when R is full row rank. Thus, there exists only one trivial solution
13
∆ θ
1
∆ θ
2
Strain energy
E
1
E
3
E
2
min E
k
Figure 2.2: A sketch of the static equilibrium and energy optimization problem: for each
sample posture, the static equilibrium condition can be satised at various combinations
of stiness values, corresponding to the minimum of each of many energy paraboloids.
The optimization problem is solved by choosing the stiness combination that minimizes
the strain energy. The green plane is an example of the set of possible postures, and the
blue surface represents a schematic illustration of the minimum strain energy manifold
over this range of postures.
at the reference posture, and other stable postures are not achievable by varying K. In
contrast, when the tendons are pre-tensioned, one has
RK(R
T
+ l
o
) = 0 : (2.9)
The system is thus unconditionally controllable and equilibrium postures are achievable
by proper choice of stiness K.
More generally, when the limb is subject to external load (non-zero ()), the equilib-
rium condition (2.7) takes the form RKR
T
( + l
o
) +() = 0: Here, in the absence
of muscle pre-tensioning (l
o
= 0), for each posture, i.e., for each value of , one has one
and only one solution where the internal tendon-driven torques equal the external
torques. This same posture can be achieved at various levels of muscle co-contraction, i.e.,
14
various values of K, satisfying the equilibrium condition. However, posture controllability
is limited to the solution of this equality. In contrast, when the tendons are pre-tensioned
(l
o
6= 0), the system is again unconditionally controllable by proper choice of stiness
parameters K. Since the external load () does not qualitatively aect the system con-
trollability, we hereafter set to zero to model a limb with no external load and low
inertia such as a nger.
It should be noted here that our work relates to the situation where muscle stiness
is arbitrarily controlled, where even if a joint rotation lengthens a muscle, the nervous
system can command that muscle to reduce or increase its stiness. Mathematically
speaking, for a given reference posture and pre-tensioning level l
o
, there exist several
combinations of stiness K that satisfy the static equilibrium condition. We choose the
pre-tensioning values l
o
to ensure that the tendons are nowhere slack in the workspace.
Thus, equilibrium postures can be computed without imposing additional restrictions
on l
o
.
Our quasi-static formulation implies that, starting from a given initial posture, the
limb has to transition to \nearby" equilibrium postures, that is to say, must be
small. A large value of means that the limb would instantaneously jump from its
current posture to a far away equilibrium posture, which violates the slow movement
assumption. Therefore, in all subsequent analyses, we consider quasi-static transitions to
only nearby postures. In particular, we formulate: (I) an optimality problem where for a
given equilibrium posture, we solve for optimal stiness values that minimize the strain
energy function associated with transitions to nearby postures, and (II) a reachability
problem where we explore reachable equilibrium postures for given stiness values.
15
Optimal equilibrium postures
Starting at a posture
o
, a desired limb posture near
o
can be achieved by tuning
the tendon stiness parameters so that the vector (R
T
+ l
o
) in equation (2.9) lies
in the null space of thenm matrix RK. The resulting stiness values are not unique, in
general. For each limb posture, equation (2.9) yields a family of solutions K, at multiple
strain energy levels, for which is achievable (Figure 2.2). This redundancy leads us to
look for optimal stiness values K
opt
that minimize the strain energy function while
satisfying the equilibrium constraint and inequality constraints on permissible stiness
values, namely,
min
K

E =
1
2
(R
T
+ l
o
)
T
K(R
T
+ l
o
)

subject to RK

R
T
+ l
o

+() = 0 ;
and k
min
k
i
k
max
;
k
min
> 0 :
(2.10)
where we set () = 0 and dene practical limits on the values of the stinesses k
i
. This
optimization problem is linear in the sense that the strain energy (cost) function and the
constraints are all linear in K. This suggests the possibility of degenerate minima if either
the strain energy or constraints functions are linearly dependent. We posit that this is
an unlikely design in biological and man-made systems; the numerical results presented
in this study are not degenerate.
16
Reachable equilibrium postures
Given a specic combination of stiness parameters K, the achievable equilibrium
posture can be found by solving the \forward problem" in (2.9), which we rewrite
as RKR
T
= RKl
o
. That is, given an initial posture
o
and muscle stiness K,
we compute a new posture =
o
+ that satises the equilibrium condition (2.7).
Given that the matrix RKR
T
is invertible, solutions to this forward problem are unique:
for each K, there exists one and only one , which in turn gives a uniquely-dened
new equilibrium posture . While the new posture is at equilibrium in the sense
that it corresponds to a local minimum of the strain energy paraboloid E, it does not
necessarily correspond to the minimum of all possible energy levels. That is to say, no
particular eort is made here to ensure optimality of the equilibrium postures. But as
in the optimality problem, we only examine small changes of equilibrium posture .
We systematically sample the stiness space at dierent precision levels, i.e., degrees
of discretization granularity k. For each precision level k, we calculate the set of
admissible departures from a reference posture, and consequently the set of reachable
(feasible) postures. The goal here is not to minimize the strain energy at each reachable
posture because the forward problem (computing given K) admits a unique solution.
Rather, our goal is to nd the set of all nearby reachable postures as a function of the
precision with which individual stinesses can be controlled.
Parameter values
A few remarks on the parameter values of k
min
and k
max
are in order here. k
min
is
always greater than zero ensuring that the optimization problem is well-posed. Here, we
17
initial pre-tensioning tendon slack length
lengthened
amount
new pre-tensioning
initial posture
after
posture change
new
initial posture
initial pre-tensioning tendon slack length
new pre-tensioning
shortened
amount
(a) (b) tendon lengthening tendon shortening
Figure 2.3: Variable pre-tensioning update rule: the pre-tensioning is history-dependent,
taking into account the tendon's lengthening (a) or shortening (b) that took place in the
previous posture change.
choose a range of stiness valuesk
min
= 100 N/m tok
max
= 1000 N/m, in agreement with
experimental data that measured stiness in a variety of muscles ranging from those that
actuate a nger to those of a leg (R atsep & Asser 2011, Chuang, Wu & Lin 2012, Bizzini
& Mannion 2003, Mustalampi, H

'akkinen, Kautiainen, Weir & Ylinen 2013). We now
non-dimensionalize all parameters using a characteristic length scale l

= 10 cm and a
characteristic force f

= 100 N. Non-dimensional parameters enable us to obtain generic
results, which we can later easily scale up or scale down to dierent limb sizes. Using
the characteristic length and force scales, the dimensionless range of stiness values is
k
min
= 0:1 to k
max
= 1 { making k = 0:05 equal to 5 % of the range. The normalized
values of the pre-stretched lengths l
i
are all set to 0:1 in the nger model and to 0:5 in
the leg model.. Note that in addition to constant pre-stretch, we consider pre-stretched
lengths that vary with posture.
Pre-tensioning of the tendons
We investigate the eect of constant versus variable tendon pre-tensioning. Constant
18
pre-tensioning means that the system is \memoryless" and resets to the initial l
o
because
pre-tensioning does not depend on the previous tendon lengths. Variable pre-tensioning
is dened according to the update rule depicted in Figure 2.3, where the pre-tensioning
value at each posture is equal to that at the previous posture plus the cumulative change
in tendon length that occurred as a result of moving from the previous to the current
posture; thus the pre-tensioning has \memory."
Workspace of the nger endpoint
We now focus on the nger model presented in Figure 2.1. Given = (
1
;
2
)
T
, the
nger's posture in the (x;y)-plane follows in a straightforward way
x
w
=l
1
cos
1
+l
2
cos(
2

1
); y
w
=l
1
sin
1
+l
2
sin(
2

1
):
(2.11)
The nger workspace W is the set of all points (x
w
;y
w
) of the (x;y)-plane accessible by
the nger's endpoint for all admissible
1
and
2
, A depiction of the joint angles space
and the workspace of the nger model is shown in Figure 2.4. Here, the admissible joint
angles are taken to lie in the range
1
2 [80

; 10

],
2
2 [0

; 90

]. Initial postures
o
in
this admissible range are discretized using a 9 9 regular grid, marked by. In the leg
model system, initial postures are taken in the range
1
2 [20

; 100

],
2
2 [0

; 100

] and

3
2 [20

; 30

], and twenty-seven equilibrium postures are computed in the vicinity of
each initial posture. For each
o
, we explore optimal transitions to equilibrium postures
in a small square neighborhood centered at
o
and of side length 10
o
. More specically,
we consider a total of nine nearby equilibrium postures, including staying at the initial
19
(a)
x
-0.2 0 0.2 0.4 0.6 0.8
y
-0.8
-0.6
-0.4
-0.2
0
θ
1
-80˚ -70˚ -60˚ -50˚ -40˚ -30˚ -20˚ -10˚ 0˚ 10˚
θ
2
-90˚
-80˚
-70˚
-60˚
-50˚
-40˚
-30˚
-20˚
-10˚
0˚
(b)
finger posture
finger
workspace
Figure 2.4: Finger model: (a) conguration space and (b) the corresponding endpoint
space which deforms the squares into "diamonds."
posture
o
. We compute the optimal stiness values for the change in joint angles
associated with a transition from the initial posture to each of these nine postures. We
linearly interpolate these values in the neighborhood of each initial posture to construct
a complete map of optimal stiness values.
Quasi-static trajectories
We consider quasi-static trajectories where the limb endpoint is required to slowly trace,
back and forth, curved and straight lines in the workspace as depicted in Figure 2.5.
We solve for the optimal stiness combinations required by the limb to perform this
back and forth motion. To this end, the quasi-static trajectories are discretized, and
the optimization problem in equation (2.10) is solved sequentially along the discrete
trajectory.
Reachable sets in the nger workspace
We discretize the range of stiness values fromk
min
tok
max
using constant increments
k, producingM = (k
max
k
min
)=k+1 discrete stiness values for each tendon. Form
20
Figure 2.5: Endpoint trajectories for (a) and (b) nger model and (c) leg model: nger
endpoint is required to trace (a) a curved trajectory, (b) a straight line, and the leg model
a curved trajectory (c). Subgures depict the prescribed values of (x
w
;y
w
) and associated
(
1
;
2
[;
3
]) and (
1
;
2
[;
3
]).
tendons, this discretization generates M
m
dierent stiness matrices K, which amounts
essentially to a uniform sampling of the m-dimensional stiness space. For example,
setting k = 0:05 in the nger model amounts to 19
4
distinct K matrices. For each
K, we solve equation (2.9) to compute starting from a reference posture
o
, and
use (2.11) to compute the reachable endpoint location (x
w
;y
w
). We constrain the full set
of M
m
reachable points to those that are in the vicinity of the reference posture.
21
2.3 Results
Quasi-static trajectories
Figure 2.6 depicts the optimal stiness values for the nger model tracing the trajecto-
ries of Figure 2.5(a,b) for all combinations of moment arms and pre-tensioning scenarios,
whereas Figure 2.7 shows similar results for the leg model tracing the curve of Fig-
ure 2.5(c). The optimal stiness values are characterized by a local jump as the limb
reverses its motion to trace the trajectory backward. The jump in the stiness of certain
tendons gets attenuated in the case of variable moment arms, but all jumps get smoothed
out completely when the pre-tensioning is changing continuously with posture. In these
cases, one observes the relaxation of the some tendons during the activity of the opposing
tendons, and the opposite eect is clearly seen during the return portion of the movement.
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
% trajectory
0 80
stiffness
20 40 60 100
0.1
0.2
0.3
0.4
0.5
0.6
0.7
% trajectory
0 80
stiffness
20 40 60 100
0.1
0.2
0.3
0.4
0.5
0.6
0.7
m
1
m
2
m
3
m
4
(b)
(a)
constant ∆l
0
and R posture-dependent R posture-dependent ∆l
0
posture-dependent ∆l
0
and R
m
1
m
2
m
3
m
4
Figure 2.6: Finger model: optimal stiness values for the quasi-static trajectories shown
in Figure 2.5(a,b).
22
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
% trajectory
0 80 20 40 60 100
0
0.2
0.4
0.6
0.8
1
% trajectory
0 80
stiffness
20 40 60 100
constant ∆l
0
and R variable R variable ∆l
0
variable ∆l
0
and R
m
7
m
8
m
17
m
19
m
21
m
22
remaining
muscles
Figure 2.7: Leg model: optimal stiness values for the quasi-static curved trajectory
shown in Figure 2.5(c). The most active muscles during this trajectory are the gluteus
maximus (m
7
and m
8
), rectus femoris (m
17
), biceps femoris long head (m
19
), tensor
fasciae late (m
21
) and tibialis anterior (m
22
).
Optimal equilibrium postures over the whole workspace
Figure 2.8 depicts the optimal stiness values K
opt
for the nger model with constant
moment arms R and constant pre-tensioning l
o
over the whole workspace. Note that
the workspace is discretized as depicted in Figure 2.4 and described in the previous
section. Tendons m
2
and m
4
exhibit higher stiness values than m
1
or m
3
A closer look
at the stiness variation within each neighborhood shows that, in the context of constant
moment arms, all neighborhoods are identical, that is to say, optimal stiness values
depend only on the change in joint angles , or posture change, and not on the initial
posture. Further, the stiness variation within each neighborhood is such that m
2
and
m
4
counter-balance each other in the
1
-direction { as the stiness of m
2
increases, that
of m
4
decreases and vice-versa.
Figure 2.9 depicts the optimal stiness values for the case when the moment arms vary
with joint angles. As in the case of constant moment arms, m
2
and m
4
require higher
23
stiness values, and the neighborhood around each initial posture indicates that each
tendon has a \preferred" direction of higher stiness, which is the same as in Figure 2.8.
However, the optimal stiness values now depend on moment arms, and thus on initial
posture. Interestingly, the stiness variation across the whole workspace is also charac-
terized bym
2
andm
4
counter-balancing each other, indicating an obligatory co-variation
in the way this pair of muscles acts in transitions to nearby equilibrium postures and to
postures across the whole workspace.
θ
1
-70° -50° -30° -10° 10°
θ
2
-90°
-70°
-50°
-30°
-10°
θ
1
-70° -50° -30° -10° 10°
θ
1
-70° -50° -30° -10° 10°
θ
1
-70° -50° -30° -10° 10°
x
-0.2 0 0.2 0.4 0.6 0.8
y
-0.8
-0.6
-0.4
-0.2
0
x
-0.2 0 0.2 0.4 0.6 0.8
x
-0.2 0 0.2 0.4 0.6 0.8
x
-0.2 0 0.2 0.4 0.6 0.8
stiffness stiffness
m
1
m
2
m
3
m
4
(a)
(b)
0.12
0.16
0.20
0.24
0.28
0.10
0.32
0.12
0.16
0.20
0.24
0.28
0.10
0.32
Figure 2.8: Finger model with constant moment arms: optimal stiness parameters
needed to achieve each posture in (a) conguration (joint angle) space, and (b) end-
point space. The plots show the stiness levels of each tendon m
i
, i = 1;:::; 4, that
minimize the strain energy function.
24
θ
1
-70° -50° -30° -10° 10°
θ
2
-90°
-70°
-50°
-30°
-10°
θ
1
-70° -50° -30° -10° 10°
θ
1
-70° -50° -30° -10° 10°
θ
1
-70° -50° -30° -10° 10°
x
-0.2 0 0.2 0.4 0.6 0.8
y
-0.8
-0.6
-0.4
-0.2
0
x
-0.2 0 0.2 0.4 0.6 0.8
x
-0.2 0 0.2 0.4 0.6 0.8
x
-0.2 0 0.2 0.4 0.6 0.8
stiffness
0.1
0.3
0.5
0.7
0.9
stiffness
m
1
m
2
m
3
m
4
(a)
(b)
0.1
0.3
0.5
0.7
0.9
Figure 2.9: Finger model with posture-dependent moment arms.
Figure 2.10(a) shows the optimal stiness combinations in joint angle space for con-
stant moment arms, while Figure 2.10(b) for variable, posture-dependent moment arms.
For clarity, only the leg muscles that manifest the highest stiness values are shown,
namely, muscles gluteus maximus, rectus femoris, tensor fasciae latae and tibialis ante-
rior, labeled m
8
, m
17
, m
21
and m
22
, respectively. The remaining muscles, while they
individually have low stiness values, contribute collectively to the equilibrium postures.
As in the case of the nger model, optimal stiness combination K
opt
depends on the
change but is independent of the initial joint angles in the case of constant moment
arms but depends on both in the case of posture-dependent moment arms.
25
100
80
60
θ
2
(knee)
40
20
100
80
60
θ
1
(hip)
40
20
0
20
0
-20
θ
3
(ankle)
m
8
100
80
60
40
20
100
80
60
40
20
0
m
17
100
80
60
40
20
100
80
60
40
20
0
m
21
100
80
60
40
20
100
80
60
40
20
0
m
22
0.1
0.3
0.5
0.7
0.9
stiffness
100
80
60
40
20
100
80
60
40
20
0
m
17
100
80
60
40
20
100
80
60
40
20
0
m
21
100
80
60
40
20
100
80
60
40
20
0
m
22
100
80
60
θ
2
(knee)
40
20
100
80
60
θ
1
(hip)
40
20
0
20
0
-20
θ
3
(ankle)
m
8
stiffness
0.1
0.2
0.3
0.4
0.5
0.6
(a)
(b)
Figure 2.10: Leg model: optimal stiness parameters needed to achieve each posture with
(a) constant and (b) posture-dependent moment arms. The plots show the activation level
of the four highly active muscles,m
8
m
17
,m
21
andm
22
(gluteus maximus, rectus femoris,
tensor fasciae latae and tibialis anterior, respectively), that minimize the strain energy
function.
Reachable equilibrium postures
We explore the reachability or feasibility of tracing certain trajectories within the limb's
workspace given constraints on the precision in the control of stiness. For illustration
purposes, we investigate the ability of the nger model to trace multiple trajectories for a
given k. Starting from the initial posture, we sequentially transition to the next point
along the trajectory by locating the nearest admissible reachable posture. Figure 2.11
shows the desired versus reachable endpoint trajectories for three precision levels k =
0:1; 0:08; 0:05, superimposed on the set of all reachable points in the limb workspace. The
reachable paths show discontinuities. These discontinuities decrease as the resolution
of stiness increases, going from k = 0:1 and k = 0:05, and eventually disappear
26
x
-0.2 0.4 0.6 0.8
y
-0.8
-0.6
-0.4
-0.2
0
0 0.2
achieved
desired
scaled 30% from original - 30 points, 0.15 pret, then 110%
∆ k = 0.1 ∆ k = 0.08 ∆ k = 0.05
x
-0.2 0.4 0.6 0.8 0 0.2
achieved
desired
x
-0.2 0.4 0.6 0.8 0 0.2
achieved
desired
Figure 2.11: Finger model: achievable versus desired trajectories overlaid on top of the
reachable postures. The level of stiness precision aects the reachability. As precision
increases (from left to right), the endpoint trajectories become smoother and closer to
the desired trajectories. Discontinuities in the set of reachable points vary nonlinearly in
the workspace, e.g., regions highlighted in red.
as k is varied continuously. The discontinuous behavior arises from unreachable nearby
desired postures along the discretized path, in which case the limb settles at the nearest
admissible posture from the set of reachable nearby postures. Further, we observe that
discontinuities exhibit a nonlinear structure in the limb workspace, implying that the
demands on the precision in stiness varies nonlinearly in space.
2.4 Discussion
We developed a novel mathematical formulation for slow and accurate tendon-driven
limb movements using the framework of driftless mechanical systems where movement is
generated by controlling the tendons' stiness parameters. We used this framework to
address three questions in slow limb movements: (1) can slow limb movements be achieved
by modulating the stiness parameters only? (2) what is the role of pre-tensioning? and
(3) how are these movements aected by limits on the precision with which stiness can
be controlled? We demonstrated that stiness control is sucient to produce accurate
27
and slow limb movements, but only when the tendons are pre-stretched. We then probed
the limb's ability to trace prescribed trajectories. When the moment arm R or muscle
pre-tensioning l
o
are held constant, the stiness values required to track the prescribed
trajectories are characterized by sharp \kinks" or discontinuities, which get smoothed out
for posture-dependent R and l
o
.
We analyzed the reachability of the limb's workspace by exploring the limb's ability
to trace prescribed trajectories under constraints on stiness precision. Here, we did
not impose any restriction on the strain energy level. By comparing dierent levels of
precision (i.e., resolution) to control stiness, we found that high resolution of stiness
is key for the limb to improve its accuracy in tracking a trajectory. More specically,
we identied discontinuous patterns in the limb workspace at low muscle precision that
almost disappeared at higher precision. Our results also revealed that the required sti-
ness resolution is not uniform throughout the limb's workspace in the sense that lower
muscle resolutions induce discontinuities that are non-uniform and are localized to certain
regions of the workspace.
In addition to their direct relevance to understanding the necessary conditions and
limitations of bio-inspired, tendon-driven robotic limbs, our results bear relevance to
the neural control of movement. In that literature, slow limb movements are, to our
knowledge, not addressed. Limb movements in general have at times been considered
from the perspectives of impedance control, force-eld primitives, spinal feedback sys-
tems (Giszter, Patil & Hart 2007, Giszter 2015), or the framework of the equilibrium
point hypothesis (EPH), which states that limb movements are generated by a sequence
of equilibrium points along a desired trajectory (Asatryan & Feldman 1965, Feldman &
28
Levin 2009, Cooke 1980). Our mathematical formulation is reminiscent of but not identi-
cal to these approaches. EPH, for example, considers that the magnitude of the \control"
force exerted on the limb, at any time, depends on the dierence between the \actual" limb
dynamics and the desired equilibrium point, and derives a second-order dierential equa-
tion for the \error" between the actual and desired dynamics (Shadmehr 1998, Asatryan
& Feldman 1965, Feldman 1966). Here, we were not concerned with the relaxation dy-
namics. Our working hypothesis is that this relaxation time scale is much faster than
the time scale associated with slow limb movements, and the limb reaches a state of
equilibrium instantaneously in response to changes in muscle stiness. Our formulation
focused on determining whether equilibrium postures of multi-joint limbs and quasi-static
transitions among them can be achieved by direct control of muscle stiness values, and
the eect of stiness resolution on achievable postures.
The direct control of muscle stinesses can be viewed as an abstraction, or `meta-
model', of muscle control. While the fact that muscles have stiness (and strain energy)
{ and that the nervous system can regulate them { is well established (e.g., Mussa-Ivaldi
& Hogan (1991),Burdet, Osu, Franklin, Yoshioka, Milner & Kawato (2000),Burdet, Osu,
Franklin, Milner & Kawato (2001)), our main goal was to present a general mathemati-
cal framework for analyzing slow limb movements that arise from the control of muscle
stiness, rather than to focus on the details of the production and control of muscle sti-
ness. In addition, our use of high and low-resolution regulation of stiness is not meant
to indicate a discontinuity in the time history of the neural control signal. Rather, it is
simply a means to indicate a lack of accuracy in the production of a given level of muscle
29
stiness at a given posture. Signal dependent noise and motor unit recruitment and rate
coding are well known sources of such inaccuracies and high-frequency variability.
Experimental work on muscle force has shown that if a muscle is clamped at a con-
stant length and stimulated (an isometric contraction), its active force depends on the
length in a non-monotonic way, with lower force at short or long lengths and a maximum
force at some length in between. This eect is called the force-length relationship. Im-
portantly, the force-length relationship is not like a Hook's law for muscle (i.e., a stress
vs. strain relationship). It simply represents the amount of active isometric contrac-
tile force a muscle can produce at a given lengthbut not directly its passive stiness.
In addition, if a muscle is allowed to shorten (a concentric contraction), the force will
also depend on the shortening speed, decreasing as the shortening speed gets higher.
If the muscle is lengthened by some external force, even as it produces force to resist
the lengthening (an eccentric contraction), the active force increases above the isomet-
ric value (McMahon 1984, Chen, Tian, Iwasaki & Friesen 2011). Together, these ef-
fects are called the force-velocity relationship and can be modeled succinctly as done in
Hill (Hill 1938) and Williams (Williams 2010) using a mathematical product of these two
eects and an activation function. Given that our study is concerned with slow limb
movements, one can reasonably argue that the shortening velocity is slow and can be
considered a constant. An interesting outcome of our formulation is that, although we
use a spring model for the muscles, the muscle forces required to trace a trajectory are not
necessarily linear. Indeed, a re-drawing of the muscle stinesses k
i
reported in Figure 2.6
as functions of the joint angles and thus of muscle lengths would reveal a highly nonlinear
30
dependence of k
i
on the change in muscle length s
i
. A detailed analysis of this non-
linearity, as well as how to incorporate the Hill-type force model into our formulation,
will be the topic of future work.
Another future direction is to investigate the similarities between the discontinuities
observed in Figure 2.11 and the discontinuous nature of slow nger movements reported
in numerous experimental studies (Vallbo & Wessberg 1993, Darling, Cole & Miller 1994).
Several studies attributed these discontinuities to neural sources while others held respon-
sible peripheral stretch re
exes (Gross, Timmermann, Kujala, Dirks, Schmitz, Salmelin
& Schnitzler 2002, Williams, Soteropoulos & Baker 2009, Evans & Baker 2003). Our
formulation oers a novel framework that can be used in conjunction with further exper-
imentation to decipher the origins of these discontinuities and their dependence on the
nature and integrity of the neural control of musculature in health and disease.
31
Chapter 3
Experimental evidence of non-smoothness in nger motion
3.1 Introduction
Our previous mathematical modeling of tendon-driven limbs suggested that low resolu-
tion of muscle stiness control may contribute to discontinuities in slow limb movements
(Figure 2.11). In this experimental study, we focus on the kinematics of slow nger
movements and investigate a type of nger tasks that is known to manifest oscillations
in healthy adults: slow nger movements. Slow nger movements are characterized by
discontinuities in trajectory, within a dominant frequency range of 6-13 Hz (Vallbo &
Wessberg 1993, Wessberg & Vallbo 1995, Wessberg & Vallbo 1996). While the frequency
content of discontinuities during slow nger movements is well known, whether or not
the locations of the discontinuities are locked to joint angles during the trajectory has
not been explored. The aim of this study is to determine if there exist location-specic
discontinuities during slow nger movements in a planar trajectory.
32
3.2 Methods
Overview
We instructed eight participants (mean age 25:54 years, 6 male) to rest their right arm
on a table. The proximal and distal interphalangeal joints (PIP and DIP, respectively) of
the index nger were splinted with a wooden stick. The remaining digits held a vertical
handle xed to the table, so only the metacarpophalangeal (MCP) joint of the index
nger was free to move in
exion-extension on the horizontal plane (Figure 3.1).
Figure 3.1: The setup showing a 3-axis accelerometer attached at the tip of a wooden
stick, and ve re
ective markers on the index nger and hand: ngertip, distal inter-
phalangeal (DIP), proximal interphalangeal (PIP), metacarpophalangeal (MCP) and the
second carpometacarpal (CMC) joints.
33
Acceleration recordings
We measured ngertip accelerations using a tri-axial accelerometer (LIS331AL, STMi-
croelectronics, Switzerland) attached to the nail of the index nger, at a sampling rate
of 1000 Hz. The accelerometer was positioned such that the z-direction captured the
acceleration tangential to the movement. We simultaneously recorded 3D joint positions
at 200 Hz by a motion capture system (Vicon Motion Systems Ltd. UK) comprising of
four infrared cameras and ve re
ective markers placed on the tip of the index nger,
the proximal and distal interphalangeal (PIP and DIP, respectively), the metacarpopha-
langeal (MCP) and the second carpometacarpal (CMC) joints (Figure 3.1).
Finger motion task Participants performed 30-second trials consisting of repetitions
of
exion and extension motions at the MCP joint by following a metronome (one
exion
or extension per beat). We instructed them to perform quasi-sinusoidal
exion-extensions
(rather than brisk movements between beats), from maximum
exion to maximum ex-
tension: a range of motion of around 70 degrees. We randomly chose the speed of the
metronome in each trial from very slow, slow, medium to fast: 12, 30, 50 or 70 beats per
minute (BPM), respectively. We recorded a total of 100
exions and 100 extensions for
each speed.
Analysis
We used the position data from the motion capture system to calculate the MCP joint
angles, and then to separate the position and acceleration data into
exion and extension
trials. We band-pass ltered (5 Hz to 50 Hz) the accelerations and mapped them to
the calculated MCP joint angles. We then separated them into eight conditions (four
34
speeds and two movement directions). We discarded about 15 degrees of motion at the
very beginning and very end of the movements to remove inertial eects and changes in
velocity during motion reversal.
For each of the eight conditions, we averaged the accelerations across joint angles
throughout the trajectory of interest. This was done by locating the acceleration value at
each joint angle and averaging the accelerations point-by-point over all repeated trials.
Power Spectral Density
The power spectral density (PSD) is a representation of the power of a signal at each
frequency. It is the Fourier transform of the autocorrelation of a signal. PSD is useful for
identifying oscillatory signals in time series data and knowing their amplitude. A higher
value of power at a specic frequency or range of frequencies indicates that variations or
oscillations at those frequencies are strong. We computed the PSD of motion trials in an
attempt to compare their frequency contents.
3.3 Results
We dene a discontinuity in the trajectory as an angular position where the ngertip
acceleration exhibits a sudden increase, which would indicate that the nger jumped over
that angle rapidly (Vallbo & Wessberg 1993). Figure 3.2 shows the averaged accelerations
over approximately 80 trials for each of eight conditions for one participant. There exist
structures in the acceleration which we classify as discontinuities. For example, there
is a discontinuity at 160 degrees, as indicated by a large amplitude in acceleration that
persists throughout speeds.
35
130 140 150 160 170 180 190
acceleration
-0.2
0
0.2
130 140 150 160 170 180 190
acceleration
-0.4
-0.2
0
0.2
0.4
130 140 150 160 170 180 190
acceleration
-0.2
0
0.2
130 140 150 160 170 180 190
acceleration
0
0.2
0.4
130 140 150 160 170 180 190
acceleration
-0.2
0
0.2
130 140 150 160 170 180 190
acceleration
-0.2
0
0.2
angle (deg)
130 140 150 160 170 180 190
acceleration
-0.2
0
0.2
angle (deg)
130 140 150 160 170 180 190
acceleration
-0.2
0
0.2
12 BPM 12 BPM
30 BPM 30 BPM
50 BPM 50 BPM
70 BPM 70 BPM
extension flexion
Figure 3.2: Averaged accelerations across all
exion and extension trials for one subject,
for slow, medium and fast speeds (an average of 80 trials per speed). Maximal extension
is represented by higher angles, and maximal
exion by lower angles.
Figure 3.3 shows the power spectral densities of the mean acceleration computed for
all eight conditions for the same subject as in Figure 3.2. The slower
exion speed (12
BPM) contains higher frequencies at 28-30 Hz, which disappear at higher speeds.
36
0 10 20 30 40 50
0
0.2
0.4
0.6
0.8
1
0 10 20 30 40 50
0
0.2
0.4
0.6
0.8
1
extension
flexion
12 BPM
30 BPM
50 BPM
70 BPM
12 BPM
30 BPM
50 BPM
70 BPM
Frequency (Hz)
Power Power
Figure 3.3: Power spectral densities of motion trials for all four speeds and two types of
motion (normalized across all speeds).
3.4 Discussion
In this study, we investigated the spatial characteristics and frequency content of slow
nger movement trajectories. We found dierences in accelerations at specic angles,
from which discontinuities could be determined. While some discontinuities persisted
throughout speeds, others either vanished at higher speeds or combined into one larger
discontinuity (Figure 3.2). Furthermore, our power spectral densities showed that slower
movements contain higher frequencies. This was evident in the slowest trials of Figure 3.3,
where higher frequencies of acceleration were observed. These frequencies matched those
seen in physiologic tremor of the nger and became less apparent in faster trials. We
speculate that this occurs because motor neuron signals during slow movements do not
37
reach the minimum threshold that produces muscle contraction, rather, they oscillate
around the thereshold resulting in discontinuous and uncoordinated activity in muscles.
38
Chapter 4
Goal-equivalent nger force variation in EEG activity
4.1 Introduction
A key feature of human movement is that task goals can be repeatedly achieved even
though the ne details of the movement are never repeated exactly (Bernstein 1967).
Increasing evidence suggests that the lack of detailed repeatability in movement (inter-
changeably called goal-equivalent or task-irrelevant variation), rather than being noise,
is actually an intended feature that may promote health by distributing tissue loading
and/or minimizing muscle fatigue (Hamill et al. 2012, Lipsitz 2002). Goal-equivalent
variation has been referred to as good variation and contrasted with bad variation, which
would lead to motor errors not consistent with the goal of the task (Latash & Anson 2006).
Goal-equivalent variation has been observed in many types of movements and tasks.
Such variation has been studied using Uncontrolled Manifold (UCM) analysis in such
varied tasks as multi-nger force production, multi-joint reaching, standing, swaying and
stepping (see Latash (2008) for review). These studies mostly focus on comparing vari-
ance across task repetitions in goal-equivalent and non-goal equivalent directions at the
39
level of motor output (e.g. kinematics or endpoint force). Goal-equivalent variation has
been shown using electromyography (EMG) to be present at the level of muscle com-
mands (Krishnamoorthy et al. 2003). However, the level in the nervous system at which
goal-equivalent variation appears has not been established. A subcortical-level hypoth-
esis for goal-equivalent variation is that descending motor commands from the cortex
are constant across repetitions of the same task, and goal-equivalent variation across
repetitions happens at the level of subcortical or spinal circuitry that coordinates mus-
cle activity (Giszter, Mussa-Ivaldi & Bizzi 1993, Saltiel, Wyler-Duda, D'Avella, Tresch
& Bizzi 2001). An alternative cortical-level hypothesis is that goal-equivalent variation
is directly formed at the level of the cortex and is present in descending motor com-
mands (Latash & Huang 2015).
Here we take the novel step of testing these hypotheses and thus localizing goal-
equivalent variation in the nervous system; we specically aim to determine if goal-
equivalent variability can be predicted by cortical-level signals. We made synchronous
recordings of bi-lateral nger force and whole-scalp electroencephalographic (EEG) signals
while participants repeatedly performed a task that allowed for goal-equivalent variation.
We reasoned that if we found cortical-level signals that could predict goal-equivalent vari-
ation prior to the start of each repetition, we could reject the subcortical-level hypothesis
that descending cortical signals were constant and support the cortical-level hypothesis
for goal-equivalent variation in motor output.
40
4.2 Methods
Participant Population
We recruited 18 healthy right-handed participants (9 males and 9 females) with a mean
age 27.4 4.5 years (range, 23 41). The studies we describe here were performed at the
University of Southern California and approved by the University of Southern California
Institutional Review Board. All participants provided informed consent.
Overview of experiments
We used a bilateral nger task where participants applied simultaneous left and right
index nger forces to repeatedly match the sum of forces to a target feedback force. Par-
ticipants always matched the target, but there was a variation in left and right forces
across repetitions: some repetitions were achieved with higher left force than right, and
vice-versa. We selected this task because it has been well-studied in the literature on goal-
equivalent motor output variation (Latash, Li & Zatsiorsky 1998, Latash, Scholz, Danion
& Sch oner 2001, Latash, Yarrow & Rothwell 2003). We used electroencephalography
(EEG) to measure cortical signals during repetitions of the task, and employed a machine
learning technique to nd out if EEG signals just prior to the start of each repetition were
correlated with goal-equivalent force output variations. We selected EEG to leverage prior
studies that have developed algorithms to decode the intention of explicit left vs. right
movements (Blankertz, Tomioka, Lemm, Kawanabe & Muller 2008, Pfurtscheller, Neu-
per, Flotzinger & Pregenzer 1997). Our analysis was based on classifying goal-equivalent
force output variation (two classes, left>right or right>left) using EEG signals only, and
41
comparing the EEG-based classication with the force-based classes. A successful classi-
cation would imply that brain signals associated with each class are in fact distinguishable,
and support the brain-level hypothesis for goal-equivalent motor output variation.
Finger force task
Participants performed a bimanual nger force task that involved abducting both index
ngers simultaneously against load cells (Measurement Specialties, Hampton, VA) (Figure
4.1a). We asked participants to repeatedly use both ngers to match a target total
nger force (left force [F
L
] + right force [F
R
]) of 6 Newtons (N). For simplicity in this
study, we chose a single target force level that would be easy to achieve for all healthy
young participants, but not so low that it would be dicult to perform reliably in many
repetitions: in a preliminary set of three participants (1 female, 2 males, mean age 31),
we identied that 6N was an appropriate force level, and used this force level in the main
study cohort of 18 participants, none of which reported discomfort or fatigue. Participants
received continuous visual feedback of their current total nger force as well as the xed
6N target; however, individual nger force was never provided to the participant as
feedback. Each repetition of the task included the following segments chosen based on
previous literature (Blankertz, Dornhege, Lemm, Krauledat, Curio & M uller 2006, M uller-
Gerking, Pfurtscheller & Flyvbjerg 1999): 1) 2 seconds of rest where nothing was shown
on the feedback screen, 2) 2 seconds of a preparation period in which participants were
cued to prepare for the task, and 3) 5 seconds of force hold. Participants performed a
total of 100 repetitions.
42
We show an example of one participant's force repetitions in Figure 4.1b. During
each repetition of the task, total force was kept at 6N (top trace), while individual
forces
uctuated. Some repetitions of the task were performed with more left nger
force (F
L
>F
R
), and some were performed with more right nger force (F
R
>F
L
). To
summarize the combination of forces used by the participant on each repetition of the task,
we computed separate averages for F
R
and F
L
across time starting one second after the
target was displayed (to allow for reaction time and force stabilization) and ending half a
second before the target disappeared (to exclude the return to baseline force). A total of
100 combinations ofF
R
versusF
L
along with a 95% covariance ellipse for one participant
are shown in Figure 4.1c. The major axis of the covariance ellipse represents the goal-
equivalent direction along which all combinations do not aect the task of producing 6N
of total nger force. The orthogonal direction is the non-goal-equivalent direction. To
identify repetitions of the task that were clearly dierent in the combination of F
R
and
F
L
used by the participant, we removed force repetitions belonging to the middle third
portion of the covariance ellipse (shown in black in Figure 4.1c), to discard right and
left forces that were almost equal. We identied the two remaining groups of force hold
repetitions as two classes, F
L
>F
R
(shown in blue) and F
R
>F
L
(shown in red).
Electroencephalography
Along with force signals described above, we simultaneously recorded 64-channel whole-
scalp electroencephalography (EEG) using a
exible cap according to the International
10-20 system (ANT Neuro, Enschede, The Netherlands). EEG signals were preprocessed
to remove artifacts related to eye blinks and muscle activity using independent component
43
0 2 4 6
0
2
4
6
0
6
F
R
+ F
L
(N)
0 50 100 150
Time (sec)
0
3
F
R
, F
L
(N)
Task:
...
...
a b c
F
R
(N)
F
L
(N)
F
R
+ F
L
= 6N
F
L
F
R
F
L
> F
R
F
R
> F
L
F
L
> F
R

F
R
> F
L
Figure 4.1: Experiment setup and force output. a. The participant places their hands

at on a table and repeatedly exerts left (F
L
) and right (F
R
) index nger forces simul-
taneously, while keeping the sum of the two constant at a target force of 6 Newtons (N).
Only the sum of forces is shown to the participant as feedback on a screen, along with the
target force level; b. example of a representative participant's repeated force holds. The
top trace shows the participant maintained the sum of the two forces at 6N, while the
bottom traces show the left (in blue) and right (in red) forces for the same repetitions.
Some repetitions of the task were accomplished with a relatively higher right nger force,
and others with a relatively higher left nger force. The interval shaded in gray (1 second
after the start and 0.5 second before the end of the hold period) is used to determine the
time average of individual forces in each repetition; c. average right versus left forces dur-
ing all task repetitions for the same participant with a 95% covariance ellipse representing
the variance of the data. The major axis represents the goal-equivalent direction, where
the sum of forces is always 6N, while the minor axis is the non-goal-equivalent direction
where the sum of forces would deviate from the target of 6N. The highest variation is
seen in the goal-equivalent direction, while the orthogonal non-goal-equivalent direction
has much lower variance. In our analysis we discard repetitions where forces belonged
to the middle third portion of the ellipse, shown in black, to only keep repetitions where
left and right forces were clearly dierent. Two classes of repetitions emerge: F
L
> F
R
(in blue) and F
R
>F
L
(in red).
analysis implemented in EEGLAB (www.sccn.ucsd.edu/eeglab). Subsequent analyses
were performed on a subset of 27 EEG electrodes that span the sensorimotor area of
the cortex, which is common practice in EEG-based classication studies (Blankertz
et al. 2008, M uller-Gerking et al. 1999, Ramoser, Muller-Gerking & Pfurtscheller 2000).
Unilateral left versus right hand movements have been successfully predicted based on
EEG signals from the 500 milliseconds (ms) prior to movement onset (Blankertz et al.
44
2006, M uller-Gerking et al. 1999), with EEG signals band-pass ltered from 7-30 Hz to
include both alpha and beta bands. We applied this same algorithm to predict whether
a given repetition of our bi-manual nger task belonged to the F
L
> F
R
class or the
F
R
>F
L
class using only the EEG signals during this 500 ms pre-movement period.
Machine learning for EEG analysis
Numerous studies have successfully predicted unilateral right hand versus left hand
movement intention from EEG signals using a machine learning technique combining
the signal processing method called Common Spatial Patterns (CSP) (Dornhege 2007,
Fukunaga 1990, Lemm, Blankertz, Curio & Muller 2005), and the classication method
linear discriminant analysis (LDA) (Blankertz et al. 2006, Blankertz, Kawanabe, Tomioka,
Hohlefeld, M uller & Nikulin 2007, Blankertz et al. 2008, Dornhege, Blankertz, Krauledat,
Losch, Curio & M uller 2006, Ramoser et al. 2000). Pre-movement EEG signals of hand
movement intention exhibit an event-related desynchronization, i.e., a decrease in power,
of alpha and/or beta rhythms on the contralateral hand motor area. That is, the power (or
equivalently, the variance) of the band-passed EEG signals is lower on the contralateral
side, and consequently relatively higher on the ipsilateral side. Therefore the features
used in the classication of left versus right hand movement intention are comprised of
the EEG variances (powers) at the electrodes.
Following this technique, we used only EEG signals to predict whether each repetition
of the task belongs to the class F
L
> F
R
or F
R
> F
L
. For each participant, we used
leave-one-out cross-validation (LOOCV) as an initial test of classier validity, forming
a training set on all but one task repetition, testing on the remaining repetition, and
45
repeating the process so that all task repetitions are held out. The CSP technique uses
EEG signals from the training repetitions, and determines spatial lters (we used two
spatial lters per class as in prior work (Blankertz et al. 2008)) such that the ltered
EEG signals from the two classes are optimally distinguishable based on EEG variance.
Then, the class of the held out test repetition, whose EEG signals are projected with
the above lters, is predicted from EEG data only using an LDA classier function that
outputs a numerical value associated with either one of the classes. Spatial patterns from
each class can be visualized on the scalp and represent a measure of brain activation
during repetitions of the task belonging to each class.
Overview of CSP and LDA
EEG data from a single task repetition is encoded in a CT matrix X = [x(t);x(t +
1);;x(t +T 1)] by concatenating C 1 vectors x(t) which are the multi-channel EEG
data at a specic time point t, T is the total number of time points in one repetition,
and C is the number of EEG channels. The CSP algorithm is performed on Cn
1
T
matrixS
1
andCn
2
T matrixS
2
formed by concatenating alln
1
training examples of the
F
L
>F
R
class, and all n
2
training examples of the F
R
>F
L
class, respectively. Matrices
S
1
andS
2
are used in the CSP algorithm to determine a single matrixW (whose columns
are spatial lters) such that W
T
S
1
and W
T
S
2
are uncorrelated across channels, which
leads to easier distinguishing of the two classes. W is arranged such that the columns
are ordered according to eigenvalue; lters on the left end ofW correspond to those most
important for predicting theF
L
>F
R
class, and lters on the right end ofW correspond
to those most important for predicting the F
R
> F
L
class. W
T
S
1
and W
T
S
2
are then
46
used to calculate linear weights
j
and bias
0
of the LDA classier using Fisher's Linear
Discriminant Analysis (Hastie et al. 2005). A test sample X from a repetition not used
in training can then be classied using LDA:
f =
J
X
(j=1)

j
log(w
T
j
XX
T
w
j
) +
0
; (4.1)
wherew
j
are the lters determined earlier with CSP, and J is the total number of lters
used. We follow common practice of feature extraction in EEG classication studies and
extract the two most important spatial lters from either end of the W matrix (a total
of four) (Blankertz et al. 2008) to be used in the projection of signals in X; therefore, in
our case J = 4.
Appendix A provides a more detailed explanation of the CSP and LDA methods;
along with their mathematical derivations and two-dimensional toy examples that help
visualize how these learning and prediction algorithms are performed.
In contrast to prior studies that use EEG to construct brain-machine interfaces (BMI)
and assess single-subject performance by quantifying classication accuracy, we wanted
simply to ensure that there was a statistically signicant relationship (across multiple
participants) between the EEG-based classier output and the force-based class for the
repetition, regardless of classication accuracy. We therefore used a linear mixed-eects
(LME) model to quantify the association between force-based repetition class and EEG-
based classier output. The LME model of the EEG-based classier output (y) contained
a xed-eects term for force-based repetition class (x
i
;i = 1 for F
L
> F
R
and 2 for
47
F
L
>F
R
), as well as for participant j a random intercept (u
0j
) and random slope (u
1j
)
to capture participant- specic dependencies:
y
j
=
0
+
1
x
ij
+u
0j
+u
1j
x
ij
+
ij
: (4.2)
We examined the p-value for the xed-eect of force-based repetition class on EEG-
based classier output, and considered p< 0:05 as evidence for a statistically signicant
relationship between cortical activity and goal-equivalent variation in nger forces.
In addition to the LME model, we implemented a receiver operating characteristic
(ROC) analysis of test data in order to assess the performance of the classiers at various
threshold settings. For each participant, we created an ROC curve using all leave-one-out
classier outputs for that participant, and estimated the area under the curve (AUC) as
a measure of classier performance. We then used all participants' classier outputs to
plot the overall ROC curve and obtain an estimate of the overall AUC.
4.3 Results
We found that all 18 participants achieved the task with a signicantly higher variance of
repeated forces along the goal-equivalent direction than the non-goal-equivalent direction.
In Figure 4.2a, we show the 95% covariance ellipses of all 18 participants, along with the
goal-equivalent line. Some participants have more elongated ellipses than others, but
the relative variance along the goal-equivalent direction is always greater than in the
orthogonal direction, as illustrated in Figure 4.2b. The ratios of goal-equivalent variance
48
0 2 4 6
0
2
4
6
0.5
0.6
0.7
0.8
0.9
1
Variance Ratio
****
F
R
(N)
F
L
(N)
b a
Figure 4.2: Covariance ellipses and force variance ratio of 18 participants. a. 95%
covariance ellipses of all participants. The covariance ellipse of each participant was
constructed using the average left and right forces from all repetitions of the task (see
Figure 4.1c). The overlapped ellipses show that all participants maintained a relatively
higher variation in the goal-equivalent direction than in the non-goal-equivalent direction;
b. boxplot of the ratio of goal-equivalent force variance to the sum of goal-equivalent and
non-goal-equivalent variance. The ratio was signicantly greater than 0.5 (p < 0:0001,
****).
to the sum of both non-goal-equivalent and goal-equivalent variances in all participants
for the chosen force level of 6N were signicantly greater than 0.5 (p< 0:0001).
Our LME model indicated that force-based repetition class had a signicant eect on
EEG-based classier output (p< 0:0001). Figure 4.3a illustrates the EEG-based classier
output of each participant averaged across repetitions from each class, compared with
actual force-based classes. The x-axis represents the label of force-based classes, while
the y-axis is the LDA classier outputf, a linear combination of the natural log of channel
variances of spatially-ltered EEG data. The EEG-based classier was trained to produce
lower values for training repetitions for which F
L
>F
R
compared to training repetitions
49
for which F
R
> F
L
. Our results showed that test repetitions for which F
L
> F
R
were
consistently associated with lower EEG-based classier output than for F
R
> F
L
. For
each participant, a positive slope from F
L
> F
R
to F
R
> F
L
indicates a successful
prediction of the two classes.
The ROC analysis indicated an overall AUC of 0:72 across all participants. Fig-
ure 4.4 shows the individual ROC curves of cross-validated test data for each participant
(individual AUC values reported in Table 4.1), as well as an overall ROC curve using
all classier outputs from all participants. The classier output from each leave-one-out
cross-validation run of each participant was rst normalized using means and standard
deviations of training data classications; to ensure classier outputs were directly com-
parable across leave-one-out classiers and across participants. In a subsequent sensitivity
analysis, we removed repetitions that belonged to only half of the middle third portion of
the covariance ellipse in Figure 4.1c, and obtained an overall AUC of 0.68, slightly lower
than our reported value for the current analysis.
To visualize patterns of brain activations during each class of repetitions, we mapped
average spatial patterns, CSPs, from each class onto the 27 electrode locations on the
scalp (Figure 4.3b). The top and bottom rows represent the most important and second
most important patterns associated with the most important and second most important
spatial lters, respectively, that were used to predict EEG-based classes. Left and right
columns represent the patterns for the two classes, F
L
>F
R
and F
R
>F
L
, respectively.
The CSP patterns illustrate the amount of EEG variance at each electrode, which played
an important role in the discrimination between classes. Both the most important and
second most important patterns displayed laterality on the sensorimotor areas, i.e., a
50
Force-based class
Predicts F
R
> F
L
Predicts F
L
> F
R
a
F
L
> F
R
F
R
> F
L
F
L
> F
R
F
R
> F
L
0 0.1 0.2
EEG-based classifier output
b
Most important CSP
Second most important CSP
0 0.1 0.2
****
Figure 4.3: EEG-based prediction of force output and averaged common spatial patterns
(CSP). a. EEG-predicted classes versus force-based classes of all 18 participants. For each
participant, the classier outputs of repetitions from each class (F
L
>F
R
and F
R
>F
L
)
were averaged and joined by a line. Most participants show a clear distinction between
the two classes as indicated by a positive slope of the line. Linear mixed-eect modeling
demonstrated a signicant eect of force-based repetition class on EEG-based classier
output (p< 0:0001, ****); b. the most important (top row) and second most important
(bottom row) common spatial patterns (CSP) used for the EEG-based classication of
repetitions whereF
L
>F
R
(left) and repetitions where F
R
>F
L
(right), averaged across
18 participants. Enhanced EEG amplitudes are seen in the left hand motor area when
F
L
leads (top left). Similarly, EEG amplitudes are higher in the right hand motor area
when F
R
leads.
relatively higher EEG variance on the left hemisphere during repetitions belonging to
class F
L
>F
R
, and vice-versa.
51
0 0.2 0.4 0.6 0.8 1
False positive rate
0
0.2
0.4
0.6
0.8
1
True positive rate
Figure 4.4: Individual ROC curves of all participants and overall ROC curve across all
participants. The colored lines represent each individual participant's ROC curve derived
from the normalized outputs of all leave-one-out cross-validated classiers. The bold black
curve is the overall ROC across all participants, and is determined using all participants'
normalized classier outputs. The area under the curve (AUC) of the overall ROC curve
is 0.72.
4.4 Discussion
We observed a goal-equivalent variability in individual forces during repetitions of the
nger force task in all of our participants. This result is consistent with previous studies
with similar multi-nger tasks (Latash et al. 1998, Latash et al. 2001): when two or
more ngers of a hand applied forces that matched to the sum of forces of all employed
52
Participant AUC
Number of rep-
etitions in each
class: F
L
>F
R
/
F
R
>F
L
1 0.6278 29=29
2 0.5374 21=21
3 0.8519 33=33
4 0.6367 30=30
5 0.6181 31=31
6 0.6475 27=27
7 0.6748 32=32
8 0.7586 31=31
9 0.4100 19=19
10 0.8117 33=33
11 0.6576 31=31
12 0.7931 29=29
13 0.7234 26=26
14 0.7622 30=30
15 0.6469 32=32
16 0.5956 19=19
17 0.7211 30=30
18 0.5843 26=26
Table 4.1: AUC values and number of repetitions in each class for each participant.
53
ngers, a variability in individual forces is seen over repetitions while keeping the sum of
forces constant. Other types of movement repetitions, e.g., a sit-to-stand task (Scholz &
Sch oner 1999) or walk to run transitions (Seay, Haddad, Van Emmerik & Hamill 2006),
have also been associated with goal-equivalent variationiability such as at the level of
joint congurations. We additionally found that pre-movement EEG signals could make
signicant predictions of goal-equivalent force variation. We observed spatially distinct
cortical patterns associated with the goal-equivalent variability: a relatively higher EEG
variance was localized on the left when higher left nger force was applied, and vice-
versa. We conclude that goal-equivalent variability is not entirely generated by subcortical
circuitry, and that the cortex may be involved in coordinating muscle activity in a goal-
equivalent way.
While our study is the rst to examine the predictive power of cortical signals associ-
ated with automatic variations in goal-equivalent motor output in a bimanual task, our
results are consistent with previous work predicting hand choice during voluntary uni-
manual tasks. The lateralization in brain activity that we determined to be predictive of
goal-equivalent motor output variation is consistent with previous studies of EEG-based
prediction of unilateral hand movement intention (M uller-Gerking et al. 1999, Ramoser
et al. 2000). These studies found that during the imagination and preparation of left
hand movements, CSP patterns indicate a high EEG variance localized on the left hand
motor area, with low variance on the right (M uller-Gerking et al. 1999, Pfurtscheller
et al. 1997, Ramoser et al. 2000). This pattern switches to a higher variance on the
right hand motor area during repetitions involving the right hand. The ipsilateral con-
centration of EEG variance is interpreted as arising from event-related desynchronization
54
of EEG signals (i.e., decrease in amplitude) on the contralateral sensorimotor area, of-
ten accompanied by event-related synchronization on the ipsilateral side (Pfurtscheller &
Aranibar 1979, Pfurtscheller & Neuper 1997, Pfurtscheller et al. 1997). Our CSP patterns
during the bilateral nger task also reveal a relatively higher EEG variance on the left
hemisphere when relatively higher left force than right was applied, and higher variance
on the right when relatively higher right force than left was applied. We interpret this
nding as evidence that task repetitions with relatively higher left nger force were as-
sociated with more activity in right motor cortex, and task repetitions with relatively
higher right nger force were associated with more activity in left motor cortex.
Previous research has revealed brain-level circuitry that could plausibly implement
goal-equivalent variation in motor output as either feedforward planning or a sensory
feedback strategy. A purely feedforward planning strategy for goal-equivalent motor
output variation for this task would suggest that nonlinear dynamics in intracortical
circuits during task preparation may converge on dierent solutions to the task before
movement onset (Churchland, Byron, Ryu, Santhanam & Shenoy 2006). Variation in
motor output generated during feedforward planning may arise from fundamental
uc-
tuations in neural ring rate at the level of the brain that then propagate through net-
works, as has been observed for cerebellar neurons (Joshua & Lisberger 2014, Medina
& Lisberger 2007), visual neurons (Softky & Koch 1993, Stein, Gossen & Jones 2005),
and motor neurons (Churchland et al. 2006). Alternatively, goal-equivalent motor output
variation could be generated by cortical circuits that are integrating sensory feedback
about downstream factors, such as muscle fatigue (Kouzaki & Shinohara 2006) or tissue
discomfort (Singh, Varadhan, Zatsiorsky & Latash 2010), and adjusting motor output to
55
trade-o among muscles to minimize the eects of these factors. Future studies can begin
to disambiguate these alternatives by manipulating muscle fatigue or skin nociception
(e.g. exerting nger force on a rough surface), and determining if these manipulations
aect cortical signals associated with predicting goal-equivalent variation.
Previous research has also revealed that spinal-level circuitry could plausibly gener-
ate goal-equivalent variability, and our study does not necessarily reject the possibility
that in addition to cortical signals, some of the variability may be encoded by the spinal
circuitry. Spinal re
ex circuits are unlikely to generate coordinated activity during bi-
manual goal-equivalent tasks where one hand applies higher force and the other applies
lower force (Mutha & Sainburg 2009). However, this type of opposing bilateral cou-
pling is present in central pattern generators (CPG) encoded at a spinal level (Grillner &
Wallen 1985). Therefore, goal-equivalent variation could plausibly be generated entirely
at the spinal level by the timing of a constant descending command serving to bias CPG-
like oscillatory circuitry into a right-dominant or left-dominant stable state. It is also
possible that cortical and subcortical-level processes both contribute to goal-equivalent
variation. For example, while pre-movement EEG signals predicted which hand applied
relatively higher force than the other, the origin of the actual amount of variability within
a class was not addressed in this study. While cortical signals may decide on which hand
applies relatively higher force, the variability within a class could also be produced at a
subcortical level. Additionally, while we excluded force data where the two forces were
almost equal for the purpose of the classication technique, the small variations in the
middle portion of the ellipse in Figure 4.2c are still important variations, and may pos-
sibly be generated by other cortical and/or spinal level contributions. A limitation of
56
our current approach is that we do not quantify cortical and subcortical contributions to
goal-equivalent variability relative to each other.
Pre-movement EEG was successful in classifying goal-equivalent variability using al-
pha and beta-band modulations over the sensorimotor area, but our study is limited in
that other frequency modulations of EEG signals, as well as other cortical and subcorti-
cal regions may additionally contain important predictors of goal-equivalent variability.
Pre-movement alpha and beta band desynchronizations of EEG have been very useful in
predicting left/right hand movement intention and type of intended action (Heinrichs-
Graham & Wilson 2016, Pfurtscheller & Aranibar 1979, Turella, Tucciarelli, Ooster-
hof, Weisz, Rumiati & Lingnau 2016), as well as distinguishing between hand and foot
movement intention (M uller-Gerking et al. 1999). Lower frequency modulations of pre-
movement EEG, such as the Bereitschaftspotential or readiness potentiala slow negative
shift in the DC potential on the contralateral M1may also be used in predicting laterality
in self-paced movement intentions (Blankertz et al. 2006, Shibasaki & Hallett 2006). In
addition to sensorimotor modulations of the EEG signal, other cortical and subcortical
brain areas could additionally reveal interesting features associated with variability. Fu-
ture studies can address these limitations by analyzing other frequency modulations in
EEG, and by using fMRI to assess whole-brain predictors of goal-equivalent variability.
Many studies suggest that goal-equivalent variability is important for human health
and motor performance. Goal-equivalent variability has been shown to be benecial in
avoiding overuse injuries (Hamill, Heiderscheit & Pollard 2005, Lipsitz 2002) and minimiz-
ing fatigue (Singh et al. 2010). A decrease in this variability leads to decreased functional
57
ability and overuse injuries (Hamill et al. 2005, Lipsitz 2002). Low goal-equivalent vari-
ability has been found in patients with knee pain (Hamill et al. 2005), runners with low
back pain (Seay et al. 2011) and lower extremity injuries (Hamill et al. 1999), and elderly
individuals vulnerable to falls (Lipsitz 2002). Varying motor output in the goal-equivalent
space may also be an important component of exploration during motor learning (Deeny,
Hau
er, Saer & Hateld 2009, Latash & Anson 2006, Shim, Kim, Oh, Kang, Zatsiorsky
& Latash 2005, Todorov & Jordan 2002, Yang & Scholz 2005) and adapting to changes in
task demand (de Freitas & Scholz 2010). While our study used a single target force level
to study goal-equivalent variability, future studies can explore learning and diculty of
task by mapping out cortical contributions to goal-equivalent variation as a function of
multiple levels of target force.
Goal-equivalent variability requires coordinated activity among a set of muscles. We
have recently shown that automatic coordination of muscles, even across multiple body
segments, likely involves activation and functional connectivity of specic motor corti-
cal regions (Asavasopon, Rana, Kirages, Yani, Fisher, Hwang, Lohman, Berk & Kutch
2014, Rana, Yani, Asavasopon, Fisher & Kutch 2015). We have also shown that the
motor cortical regions that may coordinate this muscle activity become dysfunctional
in chronic pain (Kilpatrick, Kutch, Tillisch, Nalibo, Labus, Jiang, Farmer, Apkarian,
Mackey & Martucci 2014, Kutch, Yani, Asavasopon, Kirages, Rana, Cosand, Labus, Kil-
patrick, Ashe-McNalley & Farmer 2015, Woodworth, Mayer, Leu, Ashe-McNalley, Nal-
ibo, Labus, Tillisch, Kutch, Farmer & Apkarian 2015). Given that chronic pain and
fatigue are common comorbidities (Clauw & Chrousos 1998), and goal-equivalent motor
variation could act as a fatigue mitigation strategy, our current results suggest future
58
research on the link between motor cortical dysfunction and the fatigue associated with
chronic pain.
59
Chapter 5
Voluntary nger force variation in EEG activity
5.1 Introduction
Our EEG study of a bimanual nger force task indicated a clear relationship between
cortical signals and automatic nger force variations. In chapter 4, we asked healthy
participants to perform a bimanual nger force production task where they simultaneously
applied left and right nger forces to repeatedly match the sum of the two forces to a
single target level. Individual nger forces were never displayed; and as a result, right
and left forces
uctuated across repetitions while the target force (sum of forces) was held
constant. Our results showed that EEG signals from the preparation period associated
with the group of repetitions having higher right nger forces than left were clearly distinct
from those associated with the group of repetitions with higher left forces than right.
Moreover, the location of the distinct signals from each group indicated a motor cortical
involvement in producing goal-equivalent variability, and a left-right neural agreement
with right and left output forces. We concluded that cortical signals were associated with
automatic, goal-equivalent nger force variations.
60
In this chapter, we hypothesize that a similar target matching task where individual
forces are voluntarily controlled would yield comparative results to our previous EEG
results of automatic nger force variation. We asked our participants to repeatedly match
their sum of index nger forces to a 6N target, but this time, individual forces were also
required to match individual targets from two possible combinations of forces, randomly
assigned across repetitions. In half of the repetitions, we asked participants to match
their sum of forces to a target with much higher right force than left, and in the other
half, with much higher left force than right. Individual left and right forces were displayed
on the feedback screen. We applied the same machine learning paradigm to analyze the
EEG signals from the 500 ms before movement onset using these two groups of force
matching repetitions (left>right or right>left), as an attempt to nd patterns in cortical
signals associated with voluntary variability. We hypothesized that these patterns would
be similar to those we found in the automatic nger force variation; that is, higher left
motor cortical activity associated with right>left repetitions, and higher right motor
cortical activity associated with left>right repetitions.
5.2 Methods
Participant Population
We recruited 18 healthy right-handed participants (9 males and 9 females) with a mean
age 27.4 4.5 years (range, 23 41). The studies we describe here were performed at the
University of Southern California and approved by the University of Southern California
Institutional Review Board. All participants provided informed consent.
61
Overview of experiments
We used a two-dimensional bilateral nger force task where participants applied simul-
taneous left and right index nger forces to match one of two possible targets shown on
the feedback screen: one where the left nger force was clearly higher than right, and
another where the right nger force was clearly higher than left. In both cases, the sum
of forces was the same, and individual left and right force feedback was displayed on the
screen. We used electroencephalography (EEG) to measure cortical signals during repe-
titions of the task, and employed a machine learning technique to nd out if EEG signals
just prior to the start of each repetition were correlated with force output. Our analysis
was based on classifying force output (two classes, left>right or right>left) using EEG
signals only, and comparing the EEG-based classication with the force-based classes.
Finger force task
Participants performed a bimanual nger force task that involved abducting both index
ngers simultaneously against load cells (Measurement Specialties, Hampton, VA) (Figure
4.1a). We asked participants to repeatedly use both ngers to match a target total nger
force (left force [F
L
] + right force [F
R
]) of 6 Newtons (N). There were two possible
targets where the exact individual force levels were to be matched: one target required
a combination of 90% left force and 10% right force summing up to 6N, and another
that required a combination of 90% right force and 10% left force. Participants received
continuous visual feedback of their current individual nger forces as a dot on the feedback
screen which moved to the right as participants applied right nger force, and moved
upwards as they applied left nger force. The target to be matched was also shown as
62
a circle. Each repetition of the task included the following segments chosen based on
previous literature (Blankertz et al. 2006; Mller-Gerking et al. 1999): 1) 2 seconds of rest
where nothing was shown on the feedback screen, 2) 2 seconds of a preparation period in
which participants were cued to prepare for the task, and 3) 5 seconds of force hold at
one of the two possible target locations, chosen randomly. Participants performed a total
of 60 repetitions.
We show an example of one participant's force repetitions in Figure 5.1b. During
each repetition of the task, total force was kept at 6N (top trace), while individual forces
were randomly assigned to either (F
R
= 0:6N;F
L
= 5:4N) or (F
R
= 5:4N;F
L
= 0:6N).
We computed separate averages for F
R
and F
L
across time starting one second after the
target was displayed (to allow for reaction time and force stabilization) and ending half a
second before the target disappeared (to exclude the return to baseline force). A total of
60 combinations of F
R
versus F
L
along with a 95% covariance ellipse for one participant
are shown in Figure 5.1d. The major axis of the covariance ellipse represents the goal-
equivalent direction along which all combinations do not aect the task of producing
6N of total nger force. The orthogonal direction is the non-goal-equivalent direction.
We identied the two groups of force hold repetitions (F
R
= 0:6N;F
L
= 5:4N) and
(F
R
= 5:4N;F
L
= 0:6N) as two classes, F
L
>F
R
(shown in blue) and F
R
>F
L
(shown
in red), respectively. We also show the rectied EMG signals associated with the task
repetitions in Figure 5.1c as evidence that relative EMG activity was consistent with
relative forces.
63
0
6
F
R
+ F
L
(N)
0 50 100 150
Time (sec)
0
6
F
R
, F
L
(N)
Task: 2D task with 2
possible targets
randomly selected from
...
...
a b
F
R
= 0.6N, F
L
= 5.4N
and
F
R
= 5.4N, F
L
= 0.6N
(F
R
+ F
L
= 6N)
F
L
F
R
F
L
> F
R
F
R
> F
L
-2 2 6
-2
2
6
c
F
R
(N)
F
L
(N)
0 50 100 150
Time (sec)
0
2
rectified EMG
...
Figure 5.1: Experiment setup and force output. a. The participant places their hands

at on a table and repeatedly exerts left (F
L
) and right (F
R
) index nger forces simul-
taneously, while keeping the sum of the two constant at a target force of 6 Newtons (N).
Only the sum of forces is shown to the participant as feedback on a screen, along with
the target force level. Applying a right force moved the circular target to the right, while
applying a left force moved it up. Two possible targets are displayed randomly selected
from (F
R
= 0:6N;F
L
= 5:4N) or (F
R
= 5:4N;F
L
= 0:6N); b. example of a representa-
tive participant's repeated force holds. The top trace shows the participant maintained
the sum of the two forces at 6N, while the middle traces show the left (in blue) and
right (in red) forces for the same repetitions indicating the two possible targets. The
interval shaded in gray (1 second after the start and 0.5 second before the end of the
hold period) is used to determine the time average of individual forces in each repetition.
The bottom plot shows the rectied EMG of right (red) and left (blue) FDI muscles; c.
average right versus left forces during all task repetitions for the same participant with a
95% covariance ellipse t.
5.3 Results
We found higher variance along the goal-equivalent direction than the non-goal-equivalent
direction in all 18 participants. The design of our task was such that a goal-equivalent
variation was expected for all participants. In Figure 5.2, we show the 95% covariance
ellipses of all 18 participants, along with the goal-equivalent line. As expected, igure 4.2b
shows that the ratios of goal-equivalent variance to the sum of both non-goal-equivalent
64
and goal-equivalent variances in all participants for the force level of 6N were signicantly
greater than 0.9 (p< 0:0001).
Our LME model indicated that force-based repetition class had a signicant but
weak eect on EEG-based classier output (p = 0:035). Figure 5.3a illustrates the EEG-
based classier output of each participant averaged across repetitions from each class,
compared with actual force-based classes. The x-axis represents the label of force-based
classes, while the y-axis is the LDA classier output f, a linear combination of the
natural log of channel variances of spatially-ltered EEG data. The EEG-based classier
was trained to produce lower values for training repetitions for whichF
L
>F
R
compared
to training repetitions for which F
R
> F
L
. Our results showed that test repetitions for
which F
L
>F
R
were not always associated with lower EEG-based classier output than
for F
R
>F
L
. For each participant, a positive slope from F
L
>F
R
to F
R
>F
L
indicates
a successful prediction of the two classes. Some participants had a positive slope, while
others a negative slope, resulting in a weak association between force-based class and
EEG-based classier output.
To visualize patterns of brain activations during each class of repetitions, we mapped
average spatial patterns, CSPs, from each class onto the 27 electrode locations on the
scalp (Figure 5.3b). The top and bottom rows represent the most important and second
most important patterns associated with the most important and second most important
spatial lters, respectively, that were used to predict EEG-based classes. Left and right
columns represent the patterns for the two classes, F
L
>F
R
and F
R
>F
L
, respectively.
The CSP patterns illustrate the amount of EEG variance at each electrode, which played
an important role in the discrimination between classes. The averaged patterns did not
65
indicate clear laterality on sensorimotor areas, i.e., a relatively higher EEG variance on
the left hemisphere during repetitions belonging to class F
L
>F
R
, and vice-versa.
Variance Ratio
****
F
R
(N)
F
L
(N)
b a
0.9
0.92
0.94
0.96
0.98
1
-2
2
6
-2 2 6
Figure 5.2: Covariance ellipses and force variance ratio of 18 participants. a. 95%
covariance ellipses of all participants. The covariance ellipse of each participant was
constructed using the average left and right forces from all repetitions of the task (see
Figure 5.1c). b. boxplot of the ratio of goal-equivalent force variance to the sum of goal-
equivalent and non-goal-equivalent variance. The ratio was signicantly greater than 0.9
(p< 0:0001, ****).
66
Force-based class
Predicts F
R
> F
L
Predicts F
L
> F
R
F
L
> F
R
F
R
> F
L
F
L
> F
R
F
R
> F
L
EEG-based classifier output
b
Most important CSP
Second most important CSP
0 0.1
a
Figure 5.3: EEG-based prediction of force output and averaged common spatial patterns
(CSP). a. EEG-predicted classes versus force-based classes of all 18 participants. For each
participant, the classier outputs of repetitions from each class (F
L
>F
R
and F
R
>F
L
)
were averaged and joined by a line. The participants do not all show a clear distinction
between the two classes, as indicated by a mix of positive and negative slopes of the line.
Linear mixed-eect modeling demonstrated a signicant eect of force-based repetition
class on EEG-based classier output (p = 0:035); b. the most important (top row) and
second most important (bottom row) common spatial patterns (CSP) used for the EEG-
based classication of repetitions where F
L
> F
R
(left) and repetitions where F
R
> F
L
(right), averaged across 18 participants.
67
5.4 Discussion
We designed the current task to mimic goal-equivalent variability in individual forces
during repetitions of the nger force task, but this time, in a more controlled manner.
We hypothesized that this voluntary variation in individual nger forces would yield
comparative cortical signals to the automatic nger force variation observed in Chapter 4.
While all of our participants had higher force variation along the goal-equivalent direction
than the non-goal-equivalent direction (as expected), our pre-movement EEG signals
did not necessarily make strong predictions of voluntary goal-equivalent force variation.
The EEG patterns associated with the predictions did not indicate clear motor cortical
laterality associated with F
L
> F
R
and F
R
> F
L
. Although our LME model yielded
signicant association between EEG prediction and repetition-based force output, the
patterns associated with the prediction did not show clear distinction between left and
right patterns.
Previous research has extensively investigated EEG predictions of voluntary uniman-
ual left versus right nger force production. Many Brain-Computer Interface (BCI) stud-
ies aim to automatically predict movement intention from EEG signals alone (M uller-
Gerking et al. 1999, Ramoser et al. 2000, Blankertz et al. 2006, Blankertz et al. 2008).
These studies found that during the imagination and preparation of left hand movements,
CSP patterns indicate a high EEG variance localized on the left hand motor area, with
low variance on the right (M uller-Gerking et al. 1999, Pfurtscheller et al. 1997, Ramoser
et al. 2000). This pattern switches to a higher variance on the right hand motor area dur-
ing repetitions involving the right hand. The ipsilateral concentration of EEG variance is
68
interpreted as arising from event-related desynchronization of EEG signals (i.e., decrease
in amplitude) on the contralateral sensorimotor area, often accompanied by event-related
synchronization on the ipsilateral side (Pfurtscheller & Aranibar 1979, Pfurtscheller &
Neuper 1997, Pfurtscheller et al. 1997). Our CSP patterns during the bilateral voluntary
nger force task, where F
L
>F
R
or F
R
>F
L
, did not indicate this type of lateralization
in EEG activity (Figure 5.3b).
We believe that a number of issues were responsible for the non-agreement between
the two. First, the bimanual nature of our voluntary task may have had dierent and
more complex neural signatures in the EEG signals observed, compared to the unimanual
and more distinct task often employed in BCI studies. Second, in terms of comparing
our results to our bimanual automatic goal-equivalent variation results of Chapter 4, we
believe that the voluntary bimanual task was much more complex and unnatural than the
goal-directed task of automatic nger force variation. That is, individually controlling
each index nger is not a common task, and learning may need to be involved. And
nally, the 2-dimensional display of the visual feedback may have been complicated to
follow, as applying right force moved the target to the right but applying left force moved
the target up. This could have confused our participants and resulted in the unexpected
results we found.
69
Chapter 6
Goal-equivalent nger force variation in fMRI brain activity
6.1 Introduction
Our EEG study of a coordinated bimanual nger force task of Chapter 4 demonstrated
distinct cortical maps associated with automatic variable force combinations. While EEG
signals contain high temporal resolution useful to nd oscillatory features in the machine
learning technique we used, their spatial resolution is insucient to make inferences about
the exact location(s) of cortical activity associated with goal-equivalent variability. In this
chapter, we verify our EEG results of goal-equivalent variability, and use functional Mag-
netic Resonant Imaging (fMRI) to scan whole brain activity while participants repeatedly
perform the same bimanual task. The aim of this study is to determine if Blood Oxy-
gen Level Dependent (BOLD) signals may detect motor cortical activity correlated with
goal-equivalent variability, thus supporting our previous results using a dierent modality.
70
6.2 Methods
Participant population
We recruited 9 healthy participants (5 males and 4 females, 7 right-handed, 1 left-
handed and 1 ambidexterous) with a mean age 27:3 3:4 (range, 25 { 36). The studies
we describe here were performed at the University of Southern California and approved
by the University of Southern California Institutional Review Board. All participants
provided informed consent.
Overview of experiments
We used a task similar to our EEG experiment (Babikian, Kanso & Kutch 2017)
where participants applied simultaneous left and right index nger forces to repeatedly
match the sum of their forces to a target feedback force. We used functional Magnetic
Resonance Imaging (fMRI) to simultaneously measure brain activation associated with
goal-equivalent variability in output forces. Our analysis was based on nding a relation-
ship between the ratio of left index nger M1 activity (right M1) to right index nger M1
activity (left M1), and the ratio of left index nger force to right index nger force. A
positive correlation would imply that the automatic variation in motor output is encoded
in the primary motor cortex.
Finger force task
Participants performed a bimanual nger force task that involved abducting both index
ngers simultaneously against pressure pads (Biopac Systems, Inc., California, USA). We
used pressure pads to ensure compatibility with MRI. Participants lied in supine position
71
in the scanner and held to a custom-made handle bar where they could extend their
index ngers to abduct against the pressure pads. We chose a single target force level of
2 Newtons (N) that was easy to achieve for all. Participants received continuous visual
feedback of their current total nger force as well as the xed target, without individual
force feedback. Each repetition of the task consisted of 5 seconds of rest where a xation
crosshair was shown, and 10 seconds of target hold with active feedback of the sum of
forces. Participants performed a total of 26 repetitions each.
An example of one participant's force repetitions is shown on Figure 6.1b. During
each repetition of the task, the sum of forces was kept at 2.5N (top trace), while individual
forces
uctuated.The participant performed some repetitions with higher right nger force
than left (F
R
>F
L
), and other repetitions with higher left nger force than right (F
L
>
F
R
). To summarize the combination of forces used by the participant on each repetition of
the task, we computed separate averages for F
R
and F
L
across time starting one second
after the target was displayed (to allow for reaction time and force stabilization) and
ending half a second before the target disappeared (to exclude the return to baseline
force). A total of 26 combinations of F
R
versus F
L
along with a 95% covariance ellipse
for one participant are shown in Figure 6.1c. The major axis of the covariance ellipse
represents the goal-equivalent direction along which all combinations do not aect the
task of producing 2.5N of total nger force. The orthogonal direction is the non-goal-
equivalent direction.
fMRI acquisition
We measured brain activation associated with goal-equivalent variability using fMRI.
72
0
2.5
0 100 200 300 400
Time (s)
0
1.25
0 2.5
0
2.5
F
R
+ F
L
(N) F
R
, F
L
(V)
F
L
> F
R

F
R
> F
L
F
R
(N)
F
L
(N)
Task:
a b
F
R
+ F
L
= 2.5N
F
L
F
R
c
Figure 6.1: Experiment setup and force output. a. The participant lies down in the MRI
scanner and repeatedly exerts left (F
L
) and right (F
R
) index nger forces simultaneously,
while keeping the sum of the two constant at a target force of 2.5 Newtons (N). Only
the sum of forces is shown to the participant as feedback on a screen, along with the
target force level; b. example of a representative participant's repeated force holds. The
top trace shows the participant maintained the sum of the two forces at 2.5N, while the
bottom traces show the left (in blue) and right (in red) forces for the same repetitions.
The interval shaded in gray (1 second after the start and 0.5 second before the end of the
hold period) is used to determine the time average of individual forces in each repetition;
c. average right versus left forces during all task repetitions for the same participant with
a 95% covariance ellipse t.
We used a 3-Tesla (GE Signa Excite) with an eight-channel head coil. We positioned
participants supine while viewing a xation crosshair, and placed foam pads to limit
head motion. We collected T2-weighted echo planar image volumes with blood oxygen
level- dependent (BOLD) contrast (echo time, 34.5 ms;
ip angle, 90; eld of view, 220
mm; pixel size, 3.43 mm) continually every 2.5 s during the imaging run. Each volume
consisted of 37 axial slices (3 mm slice thick- ness, 0.5 mm interslice gaps) that covered the
brain from vertex to cerebellum. We additionally acquired a T1-weighted high-resolution
anatomical image from each participant.
We preprocessed each participants fMRI data using the Functional Magnetic Reso-
nance Imaging of the Brain Expert Analysis Tool (FMRIB FEAT, http://fsl.fmrib.ox.ac.uk/fsl/),
73
which included skull extraction using the brain extraction tool (BET) in FSL (FMRIB
Software Library), slice timing correction, motion correction, and spatial smoothing us-
ing a Gaussian kernel with full-width half-maximum of 5 mm and nonlinear high-pass
temporal ltering (100 s). We extracted appropriate signals given the hemodynamics
response time.
Correlation of force output with primary motor cortex activity
After extracting the BOLD signals, we dened regions of interest (ROIs) consisting of
the left and right Primary Motor Cortices (M1s) corresponding to the voluntary activation
of right and left index ngers (from http://neurosynth.org). The ROI voxel locations were
in normalized units (MNI); and we used FSL along with each participant's structural
brain image to map the voxel locations contained in each ROI from standard MNI space
to participant-specic coordinates.
The analysis of BOLD signals across left and right M1s was carried out at the partici-
pant level. For each participant, the BOLD signals from the voxels corresponding to each
M1 were extracted. Then, for each M1, the BOLD signals were averaged across its voxels
to obtain one time series for M1-left and another time series for M1-right. The BOLD
signals during target matching repetitions were extracted from the timing-corrected sig-
nals.
To understand the relationship between left and right M1 BOLD signals and right and
left forces, for each participant, we computed the following ratios each repetition:
F
R
F
L
and
BOLD-M1
L
BOLD-M1
R
. We then examined the relationship between these values across all repetitions
by tting a line to
F
R
F
L
versus
BOLD-M1
L
BOLD-M1
R
: a positive slope would indicate that as the right
74
nger applied relatively higher force than left, the left M1 had relatively higher neural
activity with respect to right M1.
For group analysis, we used a linear mixed-eects (LME) model to quantify the as-
sociation between
BOLD-M1
L
BOLD-M1
R
and
F
R
F
L
. The LME model of the BOLD ratio
BOLD-M1
L
BOLD-M1
R
contained a xed-eects term for force ratio
F
R
F
L
:
BOLD-M1
L
BOLD-M1
R

F
R
F
L
+: (6.1)
Prior to tting the LME model, the distribution of
BOLD-M1
L
BOLD-M1
R
was converted to a
standard normal distribution (mean 0, std = 1) for each participant to ensure the data
was comparable across all participants. We examined the p-value for the xed-eect
of force ratio on BOLD ratio, and considered p < 0:05 as evidence for a statistically
signicant relationship between cortical activity and goal-equivalent variation in nger
forces.
6.3 Results
We found that all 9 participants achieved the task with a signicantly higher variance of
repeated forces along the goal-equivalent direction than the non-goal-equivalent direction.
In Figure 6.2a, we show the 95% covariance ellipses of all 9 participants, along with the
goal-equivalent line. The relative variance along the goal-equivalent direction is greater
than in the orthogonal direction in all participants, as illustrated in Figure 4.2b. The ra-
tios of goal-equivalent variance to the sum of both non-goal-equivalent and goal-equivalent
variances in all participants were signicantly greater than 0.5 (p< 0:0001).
75
0 1.25 2.5
0
1.25
2.5
0.5
0.6
0.7
0.8
0.9
1
Variance Ratio
b a
****
F
R
(N)
F
L
(N)
Figure 6.2: Covariance ellipses and force variance ratio of 9 participants. a. 95% covari-
ance ellipses of all participants. The covariance ellipse of each participant was constructed
using the average left and right forces from all repetitions of the task (see Figure 6.1c).
b. boxplot of the ratio of goal-equivalent force variance to the sum of goal-equivalent and
non-goal-equivalent variance. The ratio was signicantly greater than 0.5 (p < 0:0001,
****).
The relationship between force ratio and BOLD ratio for the 26 repetitions of one
participant is shown in Figure 6.3. The gure indicates that as the ratio
F
R
F
L
increased,
the BOLD ratio
BOLD-M1
L
BOLD-M1
R
also increased, as shown by the positive slope of the tted
line.
76
-0.04 -0.02 0.02
Increasing F
R
2
BOLD Ratio (M1 Left)/(M1 Right)
×10
-3
0
-4
-2
0
Figure 6.3: Relationship between force ratio and BOLD ratio for the 26 repetitions of
one participant. The x-axis represents a relative increase in right force, while the y-axis
is a normalized BOLD ratio of left M1 to right M1. The tted line has a positive slope
indicating a relatively higher right force was associated with a relatively higher left M1
activity.
Our LME model for group analysis indicated that force ratio had a signicant eect
on BOLD ratio (p< 0:05). Figure 6.4 illustrates the normalized BOLD ratios
BOLD-M1
L
BOLD-M1
R
of all participants compared with their force ratios
F
R
F
L
across all repetitions. The tted
line had a positive slope of 0.134 with p = 0:043. Our results showed that repetitions
for which F
R
> F
L
were generally associated with BOLD-M1
L
> BOLD-M1
R
across
participants.
77
-2 -1
-2
0
2
Increasing F
R
BOLD Ratio (M1 Left)/(M1 Right)
0 1 2 3
Figure 6.4: Relationship between force ratio and BOLD ratio for the all repetitions of all
participants. The x-axis represents a relative increase in right force, while the y-axis is
a normalized BOLD ratio of left M1 to right M1. The tted line has a positive slope of
0.134 (p< 0:05) indicating a relatively higher right force was associated with a relatively
higher left M1 activity during automatic goal-equivalent force variation.
6.4 Discussion
We observed a goal-equivalent variability in individual forces during repetitions of the
nger force task in all of our participants. This result is consistent with previous studies
with similar multi-nger tasks (Latash et al. 1998, Latash et al. 2001, Scholz & Sch oner
78
1999, Seay et al. 2006), and also with our goal-equivalent force output results of Chapter 4:
when two or more ngers of a hand applied forces that matched to the sum of forces
of all employed ngers, a variability in individual forces is seen over repetitions while
keeping the sum of forces constant. We additionally found that goal-equivalent nger
force variation was associated with cortical activity. In particular, we observed relatively
higher neural activity in the left M1 than right when relatively higher right force than
left was applied, and vice-versa. While fMRI studies have shown an increase in BOLD
signals associated with increasing voluntary contralateral FDI activation, we observed
an association between relative right and left BOLD signals and automatic changes in
relative left and right forces. This would conclude an agreement with our experimental
results that automatic force variability is at least in part produced in the primary motor
cortex, shown using two modalities: EEG which measures cortical activity with high
temporal resolution, and fMRI which provides high spatial resolution.
Although there was signicant association between left and right BOLD signals of
M1 and goal-equivalent nger force variation, this association was not as strong as our
EEG results suggested (p-value just below 0.05). We believe that our study had some
limitations that could not be addressed at the time of data collection because of lack of
equipment or resources; but that could be modied in future studies for a better study
design.
Some of the limitations that we believe aected our results follow. First, our resources
and timing of data collection allowed the recruitment of only 9 subjects, which may be
an insucient number for signicant results. Moreover, due to time restrictions during
scanner runs, we were only able to collect 26 repetitions per subject. The hold time of the
79
target matching portion of the task (5 seconds) may also have been too short for a very
low resolution of image samples (one 3D image every 2.5 seconds), which yielded only
one BOLD value per repetition after timing corrections. Second, our pressure pads that
measured ngertip forces were not always accurate, as their gains often
uctuated across
participants. Future studies could solve this issue by using higher quality force sensors
compatible with fMRI. Finally, due to our pressure pads saturating at higher forces,
we were obligated to use a low target (sum of forces) of 2.5 N, resulting in very small
changes in individual forces across repetitions. These changes may have been too small to
be detected in the contralateral BOLD signals. Again, better force sensors would allow
for higher force targets resulting in higher relative changes in individual nger forces,
hopefully to be better detected in BOLD signals.
80
Chapter 7
Goal-equivalent variability explained by interhemispheric
interactions
7.1 Introduction
Goal-equivalent variability is a dominant feature of human movements. In Chapter 4 , our
experimental study of a coordinated bimanual task demonstrated a correlation between
cortical activity and goal-equivalent variability in output forces (Babikian et al. 2017).
Our head maps of cortical activity showed that trial-to-trial variability may therefore be
produced, at least partly, by cortical signals in the primary motor cortices (Figure 4.3).
The hand area of the primary motor cortex (M1) controls hand movements and force
productions. Coordinated bimanual tasks, such as that we used in our experimental study,
rely heavily on callosally mediated communication (Carson 2005, Fling & Seidler 2012),
i.e., neural interactions between the left and right hemispheres. These interhemispheric
interactions may consist of facilitatory (excitatory) and/or inhibitory connections. At
the level of the primary motor cortices, the interaction is primarily inhibitory, i.e., there
is mutual inhibition between the two M1s (Ferbert, Priori, Rothwell, Day, Colebatch &
81
Marsden 1992, Kujirai, Caramia, Rothwell, Day, Thompson, Ferbert, Wroe, Asselman &
Marsden 1993, Grefkes, Eickho, Nowak, Dafotakis & Fink 2008). This interhemispheric
inhibition helps prevent interference of control processes between the two opposing cor-
tices (Fling, Peltier, Bo, Welsh & Seidler 2011).
Interhemispheric interactions have been well-studied using Transcranial Magnetic
Stimulation (TMS) (Ferbert et al. 1992, Kobayashi & Pascual-Leone 2003), a non-invasive
procedure that directly but non-invasively stimulates nerve cells in the brain using mag-
netic elds. In studies exploring interhemispheric interactions at the level of the hand, a
rst stimulus, called the \conditioning stimulus" is sent to one of the M1s, and a second
\test stimulus" is sent to the opposing M1 a few milliseconds later (Figure 7.1). The
electric potential is recorded in the (FDI) muscle contralateral to the motor cortex that
received the test stimulus. This is called the Motor-Evoked Potential (MEP). The MEP
after a pair of conditioning and test stimuli is compared to the MEP from a test stimulus
alone (i.e., without conditioning stimulus). If the MEP with conditioning is lower than
the MEP without conditioning stimulus, then the interhemispheric interaction between
the opposing sites is said to be inhibitory. Interhemispeheric interactions are often quan-
tied as the amplitude of conditioned MEP expressed as a percentage of the MEP from
test stimulus alone. Values less than 100% would imply interhemispheric inhibition, while
values greater than 100% would imply facilitation.
82
MEP
MEP
measured
in right FDI
a b
Figure 7.1: Interhemispheric inhibition measured using transcranial magnetic stimulation
(TMS). a. A conditioning stimulus is sent to the right hand motor area, followed by a
test stimulus to the left motor area a few milliseconds later. The motor evoked potential
(MEP) is measured in the right FDI muscle, contralateral to the test stimulus. b. The
measured MEP is compared to the control case where no conditioning stimulus was
applied. The MEP amplitude is lower with conditioning (bottom trace) than without
(top trace), which indicates interhemispheric inhibition (IHI) between left and right hand
motor cortices. Figures adapted and modied from Ibey et al. (2015) and Kobayashi &
Pascual-Leone (2003).
To understand the nature of interhemispheric inhibition with respect to contralateral
activity, the conditioning stimulus can be changed to simulate varying amounts of activ-
ity in the contralateral M1. For example, Ni et. al. used a conditioning-test paradigm
and measured First Dorsal Interosseous (FDI) muscle activation changes to explore inter-
hemispheric inhibition (IHI) in the hand area of the primary motor cortex (Ni, Gunraj,
83
Nelson, Yeh, Castillo, Hoque & Chen 2009). The inhibition at the left M1 as a func-
tion of right M1 conditioning intensity followed a non-linear function resembling a sig-
moid; it increased with conditioning stimulus intensity following a weak non-linearity (see
Figure 7.2). Other studies have also shown non-linear inhibition at the primary motor
cortices. For instance, Harris-Love et. al. and Morishita et. al. both found a similar non-
linear relationship between conditioning stimulus intensity and amount of inhibition at the
ipsilateral FDI muscle of the hand (Harris-Love, Perez, Chen & Cohen 2007, Morishita,
Uehara & Funase 2012). Similarly, other hand or arm muscles also exhibited non-linear
inhibition as a function of conditioning stimulus intensity; for example, wrist muscles
(Flexor Digitorum Radialis and Extensor Digitorum Radialis Brevis) were investigated
at dierent conditioning stimulus intensities and at dierent inter-stimulus intervals (ISI,
i.e., the time interval between conditioning stimulus and test stimulus) (Ibey et al. 2015).
The inhibition function was again non-linear.
84
a b
Figure 7.2: Interhemispheric inhibition as a function of conditioning intensity in the rst
dorsal interosseous (FDI) muscle (gure from Ni et al. (2009)). The conditioning stimulus
(CS) was delivered to the right M1, and a test stimulus (TS) was delivered to the left M1.
The MEP (motor-evoked potential) was measured in the right FDI muscle at dierent
conditioning stimulus intensities (x-axis). The y-axis indicates the percent change in the
amplitude of muscle activity compared to the control case with test stimulus alone (i.e.
without conditioning). A decrease in %TS indicates inhibition. a. TS delivered 10 ms
after CS, b. TS delivered 50 ms after CS.
Based on our previous experimental results indicating cortical involvement and left/right
laterality in the production of goal-equivalent variability, as well as TMS studies inves-
tigating the nature of IHI in the hand M1, in this chapter, we develop a mathematical
85
model of interhemispheric inhibition as a mechanism to produce goal-equivalent variabil-
ity. Our proposed model is a two-dimensional non-linear system of dierential equations
describing the amounts of activity in either M1 (left and right), where each M1 sends an
inhibitory signal to the opposing side in the form of a sigmoid. A common input to both
M1s is varied as an on/o repetition scheme representing a task similar to our experimen-
tal protocol in Chapter 4. This system belongs to the category of bistable systems, which
comprise of non-linear dierential equations with two distinct stable xed points; and
convergence to either is dictated by initial conditions (Morishita & Yajima 1972, Trotta,
Bullinger & Sepulchre 2012). These types of systems have in common a non-linear feed-
back, and have been employed in modeling various types of biological phenomena, namely,
genetic control (Grith 1968), neural networks (Hopeld 1982, Myre & Woodward 1993),
and species populations dynamics (Lotka 1925, Volterra & Brelot 1931).
We suggest that this mutually-inhibitory model with the addition of independent
noise will allow for a spontaneous switching between the two stable states at dierent
repetitions of the task. Moreover, we nd that the shape of the mutual inhibition plays
a role in making the output of the system reproduce our experimental results of goal-
equivalent variability in output forces. Therefore, we explore the space of parameters
that form the input and inhibition function, with the goal of dening conditions on the
parameters that would result in valid solutions comparable to our experimental results.
86
7.2 Methods
Mutually-inhibiting bistable system
Inspired by a mutually-inhibiting neuronal model (Myre & Woodward 1993), we con-
sidered a system of two variables, x
L
and x
R
, representing the activity in left and right
M1, respectively. There was common input to each x
L
and x
R
, and each M1 included
linear self-inhibition, in addition to non-linear inhibition originating from the opposing
M1. Figure 7.3 shows a sketch of the connections.
M1
left
M1
right
x
L
x
R
F
R
F
L
common
input u
Figure 7.3: A sketch of the connections in the M1-M1 model. The colors represent the
type of interaction: green implies an excitatory input, while red implies inhibition. There
is a common input to both M1s, and mutual inhibition between the two. The activity at
each M1 (x
L
and x
R
) would eventually produce force at the contralateral hand (F
R
and
F
L
, respectively).
87
A linear mutual inhibition does not allow convergence to dierent solutions, because
the system would have one stable xed point, and therefore would always converge to the
same solution. Therefore, we used non-linear functions to represent the mutual inhibition
between x
L
and x
R
. We adapted a neural network model of two mutually-inhibiting
neurons (Myre & Woodward 1993):
dx
L
dt
=x
L
c
1
1 +e
b(x
R
T)
+u
dx
R
dt
=x
R
c
1
1 +e
b(x
L
T)
+u
(7.1)
Here,x
L
andx
R
represent the activity at the left and right M1s, respectively, u is the
common input to both,b is a gain,T is a threshold andc is the strength of the inhibition.
There is linear feedback from each x
i
onto itself, and a non-linear feedback originating
from the opposing x
j
.
The non-linear inhibition
c
1+e
b(x
i
T)
is a sigmoid function of the activity of the op-
posing hemisphere. Figure 7.4 shows this nonlinearity as a function of x with parameters
such that c = 1;b = 5 and T = 0:5. This feedback model was a plausible choice be-
cause interhemispheric inhibition in the hand area has been shown to exhibit a weak
nonlinearity as a function of conditioning stimulus (Ni et al. 2009).
88
0.2 0.4 0.6 0.8 1
x
0.1
0.3
0.5
0.7
0.9
f(x)
0
Figure 7.4: Sigmoid function f(x) =
c
1+e
b(xT)
with c = 1;b = 5 and T = 0:5.
A bifurcation analysis was performed by (Myre & Woodward 1993). We include this
analysis in Figure 7.5. Depending on the parameter set c, b, T and u, the system may
either have one stable xed point, or three xed points: two distinct stable points and a
saddle. In the latter case, initial conditions dictate which xed point the system converges
to.
89
Figure 7.5: Bifurcation analysis. A. One of coordinates of the xed point(s) of sys-
tem (7.1) as a function of gain b. Continuous curves correspond to stable xed points,
dashed curves - to saddle xed points (T = 1, i = 1, c = 1). B. The same as A but
with varying the threshold T (b = 5, i = 1, c = 1). C and D. Bifurcation curves which
separate the regions in parameter space corresponding to the system with one xed point
(denoted by 1) from the regions corresponding to the systems with three xed points
(denoted by 3). Parameter values are T = 0:5, c = 1 (for panel C); T = 0:5, b = 4 (for
panel D). Figure and caption from (Myre & Woodward 1993).
For example, choosing c = 1;b = 5;T = 0:5 and u = 1 results in the phase plane in
Figure 7.6. The two stable points of the system are (0:145; 0:855) and (0:855; 0:145). As
90
long as the initial conditions are such thatx
0L
>x
0R
, the system converges to the stable
point that has x
L
>x
R
, and vice-versa for x
0R
>x
0L
.
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
x
L
x
R
Figure 7.6: Phase-plane of the bistable system in (7.1) with c = 1;b = 5;T = 0:5;u = 1.
The point in the middle is a saddle point [(0:5; 0:5)], whereas the top left and bottom
right points are the two stable states of the system [(0:145; 0:855) and (0:855; 0:145)].
As a summary, the system in (7.1) behaves the following way, depending on parameters
and initial conditions:
1. if it has one stable xed point, the solution always converges to this point.
2. if it is bistable (has two stable points and one saddle point such as in Figure 7.6):
(a) for equal initial conditions, the solution always converges to the saddle point.
91
(b) for unequal initial conditions, the solution always converges to one of the two
stable xed points. If the initial conditions are such that x
0R
>x
0L
, then the
system converges to the stable point corresponding tox
R
>x
L
, and vice-versa.
Therefore, repeated on/o inputs ofu can not result in the spontaneous switch of the
solution from x
R
> x
L
to x
L
> x
R
across repetitions (such as in Figure 4.1B) . As an
example of a repeated input scheme representing repetitions of the target matching task,
Figure 7.7 shows an on-time of 5 time points and an o-time of 4 time points (a). The
system here eventually converged to x
L
=x
R
(b).
80 100
time
0
0.2
0.4
0.6
0.8
1
input
time
0
0.2
0.4
0.6
output
x
L
x
R
0 20 40 60 80 100 0 20 40 60
a b
Figure 7.7: Example simulation of the system in 7.1 with repeated inputs. a. Input of
the system with u = 1 for 5 time points and u = 0 for 4 time points. b. Output of the
system, x
R
and x
L
, with initial conditions x
0L
= 0:1 and x
0R
= 0:11.
Adding noise
A well-studied aspect of neurons and neuronal activity is that neurons exhibit noise
92
even when tested under the same conditions. Variability in the neuronal activity is
present in the primary motor cortex (Lee, Port, Kruse & Georgopoulos 1998), and tends
to be correlated between neurons with high signal correlation (i.e., neurons that carry
redundant messages). We therefore modied the system of equations in (7.1) to add
independent random noise to each of x
L
and x
R
:
dx
L
dt
=x
L
c
1
1 +e
b(x
R
T)
+u +
1
dx
R
dt
=x
R
c
1
1 +e
b(x
L
T)
+u +
2
;
(7.2)
where
1
and
2
are added random noise, each from a gaussian distribution with mean
0 and standard deviation :N (0;
2
).
As an example of the system with added noise, we used the same input repetitions as
in Figure 7.7 applied to the system in (7.2). Figure 7.8 shows the input and simulated
output repetitions. In this case, the response exhibited a spontaneous switching between
x
R
>x
L
and x
L
>x
R
(b).
93
80 100
time
0
0.2
0.4
0.6
0.8
1
input
time
output
x
L
x
R
0 20 40 60 80 100 0 20 40 60
-0.2
0
0.2
0.4
0.6
a b
Figure 7.8: Example simulation of the system in 7.2 with repeated inputs including noise.
a. Input of the system with u = 1 for 5 time points and u = 0 for 4 time points. b.
Output of the system, x
R
and x
L
, with initial conditions x
0L
=x
0R
= 0.
Parameter space exploration
To understand the eects of the parameters on the inhibition function and the output
of the system, we explored dierent combinations of parameters and simulated the system
with repeated inputs (5 on, 4 o input repetition scheme for 200 time points). The free
parameters are c, b, T , u and . We performed 1000 simulations where at each iteration
the parameters were randomly selected from a uniform distribution with the intervals:
c : [0; 5]
b : [0; 20]
T : [1; 2]
u : [0; 7]
94
: [0; 0:2u] (upper limit of 20% of u)
These intervals were adapted from the bifurcation analysis that was performed in (Myre
& Woodward 1993).
For each of the 1000 simulations we performed, we used the following criteria to decide
if the output repetitions represented a valid solution comparable to our experimental
results:
(i) the repetition means are non-negative;
(ii) the major axis of the tted ellipse is approximately parallel to the line joining (0,1),
(1,0);
(iii) the variance ratio > 0:8, as per our experimental results in Figure 4.2.
7.3 Results
An example simulation using the system of equations 7.2 was performed using the pa-
rameter values c = 1, b = 4, T = 0:5, u = 1, and = 0:2 in Figure 7.9. The input
repetitions were "on" for 5 time points, and "o" for 4; and the simulation was run for
1000 s, yielding a total of 110 repetitions. Figure 7.2A shows a few of the simulated
repetitions, indicating the sum ofx
R
andx
L
was approximately constant at 1 across rep-
etitions (top trace in black). The addition of noise allowed for individual activity to vary
from repetition to repetition, resulting in some repetitions being achieved with x
R
>x
L
,
and others with x
L
> x
R
(bottom trace). We extracted (x
L
;x
R
) pairs averaged across
each repetition, and show them in Figure 7.9B, where a 95% covariance ellipse was t to
95
the data. The highest variation was observed across the major axis of the ellipse, where
x
R
+x
L
= 1, while the lowest variation occurred across the minor axis.
...
a b
x
R

x
L

x
R
+ x
L

0.2 0.4 0.6 0.8
0
0.2
0.4
0.6
0.8
1
0 1
...
0
1
x
R
, x
L

x
L
> x
R

x
R
> x
L
0
0.5
0 50 100 150
Time
...
Figure 7.9: Example simulation and covariance ellipse with c = 1, b = 4, T = 0:5, u = 1,
and = 0:2. a. a few of the simulated repetitions. The top trace shows x
R
+x
L
maintained at 1 across repetitions; while the bottom traces show x
L
(in blue) andx
R
(in
red) for the same repetitions. Some repetitions had a relatively higher right activity than
left (x
R
>x
L
), and others vice-versa (x
L
>x
R
). The interval shaded in gray is used to
determine the time average of individual activities in each repetition; b. average right
versus left activities during all 110 repetitions with a 95% covariance ellipse representing
the variance of the data. The major axis represents the goal-equivalent direction, where
the sum of activities is always 1, while the minor axis is the non-goal-equivalent direction
where the sum would deviate from 1. The highest variation is seen in the goal-equivalent
direction, while the orthogonal non-goal-equivalent direction has much lower variance
96
Parameter space exploration
Figure 7.10a shows all 95% covariance ellipses from 1000 simulations that were per-
formed with dierent parameter sets. Each ellipse was t to the repetition means for each
iteration. Figure 7.10b shows the 267 ellipses corresponding to simulations that were de-
termined as being valid solutions based on the criteria listed in the Results section. All
ellipses here have their major axes correspond to the goal-equivalent line where x
R
+x
L
is always constant.
-4
-2
0
2
4
6
-4 -2 8
x
L
x
R
0 2 4 6 -6 -4 -2 8
x
L
-1
1
3
5
7
0 2 4 6
a b
Figure 7.10: Covariance ellipses t to each of 1000 simulations with randomly selected pa-
rameters. a. All ellipses superimposed, b. ellipses from valid parameter sets determined
using the criteria mentioned in the text.
97
Valid inhibition functions
To understand the condition(s) that make the solution valid, we show all inhibition
functions for the valid solutions c
1
1+e
b(x
R
T)
in Figure 7.11 (left), along with all inhi-
bition functions for invalid solutions (right). The lines were color-coded based on the
amplitude ofu at each iteration. Some Inhibition functions for invalid solutions included
linear functions, whereas all inhibition functions for valid solutions were non-linear. More-
over, for valid solutions, the amplitude of u appeared proportional to c, the scaling of
the sigmoid, whereas no clear relationship was seen for invalid solutions. We show this
relationship betweenc and the amplitude ofu in Figure 7.12. In the valid set of solutions
(left), u is proportional to c, whereas no clear relationship is seen in the set of invalid
solutions (right).
98
-4
activity (x)
inhibition
activity (x)
0
1
2
3
4
5
4 8 0 -4 4 8 0
0
1
2
3
4
5
1
2
3
4
5
6
u
a b
Figure 7.11: Inhibition functions of valid (a) and invalid (b) parameter sets. Each curve
represents the inhibition function for one simulation, plotted against x using
c
1+e
b(xT)
with the set of parameters c, b and T corresponding to that simulation. Each curve is
colored with the amplitude of u for the corresponding simulation.
To further investigate the relationship between c andu in valid and invalid solutions,
we performed simulations for xed values ofu. In other words, we ran a 1000 simulations,
7 times, where all parameters except foru were randomly selected at each iteration, while
the amplitude of u was held constant at u = 1; 2; 3; 4; 5; 6; or 7. The results are shown as
box plots in Figure 7.13 for valid (a) and invalid (b) solution sets. In the valid solutions,
there was a clear linear increase ofc with the increase ofu starting atu = 2, while invalid
solutions did not show a relationship between c and u.
99
u
0
1
2
3
4
5
c
c
0 2 4 6
1
2
3
4
5
u
0 2 4 6
0
a b
Figure 7.12: Relationship between inhibition scaling c and input amplitude u for valid
(a) and invalid (b) parameter sets. 1000 simulations were carried out with randomly
selected parameters each time. Valid parameter sets were selected according to the criteria
described in the text. a. For valid parameter sets, the amplitude of u and the scaling of
inhibition c followed a linear relationship. b. Invalid parameter sets did not exhibit any
clear relationship between c and u.
100
u
0
1
2
3
4
5
c
0 1 2 3 4 5 6 7
u
0
1
2
3
4
5
c
0 1 2 3 4 5 6 7
a b
Figure 7.13: Boxplot showing the relationship between inhibition scaling c and input
amplitude u for valid (a) and invalid (b) parameter sets. For each xed u amplitude,
1000 simulations were carried out with the rest of the parameters randomly-selected
from uniform distributions. Valid parameter sets were selected according to the criteria
described in the text. a. The box plots of c with respect to u from valid parameter
sets show a clear linear increase in c with respect to the amplitude of u. b. No clear
relationship was observed between c and u for non-valid parameter sets.
7.4 Discussion
We observed a switching of the solution from x
L
>x
R
to x
R
>x
L
with the addition of
noise to the bistable system in (7.1) (Figures 7.7 and 7.9). Moreover, the means across
repetitions spanned an ellipse such that x
R
+x
L
was constant throughout repetitions
(Figure 7.9). We explored the space of parameters that dene the system in (7.1) and
identied the sets of parameters that agree with our experimental results of Figure 4.2.
101
We found that the inhibition function had an eect on the validity of the solutions. A
non-linear inhibition function was not a sucient condition for the solution to agree with
our experimental results; rather, a linear relationship between c and u was required for
the solution to be valid: the higher the input u, the higher scalingc of the inhibition was
needed.
The non-linear nature of the inhibition function in our mathematical model agrees
with interhemispheric inhibition functions investigated using TMS. These studies have
shown a clear non-linearity of the inhibition with respect to increasing activity of the
contralateral hemisphere (Ferbert et al. 1992, Ni et al. 2009). The inhibition functions
were determined by delivering conditioning and test stimuli to the M1s, each time with
a dierent conditioning stimulus intensity representing contralateral activity. As the
intensity of the conditioning stimulus on one hemisphere was increased, the inhibition at
the contralateral hemisphere where a test stimulus was delivered also increased. In the
hand area and in particular the FDI muscle, the increase in inhibition with respect to
increasing conditioning stimulus intensity was found to be non-linear, often resembling
a sigmoid (Chen, Yung & Li 2003, Ni et al. 2009, Fling & Seidler 2012, Harris-Love
et al. 2007, Morishita et al. 2012).
Our simulations with dierent parameter sets revealed that the non-linearity of the
inhibition function was not enough for the system to replicate our experimental results
of goal-equivalent variability of force output. In addition to having a sigmoidal shape,
the amplitude of the inhibition function for valid solutions scaled with the amplitude of
the input u. To our knowledge, this result has not been shown using TMS. We suggest
further TMS experiments with a conditioning-test paradigm while participants perform
102
the bimanual goal-equivalent nger force task as in Chapter 4, with a few dierent target
levels. During repetitions, a conditioning stimulus would be sent to the FDI hot-spot of
the right M1, with a test stimulus at the left M1 following a few ms after. The MEP
would be measured in the right FDI following the test stimulus and compared to the
MEP without conditioning. This paradigm would be repeated using higher conditioning
stimulus intensities; and the inhibition function would be determined with respect to
conditioning stimulus intensity. The same steps would be taken to determine the inhibi-
tion function at dierent force target levels, representing a scaling of u. The inhibition
functions at dierent target levels would be compared to determine if the scaling of the
inhibition, representing c, increases with u in agreement with our simulation results.
103
Chapter 8
Conclusions
Our results collectively indicate variation in the ne details of human movements across
repetitions of the same task, in addition to brain involvement in producing this variability.
In this dissertation, we explored the plausible solution space of muscle activation patterns
that produce the same movement or force production task. Our mathematical model of
a tendon-driven limb showed that producing smooth limb motions requires variable and
location-specic muscle control precision levels. This result and our experimental study
exploring slow and fast nger movements may explain the origins of the experimentally-
observed non-smoothness in slow nger movements in humans.
To understand the neural origins of goal-equivalent variability, we measured EEG and
fMRI signals while healthy participants repeated a bimanual nger force matching task.
All of our participants used variable relative levels of left and right forces across dierent
repetitions to match to the same target. This variability in output forces agrees with
previous research that studied numerous types of tasks such as running, walking, sit-
to-stand, and multi-nger force production tasks. Our EEG analysis of the preparation
104
period of the force matching task involved using machine learning to nd features in cor-
tical activity correlated with output force variability. Our results revealed motor cortical
involvement in producing automatic left and right nger force variations. Moreover, we
found higher left motor cortex activity when participants applied relatively higher right
nger forces, and higher right activity when they applied relatively higher left nger
forces. This result agreed with previous research involving cortical signals in voluntary
unimanual force production, but our study was the rst to directly measure brain signals
to understand the origins of (automatic) goal-equivalent variation.
We also analyzed Blood Oxygen Level Dependent (BOLD) signals in whole brain fMRI
scans of participants performing the same bimanual nger force task in the scanner. The
high spatial resolution of fMRI signals allowed us to identify the exact locations of left
and right motor cortices associated with the index ngers; and extracted BOLD signals
from the voxels belonging to these locations. We then compared the relative BOLD
amplitudes in the two locations with respect to the relative individual force levels. We
found relatively higher left motor cortex activity than right when participants applied
higher right nger forces than left, and vice-versa, which validated our EEG results of
goal-equivalent variability.
Finally, we built a mathematical model to explain the mechanism driving goal-equivalent
variation in the brain. We based our model on previous experiments that used Tran-
scranial Magnetic Stimulation (TMS) of the brain, and found that interactions between
contralateral motor cortices are primarily inhibitory and follow non-linear behavior with
respect to contralateral activity. The model of mutually-inhibiting left and right motor
cortex activities with the addition of random noise produced goal-equivalent variation
105
in left and right activities across repetitions. We tuned the model parameters for the
output to agree with our experimental results in nger force variation, and found that
the non-linearity of the inhibition function was not enough to match our experimental
results. We deduced that in addition to the non-linearity condition, the scaling (strength)
of the inhibition function must increase with the increase of the input, e.g., the target
level of (sum of) forces. This suggests further experimental studies using TMS aimed at
understanding the changes in the mutual inhibition with respect to dierent target levels
used for the task.
Reduced motor variability has been associated with injuries and chronic pain. Pre-
vious research has shown lower stride-to-stride variations in pelvis to trunk rotations in
runners with current or with a history of low back pain, and higher variability in healthy
runners (Seay et al. 2011). This suggests the possibility that insucient motor variability
may contribute to chronic pain and injury risk; moreover, that goal-equivalent variabil-
ity may be a mechanism aimed at avoiding overuse injuries by varying tissue loading.
Therefore, studying goal-equivalent variability, in addition to contributing to our un-
derstanding of motor control, helps unravel neural signatures and specic brain regions
associated with movement patterns that may predispose an individual to overuse injury
and chronic pain. Future studies can therefore manipulate goal-equivalent variation at
the brain level and the peripheral (output) level to understand how the output variabil-
ity is changed or disrupted. At the brain level, repetitive TMS can be used to reduce
excitability in the motor cortex and determine how the output variability is altered. At
the peripheral level (e.g., nger forces), variability can be manipulated by using a rough
surface at the nger force application location. By understanding where in the nervous
106
system variability is produced and altered, these studies can determine if there is a spe-
cic neural signature that predicts an individual's injury risk, and ultimately help develop
personalized therapy plans that target these brain locations to reduce injury risk and/or
chronic pain.
107
Appendix A
Common Spatial Patterns and Linear Discriminant
Analysis for the EEG prediction of movement intention
This chapter is concerned with explaining the EEG machine learning technique used in
Chapters 4 and 5. We rst explain common spatial pattern ltering for EEG as well as
linear discriminant analysis. We then present two simulated examples with step-by-step
instructions illustrating each of these two techniques.
A.1 Common Spatial Patterns (CSP)
The calibration phase of a classication task consists of extracting features of EEG sig-
nals that are characteristic of each of the two classes (M uller-Gerking et al. 1999). A
spatial ltering technique called Common Spatial Patterns (CSP) has been successful in
improving the prediction of the class of a trial (Fukunaga 1990, Koles 1991). CSP deter-
mines a matrix W that projects the multi-channel EEG data x(t)2<
c
in EEG channel
space into a surrogate space of the same dimension cc, where c is the number of EEG
channels: x
CSP
(t) = W
T
x(t). Each column vector w
j
(j = 1:::c) of the W matrix is a
108
spatial lter, and each column vector a
j
of a matrix A = (W
1
)
T
is a spatial pattern.
CSP maximizes the variance of the projected signal for one class while minimizing it
for the second class, thereby optimizing discriminability with respect to variance. The
CSP method achieves this by projecting the original data from channel space into a new
surrogate space such that the projected data is uncorrelated in both classes. This ensures
maximum dissimilarity along the new directions, making it easier to distinguish between
the two classes based on variance.
Before applying spatial lters, we bandpass lter the EEG data from 7 to 30 Hz, which
includes both alpha and beta rhythms. Then, we group the training set into two classes,
(1) F
L
> F
R
and (2) F
R
> F
L
. We then calculate the pooled covariance matrices C
1
andC
2
corresponding to classes (1) and (2), respectively. The matrixW decorrelates the
covariance matrices of the two classes. Therefore, we apply a simultaneous diagonalization
of the pooled covariance matrices C
1
and C
2
(Blankertz et al. 2008):
W
T
C
1
W =
(1)
W
T
C
2
W =
(2)
;
(A.1)
where
(1)
and
(2)
are diagonal such that
(1)
+
(2)
=I.
This is equivalent to solving the generalized eigenvalue problem:
C
1
w
j
=

(1)
j

(2)
j
C
2
w
j
; (A.2)
where the columns of W , w
j
, are the generalized eigenvectors, or spatial lters.
109
The classier function
The classier outputs a single value that predicts the class of a trial. We follow common
practice in EEG classication literature and use the two most important spatial lters
from each side of the W matrix (a total of four) to project the bandpass ltered EEG
data. The classier function projects the original EEG data using the spatial lters,
nds the log of its power, applies a linear combination and adds a bias. It is dened as
follows (Blankertz et al. 2008):
f =
J
X
j=1

j
log(w
T
j
XX
T
w
j
) +
0
; (A.3)
where w
j
are the lters determined with CSP, J is the total number of lters used
(here, 4), and the linear weights
j
are calculated statistically from the training data
using Fisher's Linear Discriminant Analysis described in section A.2 (Hastie, Tibshirani,
Friedman & Franklin 2005). The idea is that the change in EEG oscillatory activity,
which is equivalent to the band-power or variance of the EEG signal, is captured by the
log-power of the band-pass ltered signal. This way, each projection of the signal captures
dierent spatial localization. For each new repetition, f > 0 is predicted as one class,
and f < 0 is predicted as the other class.
A.2 Linear Discriminant Analysis (LDA)
LDA is a supervised learning technique that nds a linear combination of features that
distinguishes two or more classes. The goal is to maximize class separability based on a
set of features unique to the classes, using training examples. The result is a classication
110
rule, or model, based on the training examples, that best separates the classes. A new
(unseen) example is then classied into one of the groups using a linear classier that is
based on the model found earlier.
The basics of LDA for two-class discrimination will be explained here, and we will
see how it translates into distinguishing EEG class based on variance. The purpose of
LDA is to minimize the total error of classication, which is equivalent to using Bayes'
rule of assigning an object to the class with highest conditional probability. For a set of
features contained in object x, and for g possible groups (or classes), Bayes' rule is to
assign the object to group i for which: P (ijx) > P (jjx);8j6= i. This is equivalent to
assigning group i to x if the conditional probability that this object belongs to group i
given a set of features x is maximum across all groups. According to Bayes Theorem, the
conditional probability can be written as:
P (ijx) =
P (xji)P (i)
P
8j
P (xjj)P (j)
; (A.4)
where prior probabilityP (i) is the probability of groupi regardless of features. In practice,
this is unknown but can be assumed equal for all groups or be calculated based on the
number of examples in each group.
Bayes' rule becomes:
Assign object x to class i if:
P (xji)P (i)>P (xjj)P (j);8j6=i
(A.5)
111
Obtaining P (xji) based on data measurements is impractical because a lot of data
would be required. Therefore, further assumptions are made about the distribution of the
data and this probability is directly calculated. That is, an assumption is made that each
group has n-dimensional multivariate normal distribution, where n is the total number
of features in x. The conditional probability for this distribution has a known equation:
P (xji) =

1
(2)
n
2jC
i
j
1
2

exp

1
2
(x
i
)
T
C
1
i
(x
i
)

, where
i
is the vector mean and
C
i
is the nn covariance matrix of group i.
Bayes' rule then becomes:

1
jC
i
j
1
2
!
exp

1
2
(x
i
)
T
C
1
i
(x
i
)

>

1
jC
j
j
1
2
!
exp

1
2
(x
j
)
T
C
1
j
(x
j
)

;i6=j
(A.6)
Taking the natural log and multiplying both sides by2:
ln(jC
i
j) 2 ln

P (i)

+ (x
i
)
T
C
1
i
(x
i
)< ln(jC
j
j)
2 ln

P (j)

+ (x
j
)
T
C
1
j
(x
j
);i6=j
(A.7)
112
The above classier is referred to as Quadratic Discriminant Analysis (QDA). LDA
makes the further assumption of homoscedasticity, i.e., that both classes have equal co-
variance matrices, C
i
= C
j
= C. In this case, the classier can be further simplied into:
ln(P (i)) +
i
C
1
x
T

1
2

i
C
1

T
i
> ln(P (j)) +
j
C
1
x
T

1
2

j
C
1

T
j
; i6=j;
or
(
j

i
)C
1
x
T

1
2
(
j
C
1

T
j

i
C
1

T
i
) + ln

P (j)
P (i)

< 0; i6=j
(A.8)
The classication is therefore carried out using a projection along two perpendicular
directions, F
1
and F
2
:
F
i
=
i
C
1
x
T

1
2

i
C
1

T
i
+ ln(P (i))
F
j
=
j
C
1
x
T

1
2

j
C
1

T
j
+ ln(P (j))
Classify as class i if f =F
j
F
i
< 0;
Classify as class j if f =F
j
F
i
> 0:
(A.9)
Fisher's LDA does not make the assumption that covariances are equal or that classes
are normally distributed (Fisher 1936). Instead, it denes the separation S between the
two projected distributions as the ratio of between-class variance to within-class variance.
For two classes 1 and 2 having means
1
and
2
, and covariances C
1
and C
2
, the linear
combination of features (i.e., the LDA-projected features using vector w) have means
~
1
= w
T

1
and ~
2
= w
T

2
. The maximum separation occurs when S is maximized,
i.e., that the projected means are as far from each other as possible, and that the scatters
113
(variances) S
1
and S
2
of the projected classes are as small as possible. Therefore, the
problem becomes:
max
w
S =
(~
2
~
1
)
2
S
1
2
+ S
2
2
=
w
T
(
2

1
)(
2

1
)
T
w
w
T
(S
1
+ S
2
)w
; (A.10)
which yields:
(
2

1
)(
2

1
)
T
w =(S
1
+ S
2
)w; (A.11)
where =
w
T
(
2

1
)(
2

1
)
T
w
w
T
(S
1
+S
2
)w
. The above is a generalized eigenvalue problem to which
the solution can be shown to be:
w = (S
1
+ S
2
)
1
(
2

1
): (A.12)
The projection of features from each class onto the line with direction w maximizes the
separation between projected classes.
Combining CSP and LDA, and based on Fisher's LDA above (equation A.12, the
classier function used in conjunction with CSP projection is:
f =
J
X
(j=1)

j
log(w
T
j
XX
T
w
j
) +
0
; (A.13)
114
where w
j
are the spatial lters determined earlier with CSP, and J is the total number
of lters used, and = [
1
:::
j
] and
0
are dened as follows based on LDA:
= (S
1
+ S
2
)
1
(
2

1
);
and

0
=:
1
2
(
1
+
2
)
(A.14)
Toy example illustrating LDA
Here we show an example of how LDA works with a simple two-feature classication
problem (Teknomo 2015). Suppose a factory measures the quality of their product (chip
rings) by curvature and diameter. Table A.1 lists these features for each of the available
training examples along with quality control result \Passed" or \Not passed". The pur-
pose of this problem is to set up a classier to automatically, using only the two features,
decide if the product will pass or not pass quality control. We will test this classier on
a the new chip ring shown in bold letters in Table A.1. We will apply LDA and visualize
how LDA projection of features from training examples makes them distinguishable. We
will then classify an unseen example as \passed" or \not passed."
Curvature Diameter Quality control result
2.95 6.63 Passed
2.53 7.79 Passed
3.57 5.65 Passed
3.16 5.47 Passed
2.58 4.46 Not passed
2.16 6.22 Not passed
3.27 3.52 Not passed
2.81 5.46 (to be classied)
Table A.1: LDA example: product dimensions and quality control result
115
We rst show the features in this example as a plot of diameter versus curvature in
Figure A.1. The blue x's belong to class 1, or \Passed, " while the red o's come from class
2, or \Not passed." We also show the test example as a black triangle, which will not
be used in training the classier but will be classied after training. Looking at the raw
features one may not easily classify the test example as "Passed" or "Not passed." We
now apply LDA projection using the seven training examples and compare the projected
features to the raw features.
curvature
4
5
6
7
8
diameter
2 2.4 2.8 3.2 3.6
Figure A.1: Features of the training set (x: "passed" and o: "not passed") and test set
(triangle) before projection.
116
30 35 40 45 50 55 60 65
F
1
35
40
45
50
55
60
F
2
F
1
= F
2
class 1: “passed”
class 2: “not passed”
Figure A.2: Projected features with LDA
Figure A.2 shows the LDA projected features, on projected axes F
1
and F
2
, along
with the line that separates the two classes. It is now clear that points to the left of the
separation line belong to class \Not passed," and those to the right of the line come from
class \Passed." Using the same reasoning, visually, the test example can be classied as
class "Passed" since its LDA projection is to the right of the separation line. Without
visual inspection, the test example can also be classied using the classier function:
f =F
2
F
1
= 4646:5 =0:5< 0, therefore, this example belongs to class 1: \passed".
117
A.3 Classication of simulated EEG example
In this section we simulate a training and classication problem based on two channels
for illustration purposes (inspired by (Blankertz et al. 2008)). Using two channels makes
it easier to visualize the data from each class and observe how they vary across channels
before and after CSP ltering.
First, we create two randomly-generated datasets X
1
andX
2
, from two 2D Gaussian
distributions, each containing signals from two channels c
1
andc
2
as correlated variables
(Figure A.3a). We assume that each dataset belongs to one class. Figure A.3a shows
X
1
has higher variance across channel c
1
, and lower variance across channel c
2
, and X
2
vice-versa. However,X
1
andX
2
are not completely distinct; there is a strong correlation
between the two channels. This correlation makes it dicult and most often inaccurate
to classify a new example belonging to either class based on variance distribution across
channels.
We now apply CSP ltering to the datasets and visualize how the spatially-ltered dis-
tributions change after transforming the raw data. Figure A.3b shows the same datasets,
mapped with CSP lters determined using the above method, and plotted again against
now surrogate channels s
1
and s
2
. The surrogate channels are the new directions across
which variances are maximally-discernible. As such,X
1
CSP
still has higher variance across
the new, surrogate channel s
1
, and X
2
CSP
has higher variance across s
2
. However, the
training datasets are now simultaneously de-correlated. In other words, the variance of
the projected signal X
1
CSP
is maximized across the S
1
axis, while being simultaneously
minimized across the s
2
axis, and vice-versa for X
2
CSP
. This simultaneous decorrelation
118
between classes makes it much easier and more accurate to predict the class of a new,
unseen dataset after projecting it with the same spatial lter coecients found earlier
using the training data.
The maximized dierence in variance between classes may also be visualized across
time. Figure A.4a and b shows the same raw data from Figure A.3a, now plotted against
time samples and from channelsc
1
andc
2
, respectively. It can be seen that there is a slight
dierence in the variance (or overall amplitude) between the two classes: X
1
has higher
variance in thec
1
direction, andX
2
has higher variance in thec
2
direction. Now consider
the same data shown in Figure A.4c and d, where X
1
CSP
=W
T
X
1
and X
2
CSP
=W
T
X
2
were mapped using the spatial lters determined earlier. Compared to the original data,
it is clear that the transformed dataset contains much better discernibility based on
variance or amplitude.
119
-15 -10 -5 0 5 10 15
c
1
-15
-10
-5
0
5
10
15
c
2
Raw signals: before CSP
-4 -2 2 4
-4
-2
0
2
4
s
2
Spatially filtered signals: after CSP
a b
X
1
X
2
X
1
CSP
X
2
CSP
s
1
0
Figure A.3: Raw (a) and spatially ltered (b) simulated EEG signals.
120
0 100 200 300 400
time
-15
-5
5
15
amplitude
Raw signals, c
1
direction
X
1
X
2
0 100 200 300 400
time
-15
-5
5
15
amplitude
Raw signals, c
2
direction
X
1
X
2
0 100 200 300 400
time
-4
-2
0
2
amplitude
Spatially filtered signals, s
1
direction
X
1
CSP
X
2
CSP
0 100 200 300 400
time
-4
-2
0
2
4
amplitude
Spatially filtered signals, s
2
direction
X
1
CSP
X
2
CSP
Figure A.4: Raw and projected time series of simulated EEG signals.
We now attempt to predict the class of a test set that was not used in the training. We
create a new dataset from the same gaussian distribution as the training set of class 1, and
use the CSP and LDA methods to classify the unseen test case into class 1 or 2. Figure A.5
shows the results of the prediction after CSP spatial ltering and LDA classication. The
F
1
and F
2
axes represent the lines representing the projection onto class 1 and class 2,
respectively. Each blue dot represents one of the eight training sets used for classier
training. The dashed line is that of F
1
= F
2
, and separates the two classes such that
121
any projection (F
1
;F
2
)that falls to the right of the separation is classied as class 1, and
any projection that falls to its left is classied as class 2. The alternative of non-visual
classication is to calculate f = F
2
F
1
: if f < 0, the test set belongs to class 1, and
if f > 0, it belongs to class 2. In this case, the red x is a new test set to be classied.
Visually, it is easy to see the red x belongs to class 1, as it is placed to the right of the line
dividing the two classes. We can calculate f in this case as: f =129 60 =189< 0,
therefore, the set belongs to class 1.
-150 -100 -50 50 100
F
1
-120
-80
-40
0
40
F
2
80
0
F
1
= F
2
class 1
class 2
Figure A.5: Class prediction of simulated EEG signals. Visually, the test set (red x)
belongs to class 1, because it falls to the right of the separation line F
1
=F
2
.
122
Bibliography
An, K., Ueba, Y., Chao, E., Cooney, W. & Linscheid, R. (1983), `Tendon excursion and
moment arm of index nger muscles', Journal of Biomechanics 16(6), 419{425.
Anderson, F. C. & Pandy, M. G. (2001), `Dynamic optimization of human walking',
Journal of biomechanical engineering 123(5), 381{390.
Arnold, E. M., Ward, S. R., Lieber, R. L. & Delp, S. L. (2010), `A model of the lower limb
for analysis of human movement', Annals of biomedical engineering 38(2), 269{279.
Asatryan, D. G. & Feldman, A. G. (1965), `Functional tuning of the nervous system with
control of movement or maintenance of a steady posture: I. mechanographic analysis
of the work of the joint or execution of a postural task', Biophysics 10(5), 925{934.
Asavasopon, S., Rana, M., Kirages, D. J., Yani, M. S., Fisher, B. E., Hwang, D. H.,
Lohman, E. B., Berk, L. S. & Kutch, J. J. (2014), `Cortical activation associated
with muscle synergies of the human male pelvic
oor', The Journal of Neuroscience
34(41), 13811{13818.
Babikian, S., Kanso, E. & Kutch, J. J. (2017), `Cortical activity predicts good variation
in human motor output', Experimental Brain Research pp. 1{9.
Berkinblit, M., Feldman, A. & Fukson, O. (1986), `Adaptability of innate motor patterns
and motor control mechanisms', Behavioral and Brain Sciences 9(04), 585{599.
Bernstein, N. (1967), Co-ordination and regulation of movements, Pergamon Press, Ox-
ford.
Bizzini, M. & Mannion, A. F. (2003), `Reliability of a new, hand-held device for assessing
skeletal muscle stiness', Clinical Biomechanics 18(5), 459{461.
Blankertz, B., Dornhege, G., Lemm, S., Krauledat, M., Curio, G. & M uller, K.-R. (2006),
`The berlin brain-computer interface: Machine learning based detection of user spe-
cic brain states', J. UCS 12(6), 581{607.
Blankertz, B., Kawanabe, M., Tomioka, R., Hohlefeld, F., M uller, K.-r. & Nikulin,
V. V. (2007), Invariant common spatial patterns: Alleviating nonstationarities in
brain-computer interfacing, in `Advances in neural information processing systems',
pp. 113{120.
123
Blankertz, B., Tomioka, R., Lemm, S., Kawanabe, M. & Muller, K.-R. (2008), `Optimizing
spatial lters for robust eeg single-trial analysis', Signal Processing Magazine, IEEE
25(1), 41{56.
Boscain, U. & Chitour, Y. (2002), `On the minimum time problem for driftless left-
invariant control systems on so (3)', Commun. Pure Appl. Anal 1, 285{312.
Buchanan, T. S., Rovai, G. P. & Rymer, W. Z. (1989), `Strategies for muscle activation
during isometric torque generation at the human elbow', Journal of Neurophysiology
62(6), 1201{1212.
Burdet, E., Osu, R., Franklin, D. W., Milner, T. E. & Kawato, M. (2001), `The central
nervous system stabilizes unstable dynamics by learning optimal impedance', Nature
414(6862), 446{449.
Burdet, E., Osu, R., Franklin, D., Yoshioka, T., Milner, T. & Kawato, M. (2000), `A
method for measuring endpoint stiness during multi-joint arm movements', Journal
of biomechanics 33(12), 1705{1709.
Carson, R. (2005), `Neural pathways mediating bilateral interactions between the upper
limbs', Brain Research Reviews 49(3), 641{662.
Chen, J., Tian, J., Iwasaki, T. & Friesen, W. O. (2011), `Mechanisms underlying rhythmic
locomotion: dynamics of muscle activation', The Journal of experimental biology
214(11), 1955{1964.
Chen, R., Yung, D. & Li, J.-Y. (2003), `Organization of ipsilateral excitatory and
inhibitory pathways in the human motor cortex', Journal of neurophysiology
89(3), 1256{1264.
Chow, C. & Jacobson, D. (1971), `Studies of human locomotion via optimal program-
ming', Mathematical Biosciences 10(3-4), 239{306.
Chuang, L.-l., Wu, C.-y. & Lin, K.-c. (2012), `Reliability, validity, and responsiveness of
myotonometric measurement of muscle tone, elasticity, and stiness in patients with
stroke', Archives of physical medicine and rehabilitation 93(3), 532{540.
Churchland, M. M., Byron, M. Y., Ryu, S. I., Santhanam, G. & Shenoy, K. V. (2006),
`Neural variability in premotor cortex provides a signature of motor preparation',
The Journal of neuroscience 26(14), 3697{3712.
Clauw, D. & Chrousos, G. (1998), `Chronic pain and fatigue syndromes: overlapping
clinical and neuroendocrine features and potential pathogenic mechanisms', Neu-
roimmunomodulation 4(3), 134{153.
Collins, J. (1995), `The redundant nature of locomotor optimization laws', Journal of
biomechanics 28(3), 251{267.
Cooke, J. (1980), `The organization of simple, skilled movements', Tutorials in motor
behavior pp. 199{212.
124
Darling, W. G., Cole, K. J. & Miller, G. F. (1994), `Coordination of index nger move-
ments', Journal of biomechanics 27(4), 479{491.
Davy, D. & Audu, M. (1987), `A dynamic optimization technique for predicting muscle
forces in the swing phase of gait', Journal of biomechanics 20(2), 187{201.
de Freitas, S. M. S. F. & Scholz, J. P. (2010), `A comparison of methods for identifying
the jacobian for uncontrolled manifold variance analysis', Journal of biomechanics
43(4), 775{777.
De Luca, A., Oriolo, G. & Samson, C. (1998), Feedback control of a nonholonomic car-like
robot, Springer, book section 4, pp. 171{253.
Deeny, S. P., Hau
er, A. J., Saer, M. & Hateld, B. D. (2009), `Electroencephalographic
coherence during visuomotor performance: a comparison of cortico-cortical commu-
nication in experts and novices', Journal of motor behavior 41(2), 106{116.
Dornhege, G. (2007), Toward brain-computer interfacing, MIT press.
Dornhege, G., Blankertz, B., Krauledat, M., Losch, F., Curio, G. & M uller, K.-R. (2006),
`Optimizing spatio-temporal lters for improving brain-computer interfacing', Ad-
vances in Neural Information Processing Systems 18, 315.
Evans, C. & Baker, S. N. (2003), `Task?dependent intermanual coupling of 8?hz disconti-
nuities during slow nger movements', European Journal of Neuroscience 18(2), 453{
456.
Feldman, A. (1966), `Functional tuning of the nervous system with control of movement
or maintenance of a steady posture. ii. controllable parameters of the muscle', Bio-
physics 11(3), 565{578.
Feldman, A. G. & Levin, M. F. (2009), The Equilibrium-Point Hypothesis - Past, Present
and Future, Vol. 629, Springer, pp. 699{726.
Ferbert, A., Priori, A., Rothwell, J., Day, B., Colebatch, J. & Marsden, C. (1992), `In-
terhemispheric inhibition of the human motor cortex', The Journal of physiology
453(1), 525{546.
Fisher, R. A. (1936), `The use of multiple measurements in taxonomic problems', Annals
of eugenics 7(2), 179{188.
Flash, T. & Hogan, N. (1985), `The coordination of arm movements: an experimentally
conrmed mathematical model', Journal of neuroscience 5(7), 1688{1703.
Fling, B. W., Peltier, S., Bo, J., Welsh, R. C. & Seidler, R. D. (2011), `Age dierences in
interhemispheric interactions: callosal structure, physiological function, and behav-
ior', Frontiers in neuroscience 5, 38.
Fling, B. W. & Seidler, R. D. (2012), `Task-dependent eects of interhemispheric inhibi-
tion on motor control', Behavioural brain research 226(1), 211{217.
125
Frecker, M. & Snyder, A. J. (2005), `Surgical robotics: multifunctional end eectors for
robotic surgery', Operative Techniques in General Surgery 7(4), 165{169.
Fukunaga, R. (1990), `Statistical pattern recognition'.
Giszter, S. F. (2015), `Motor primitives|new data and future questions', Current opinion
in neurobiology 33, 156{165.
Giszter, S. F., Mussa-Ivaldi, F. A. & Bizzi, E. (1993), `Convergent force elds organized
in the frog's spinal cord', The journal of neuroscience 13(2), 467{491.
Giszter, S., Patil, V. & Hart, C. (2007), `Primitives, premotor drives, and pattern gen-
eration: a combined computational and neuroethological perspective', Progress in
brain research 165, 323{346.
Grefkes, C., Eickho, S. B., Nowak, D. A., Dafotakis, M. & Fink, G. R. (2008), `Dy-
namic intra-and interhemispheric interactions during unilateral and bilateral hand
movements assessed with fmri and dcm', Neuroimage 41(4), 1382{1394.
Grith, J. (1968), `Mathematics of cellular control processes i. negative feedback to one
gene', Journal of theoretical biology 20(2), 202{208.
Grillner, S. & Wallen, P. (1985), `Central pattern generators for locomotion, with special
reference to vertebrates', Annual review of neuroscience 8(1), 233{261.
Gross, J., Timmermann, L., Kujala, J., Dirks, M., Schmitz, F., Salmelin, R. & Schnitzler,
A. (2002), `The neural basis of intermittent motor control in humans', Proceedings
of the National Academy of Sciences 99(4), 2299{2302.
Hamill, J., Heiderscheit, B. C. & Pollard, C. D. (2005), `Gender dierences in lower
extremity coupling variability during an unanticipated cutting maneuver', Journal
of applied biomechanics 21(2).
Hamill, J., Palmer, C. & Van Emmerik, R. E. (2012), `Coordinative variability and overuse
injury', BMC Sports Science, Medicine and Rehabilitation 4(1), 45.
Hamill, J., van Emmerik, R. E., Heiderscheit, B. C. & Li, L. (1999), `A dynamical systems
approach to lower extremity running injuries', Clinical Biomechanics 14(5), 297{308.
Harris-Love, M. L., Perez, M. A., Chen, R. & Cohen, L. G. (2007), `Interhemispheric
inhibition in distal and proximal arm representations in the primary motor cortex',
Journal of Neurophysiology 97(3), 2511{2515.
Hastie, T., Tibshirani, R., Friedman, J. & Franklin, J. (2005), `The elements of statistical
learning: data mining, inference and prediction', The Mathematical Intelligencer
27(2), 83{85.
Hatze, H. & Buys, J. (1977), `Energy-optimal controls in the mammalian neuromuscular
system', Biological cybernetics 27(1), 9{20.
126
Heinrichs-Graham, E. & Wilson, T. W. (2016), `Is an absolute level of cortical beta
suppression required for proper movement? magnetoencephalographic evidence from
healthy aging', NeuroImage 134, 514{521.
Hill, A. (1938), `The heat of shortening and the dynamic constants of muscle', Proceedings
of the Royal Society of London B: Biological Sciences 126(843), 136{195.
Hogan, N. (1984a), `Adaptive control of mechanical impedance by coactivation of antag-
onist muscles', Automatic Control, IEEE Transactions on 29(8), 681{690.
Hogan, N. (1984b), `An organizing principle for a class of voluntary movements', Journal
of Neuroscience 4(11), 2745{2754.
Hogan, N. (1985), `Impedance control: An approach to manipulation: Part ii - imple-
mentation', Journal of dynamic systems, measurement, and control 107(1), 8{16.
Hopeld, J. J. (1982), `Neural networks and physical systems with emergent col-
lective computational abilities', Proceedings of the national academy of sciences
79(8), 2554{2558.
Ibey, R., Bolton, D., Buick, A., Staines, W. & Carson, R. (2015), `Interhemispheric
inhibition of corticospinal projections to forearm muscles', Clinical Neurophysiology
126(10), 1934{1940.
Ivaldi, F. M., Morasso, P. & Zaccaria, R. (1988), `Kinematic networks', Biological cyber-
netics 60(1), 1{16.
Jing, F. & Alben, S. (2013), `Optimization of two-and three-link snakelike locomotion',
Physical Review E 87(2), 022711.
Joshua, M. & Lisberger, S. G. (2014), `A framework for using signal, noise, and varia-
tion to determine whether the brain controls movement synergies or single muscles',
Journal of neurophysiology 111(4), 733{745.
Jung, S.-Y., Kang, S.-k., Lee, M.-J. & Moon, I. (2015), Design of robotic hand with
tendon-driven three ngers, in `Control, Automation and Systems, 2007. ICCAS'07.
International Conference on', IEEE, pp. 83{86.
Kanso, E. (2009), `Swimming due to transverse shape deformations', Journal of Fluid
Mechanics 631, 127{148.
Kanso, E., Marsden, J. E., Rowley, C. W. & Melli-Huber, J. (2005), `Locomotion of
articulated bodies in a perfect
uid', Journal of Nonlinear Science 15(4), 255{289.
Kilpatrick, L. A., Kutch, J. J., Tillisch, K., Nalibo, B. D., Labus, J. S., Jiang, Z., Farmer,
M. A., Apkarian, A. V., Mackey, S. & Martucci, K. T. (2014), `Alterations in resting
state oscillations and connectivity in sensory and motor networks in women with
interstitial cystitis/painful bladder syndrome', The Journal of urology 192(3), 947{
955.
127
Kobayashi, M. & Pascual-Leone, A. (2003), `Transcranial magnetic stimulation in neu-
rology', The Lancet Neurology 2(3), 145{156.
Koles, Z. J. (1991), `The quantitative extraction and topographic mapping of the abnor-
mal components in the clinical eeg', Electroencephalography and clinical Neurophys-
iology 79(6), 440{447.
Kouzaki, M. & Shinohara, M. (2006), `The frequency of alternate muscle activity is
associated with the attenuation in muscle fatigue', Journal of applied physiology
101(3), 715{720.
Krishnamoorthy, V., Latash, M. L., Scholz, J. P. & Zatsiorsky, V. M. (2003), `Muscle
synergies during shifts of the center of pressure by standing persons', Experimental
brain research 152(3), 281{292.
Kujirai, T., Caramia, M., Rothwell, J. C., Day, B., Thompson, P., Ferbert, A., Wroe, S.,
Asselman, P. & Marsden, C. D. (1993), `Corticocortical inhibition in human motor
cortex', The Journal of physiology 471(1), 501{519.
Kuo, A. D. & Zajac, F. E. (1993), `A biomechanical analysis of muscle strength as a
limiting factor in standing posture', Journal of Biomechanics 26, 137{150.
Kutch, J. J. & Valero-Cuevas, F. J. (2012), `Challenges and new approaches to prov-
ing the existence of muscle synergies of neural origin', PLoS computational biology
8(5), e1002434.
Kutch, J. J., Yani, M. S., Asavasopon, S., Kirages, D. J., Rana, M., Cosand, L., Labus,
J. S., Kilpatrick, L. A., Ashe-McNalley, C. & Farmer, M. A. (2015), `Altered resting
state neuromotor connectivity in men with chronic prostatitis/chronic pelvic pain
syndrome: A mapp: Research network neuroimaging study', NeuroImage: Clinical
8, 493{502.
Latash, M. & Huang, X. (2015), `Neural control of movement stability: lessons from
studies of neurological patients', Neuroscience 301, 39{48.
Latash, M. L. (2008), Synergy, Oxford University Press.
Latash, M. L. & Anson, J. G. (2006), `Synergies in health and disease: relations to
adaptive changes in motor coordination', Physical therapy 86(8), 1151{1160.
Latash, M. L., Li, Z.-M. & Zatsiorsky, V. M. (1998), `A principle of error compensation
studied within a task of force production by a redundant set of ngers', Experimental
brain research 122(2), 131{138.
Latash, M. L., Scholz, J. F., Danion, F. & Sch oner, G. (2001), `Structure of motor
variability in marginally redundant multinger force production tasks', Experimental
brain research 141(2), 153{165.
Latash, M. L., Yarrow, K. & Rothwell, J. C. (2003), `Changes in nger coordination and
responses to single pulse tms of motor cortex during practice of a multinger force
production task', Experimental brain research 151(1), 60{71.
128
Lee, D., Port, N. L., Kruse, W. & Georgopoulos, A. P. (1998), `Variability and correlated
noise in the discharge of neurons in motor and parietal areas of the primate cortex',
Journal of Neuroscience 18(3), 1161{1170.
Lee, Y.-T., Choi, H.-R., Chung, W.-K. & Youm, Y. (1994), `Stiness control of a coupled
tendon-driven robot hand', Control Systems, IEEE 14(5), 10{19.
Lemm, S., Blankertz, B., Curio, G. & Muller, K.-R. (2005), `Spatio-spectral lters for
improving the classication of single trial eeg', IEEE Transactions on Biomedical
Engineering 52(9), 1541{1548.
Lipsitz, L. A. (2002), `Dynamics of stability the physiologic basis of functional health
and frailty', The Journals of Gerontology Series A: Biological Sciences and Medical
Sciences 57(3), B115{B125.
Lotka, A. J. (1925), `Elements of physical biology'.
McMahon, T. A. (1984), `Muscles, re
exes, and locomotion', Princeton, Princeton .
Medina, J. F. & Lisberger, S. G. (2007), `Variation, signal, and noise in cerebellar sensory-
motor processing for smooth-pursuit eye movements', The Journal of neuroscience
27(25), 6832{6842.
Morishita, I. & Yajima, A. (1972), `Analysis and simulation of networks of mutually
inhibiting neurons', Kybernetik 11(3), 154{165.
Morishita, T., Uehara, K. & Funase, K. (2012), `Changes in interhemispheric inhibition
from active to resting primary motor cortex during a ne-motor manipulation task',
Journal of neurophysiology 107(11), 3086{3094.
M uller-Gerking, J., Pfurtscheller, G. & Flyvbjerg, H. (1999), `Designing optimal spatial
lters for single-trial eeg classication in a movement task', Clinical neurophysiology
110(5), 787{798.
Mussa-Ivaldi, F. A. & Hogan, N. (1991), `Integrable solutions of kinematic redundancy via
impedance control', The International Journal of Robotics Research 10(5), 481{491.
Mustalampi, S., H

'akkinen, A., Kautiainen, H., Weir, A. & Ylinen, J. (2013), `Respon-
siveness of muscle tone characteristics to progressive force production', The Journal
of Strength and Conditioning Research 27(1), 159{165.
Mutha, P. K. & Sainburg, R. L. (2009), `Shared bimanual tasks elicit bimanual re
exes
during movement', Journal of neurophysiology 102(6), 3142{3155.
Myre, C. & Woodward, D. (1993), `Bistability, switches and working memory in a two-
neuron inhibitory-feedback model', Biological cybernetics 68(5), 441{449.
Nakano, E., Imamizu, H., Osu, R., Uno, Y., Gomi, H., Yoshioka, T. & Kawato, M.
(1999), `Quantitative examinations of internal representations for arm trajectory
planning: minimum commanded torque change model', Journal of Neurophysiology
81(5), 2140{2155.
129
Ni, Z., Gunraj, C., Nelson, A. J., Yeh, I.-J., Castillo, G., Hoque, T. & Chen, R. (2009),
`Two phases of interhemispheric inhibition between motor related cortical areas and
the primary motor cortex in human', Cerebral Cortex 19(7), 1654{1665.
Nichols, T. (1994), `A biomechanical perspective on spinal mechanisms of coordinated
muscular action: an architecture principle', Cells Tissues Organs 151(1), 1{1.
Pandy, M., Garner, B. & Anderson, F. (1995), `Optimal control of non-ballistic ivius-
cuiar iviovements: A constraint-based performance criterion for rising from a cliair',
Journal of Biomechanical Engineering 117, 15.
Pfurtscheller, G. & Aranibar, A. (1979), `Evaluation of event-related desynchronization
(erd) preceding and following voluntary self-paced movement', Electroencephalogra-
phy and clinical neurophysiology 46(2), 138{146.
Pfurtscheller, G. & Neuper, C. (1997), `Motor imagery activates primary sensorimotor
area in humans', Neuroscience letters 239(2), 65{68.
Pfurtscheller, G., Neuper, C., Flotzinger, D. & Pregenzer, M. (1997), `Eeg-based discrim-
ination between imagination of right and left hand movement', Electroencephalogra-
phy and clinical Neurophysiology 103(6), 642{651.
Popovic, D., Stein, R. B., Oguztoreli, M. N., Lebiedowska, M. & Jonic, S. (1999), `Opti-
mal control of walking with functional electrical stimulation: a computer simulation
study', IEEE Transactions on Rehabilitation Engineering 7(1), 69{79.
Prilutsky, B. I. & Zatsiorsky, V. M. (2002), `Optimization-based models of muscle coor-
dination', Exercise and sport sciences reviews 30(1), 32.
Ramoser, H., Muller-Gerking, J. & Pfurtscheller, G. (2000), `Optimal spatial ltering of
single trial eeg during imagined hand movement', Rehabilitation Engineering, IEEE
Transactions on 8(4), 441{446.
Rana, M., Yani, M. S., Asavasopon, S., Fisher, B. E. & Kutch, J. J. (2015), `Brain con-
nectivity associated with muscle synergies in humans', The Journal of Neuroscience
35(44), 14708{14716.
Rasmussen, J., Damsgaard, M. & Voigt, M. (2001), `Muscle recruitment by the min/max
criterion|a comparative numerical study', Journal of biomechanics 34(3), 409{415.
R atsep, T. o. & Asser, T. (2011), `Changes in viscoelastic properties of skeletal muscles
induced by subthalamic stimulation in patients with parkinson's disease', Clinical
Biomechanics 26(2), 213{217.
Rosenbaum, D. A., Loukopoulos, L. D., Meulenbroek, R. G., Vaughan, J. & Engelbrecht,
S. E. (1995), `Planning reaches by evaluating stored postures', Psychological review
102(1), 28.
Saltiel, P., Wyler-Duda, K., D'Avella, A., Tresch, M. C. & Bizzi, E. (2001), `Muscle
synergies encoded within the spinal cord: evidence from focal intraspinal nmda
iontophoresis in the frog', Journal of Neurophysiology 85(2), 605{619.
130
Scholz, J. P. & Sch oner, G. (1999), `The uncontrolled manifold concept: identifying
control variables for a functional task', Experimental brain research 126(3), 289{
306.
Seay, J. F., Haddad, J. M., Van Emmerik, R. E. & Hamill, J. (2006), `Coordination
variability around the walk to run transition during human locomotion', MOTOR
CONTROL-CHAMPAIGN- 10(2), 178.
Seay, J. F., Van Emmerik, R. E. & Hamill, J. (2011), `Low back pain status aects pelvis-
trunk coordination and variability during walking and running', Clinical Biomechan-
ics 26(6), 572{578.
Seif-Naraghi, A. H. & Winters, J. M. (1990), Optimized strategies for scaling goal-directed
dynamic limb movements, Springer, pp. 312{334.
Shadmehr, R. (1998), `The equilibrium point hypothesis for control of movement', Balti-
more, MD: Department of Biomedical Engineering, Johns Hopkins University .
Shapere, A. & Wilczek, F. (1989), `Geometry of self-propulsion at low reynolds number',
Journal of Fluid Mechanics 198, 557{585.
Shibasaki, H. & Hallett, M. (2006), `What is the bereitschaftspotential?', Clinical neuro-
physiology 117(11), 2341{2356.
Shim, J. K., Kim, S. W., Oh, S. J., Kang, N., Zatsiorsky, V. M. & Latash, M. L. (2005),
`Plastic changes in interhemispheric inhibition with practice of a two-hand force
production task: a transcranial magnetic stimulation study', Neuroscience letters
374(2), 104{108.
Simaan, N., Xu, K., Wei, W., Kapoor, A., Kazanzides, P., Taylor, R. & Flint, P. (2009),
`Design and integration of a telerobotic system for minimally invasive surgery of the
throat', The International journal of robotics research 28(9), 1134{1153.
Singh, T., Varadhan, S., Zatsiorsky, V. M. & Latash, M. L. (2010), `Fatigue and motor
redundancy: adaptive increase in nger force variance in multi-nger tasks', Journal
of neurophysiology 103(6), 2990{3000.
Smeets, J. B. & Brenner, E. (1999), `A new view on grasping', Motor control 3(3), 237{
271.
Softky, W. R. & Koch, C. (1993), `The highly irregular ring of cortical cells is inconsistent
with temporal integration of random epsps', The Journal of Neuroscience 13(1), 334{
350.
Stein, R. B., Gossen, E. R. & Jones, K. E. (2005), `Neuronal variability: noise or part of
the signal?', Nature Reviews Neuroscience 6(5), 389{397.
Teel, A. R., Murray, R. M. & Walsh, G. C. (1995), `Non-holonomic control systems: from
steering to stabilization with sinusoids', International Journal of Control 62(4), 849{
870.
131
Teknomo, K. (2015), `Discriminant analysis tutorial.', http://people.revoledu.com/
kardi/tutorial/.
Todorov, E. (2004), `Optimality principles in sensorimotor control', Nature neuroscience
7(9), 907{915.
Todorov, E. & Jordan, M. I. (1998), `Smoothness maximization along a predened path
accurately predicts the speed proles of complex arm movements', Journal of Neu-
rophysiology 80(2), 696{714.
Todorov, E. & Jordan, M. I. (2002), `Optimal feedback control as a theory of motor
coordination', Nature neuroscience 5(11), 1226{1235.
Trotta, L., Bullinger, E. & Sepulchre, R. (2012), `Global analysis of dynamical decision-
making models through local computation around the hidden saddle', PloS one
7(3), e33110.
Turella, L., Tucciarelli, R., Oosterhof, N. N., Weisz, N., Rumiati, R. & Lingnau, A. (2016),
`Beta band modulations underlie action representations for movement planning',
NeuroImage .
Uno, Y., Kawato, M. & Suzuki, R. (1989), `Formation and control of optimal trajectory
in human multijoint arm movement', Biological cybernetics 61(2), 89{101.
Valero-Cuevas, F. J. (2009), A mathematical approach to the mechanical capabilities of
limbs and ngers, in `Progress in Motor Control', Vol. 629 of Advances in Experi-
mental Medicine and Biology, Springer US, pp. 619{633.
Valero-Cuevas, F. J., Yi, J.-W., Brown, D., McNamara, R. V., Paul, C. & Lipson, H.
(2007), `The tendon network of the ngers performs anatomical computation at
a macroscopic scale', Biomedical Engineering, IEEE Transactions on 54(6), 1161{
1166.
Valero-Cuevas, F. J., Zajac, F. E. & Burgar, C. G. (1998), `Large index-ngertip forces are
produced by subject-independent patterns of muscle excitation', Journal of biome-
chanics 31(8), 693{703.
Vallbo, A. & Wessberg, J. (1993), `Organization of motor output in slow nger movements
in man', The Journal of physiology 469(1), 673{691.
Volterra, V. & Brelot, M. (1931), Le cons sur la th eorie math ematique de la lutte pour la
vie, Vol. 1, Gauthier-Villars Paris.
Wessberg, J. & Vallbo, A. (1995), `Coding of pulsatile motor output by human muscle
aerents during slow nger movements', The Journal of physiology 485(Pt 1), 271{
282.
Wessberg, J. & Vallbo, A. (1996), `Pulsatile motor output in human nger movements is
not dependent on the stretch re
ex', The Journal of physiology 493(Pt 3), 895{908.
132
Whitney, D. E. (1969), `Resolved motion rate control of manipulators and human pros-
theses', IEEE Transactions on man-machine systems .
Williams, E. R., Soteropoulos, D. S. & Baker, S. N. (2009), `Coherence between mo-
tor cortical activity and peripheral discontinuities during slow nger movements',
Journal of neurophysiology 102(2), 1296{1309.
Williams, T. L. (2010), `A new model for force generation by skeletal muscle, incorporat-
ing work-dependent deactivation', The Journal of experimental biology 213(4), 643{
650.
Woodworth, D., Mayer, E., Leu, K., Ashe-McNalley, C., Nalibo, B. D., Labus, J. S.,
Tillisch, K., Kutch, J. J., Farmer, M. A. & Apkarian, A. V. (2015), `Unique mi-
crostructural changes in the brain associated with urological chronic pelvic pain
syndrome (ucpps) revealed by diusion tensor mri, super-resolution track density
imaging, and statistical parameter mapping: A mapp network neuroimaging study',
PloS one 10(10), e0140250.
Yang, J.-F. & Scholz, J. (2005), `Learning a throwing task is associated with dierential
changes in the use of motor abundance', Experimental brain research 163(2), 137{
158.
133

Abstract (if available)

Abstract It remains unknown how the human brain coordinates the many biomechanical degrees of freedom (e.g. muscle activations) to generate movements. Many theories explaining motor coordination have been proposed, ranging from optimizing for a unique muscle activation pattern, to exploring the subspace of plausible activation patterns that equally accomplish the same task. This dissertation focuses on investigating the feasible solution space of muscle activation patterns that produce the same movement or force production. We use both computational and experimental approaches to understand the mechanisms and neural origins of varied patterns that achieve the same task. ❧ We first present a mathematical model of a tendon-driven limb controlled using energy minimization in muscles

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Transitions in the collective behavior of microswimmers confined to microfluidic chips

PDF

Decoding limb movement from posterior parietal cortex in a realistic task

PDF

Development and control of biologically-inspired robots driven by artificial muscles

PDF

Neuromuscular dynamics in the context of motor redundancy

PDF

Passive flight in density-stratified fluid environments

PDF

A theoretical investigation in multi-dimensional and non-thermal plasma effects on combustion processes and pollutant remediation

PDF

Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response

PDF

Synchronization in carpets of model cilia: numerical simulations and continuum theory

PDF

insideOut: Estimating joint angles in tendon-driven robots using Artificial Neural Networks and non-collocated sensors

PDF

The footprint of pain: investigating persistence of altered trunk control in recurrent low back pain

PDF

Sensory acquisition for emergent body representations in neuro-robotic systems

PDF

Model-based studies of control strategies for noisy, redundant musculoskeletal systems

PDF

Task-dependent modulation of corticomuscular coherence during dexterous manipulation

PDF

Bio-inspired tendon-driven systems: computational analysis, optimization, and hardware implementation

PDF

Influence of endogenous and exogenous estrogen on skeletal muscle damage and systemic inflammation after electrically stimulated contraction in young women

PDF

Evaluating sensing and control in underwater animal behaviors

PDF

Central nervous system generation of flexor bias in the control of leg movements during locomotor-related behavior in chick embryos

PDF

Physics-based and data-driven models for bio-inspired flow sensing and motion planning

PDF

A stochastic Markov chain approach for tennis: Monte Carlo simulation and modeling

PDF

Detection and decoding of cognitive states from neural activity to enable a performance-improving brain-computer interface

Asset Metadata

Creator Babikian, Sarine (author)

Core Title Coordination in multi-muscle systems: from control theory to cortex

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Publication Date 07/10/2017

Defense Date 04/24/2017

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag brain,cortex,movement variability,muscle,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kanso, Eva (committee chair), Kutch, Jason J. (committee member), Newton, Paul (committee member), Pahlevan, Niema (committee member)

Creator Email sarine.babikian@gmail.com,sbabikia@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-397286

Unique identifier UC11263258

Identifier etd-BabikianSa-5494.pdf (filename),usctheses-c40-397286 (legacy record id)

Legacy Identifier etd-BabikianSa-5494.pdf

Dmrecord 397286

Document Type Dissertation

Rights Babikian, Sarine

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA