UC Davis
UC Davis Previously Published Works
Title
Muscle Synergies Facilitate Computational Prediction of Subject-Specific Walking Motions.
Permalink
https://escholarship.org/uc/item/84h3r42d
Journal
Frontiers in bioengineering and biotechnology, 4(OCT)
ISSN
2296-4185
Authors
Meyer, Andrew J
Eskinazi, Ilan
Jackson, Jennifer N
et al.
Publication Date
2016
DOI
10.3389/fbioe.2016.00077
Peer reviewed
eScholarship.org Powered by the California Digital Library
University of California
October 2016 | Volume 4 | Article 771
METHODS
published: 13 October 2016
doi: 10.3389/fbioe.2016.00077
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org
Edited by:
Ramana Vinjamuri,
Stevens Institute of
Technology, USA
Reviewed by:
Laurent Simon,
New Jersey Institute of
Technology, USA
Jonathan B. Shute,
University of Florida, USA
*Correspondence:
Benjamin J. Fregly
fregly@u.edu
Specialty section:
This article was submitted
to Bionics and Biomimetics,
a section of the journal
Frontiers in Bioengineering
and Biotechnology
Received: 01June2016
Accepted: 21September2016
Published: 13October2016
Citation:
MeyerAJ, EskinaziI, JacksonJN,
RaoAV, PattenC and FreglyBJ
(2016) Muscle Synergies Facilitate
Computational Prediction of
Subject-Specic Walking Motions.
Front. Bioeng. Biotechnol. 4:77.
doi: 10.3389/fbioe.2016.00077
Muscle Synergies Facilitate
Computational Prediction of
Subject-Specic Walking Motions
Andrew J. Meyer
1
, Ilan Eskinazi
1
, Jennifer N. Jackson
2
, Anil V. Rao
1
, Carolynn Patten
3,4
and
Benjamin J. Fregly
1,2
*
1
Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL, USA,
2
Department of
Biomedical Engineering, University of Florida, Gainesville, FL, USA,
3
Department of Physical Therapy, University of Florida,
Gainesville, FL, USA,
4
Neural Control of Movement Lab, Malcom-Randall VA Medical Center, Gainesville, FL, USA
Researchers have explored a variety of neurorehabilitation approaches to restore normal
walking function following a stroke. However, there is currently no objective means for
prescribing and implementing treatments that are likely to maximize recovery of walking
function for any particular patient. As a rst step toward optimizing neurorehabilitation
effectiveness, this study develops and evaluates a patient-specic synergy-controlled
neuro musculoskeletal simulation framework that can predict walking motions for
an individual post-stroke. The main question we addressed was whether driving a
subject-specic neuromusculoskeletal model with muscle synergy controls (5 per leg) facil-
itates generation of accurate walking predictions compared to a model driven by muscle
activation controls (35 per leg) or joint torque controls (5 per leg). To explore this question,
we developed a subject-specic neuromusculoskeletal model of a single high-functioning
hemiparetic subject using instrumented treadmill walking data collected at the subject’s
self-selected speed of 0.5m/s. The model included subject-specic representations of
lower-body kinematic structure, foot–ground contact behavior, electromyography-driven
muscle force generation, and neural control limitations and remaining capabilities. Using
direct collocation optimal control and the subject-specic model, we evaluated the ability
of the three control approaches to predict the subject’s walking kinematics and kinetics at
two speeds (0.5 and 0.8m/s) for which experimental data were available from the subject.
We also evaluated whether synergy controls could predict a physically realistic gait period
at one speed (1.1m/s) for which no experimental data were available. All three control
approaches predicted the subject’s walking kinematics and kinetics (including ground
reaction forces) well for the model calibration speed of 0.5m/s. However, only activation
and synergy controls could predict the subject’s walking kinematics and kinetics well for
the faster non-calibration speed of 0.8m/s, with synergy controls predicting the new gait
period the most accurately. When used to predict how the subject would walk at 1.1m/s,
synergy controls predicted a gait period close to that estimated from the linear relationship
between gait speed and stride length. These ndings suggest that our neuromusculo-
skeletal simulation framework may be able to bridge the gap between patient-specic
muscle synergy information and resulting functional capabilities and limitations.
Keywords: biomechanics, computational neurorehabilitation, direct collocation optimal control, muscle synergy
analysis, neuromusculoskeletal modeling, predictive gait optimization
2
Meyer et al. Muscle Synergies Facilitate Walking Predictions
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2016 | Volume 4 | Article 77
INTRODUCTION
Roughly one in six people worldwide will suer a stroke at some
point in their lifetime, with ~15 million people experiencing
a stroke each year (
World Stroke Organization, 2016). Due to
improvements in acute stroke management, the majority of these
individuals will survive their initial stroke, which has helped
make stroke a leading cause of serious, long-term disability in
adults worldwide (
Go etal., 2013; World Stroke Organization,
2016
). More than a third of stroke survivors experience signi-
cant physical disability (
Lloyd-Jones etal., 2010), with walking
dysfunction being among the greatest stroke-related limitations
contributing to disability. While the majority of persons who
suer a stroke regain some level of ambulatory function, their
gait is typically slow, asymmetrical, and metabolically inecient
(
Olney etal., 1986; Roth etal., 1997). Diminished walking ability
is tied to decreased quality of life, increased risk of depression,
and increased risk of serious secondary health conditions
(
Blair et al., 1989; Mutikainen et al., 2011; Ostir et al., 2013).
Restoration of walking function following a stroke is therefore
both a high priority for rehabilitation and an important public
health problem.
Despite recognition of the problem, current clinic-based
neurorehabilitation methods produce only modest improve-
ments in walking function for persons post-stroke (
States etal.,
2009
; Bogey and Hornby, 2014; Winstein etal., 2016). For this
reason, researchers and clinicians have explored a variety of
neurorehabilitation approaches in search of an eective means
to restore post-stroke walking function. ese approaches
include functional electrical stimulation (FES) (
Popovic et al.,
1999
; Kesar etal., 2009, 2010; Sabut etal., 2013; Chung etal.,
2014
; O’Dell etal., 2014; Pilkar etal., 2014; Auchstaetter etal.,
2015
; Chantraine et al., 2016), ankle–foot orthoses (AFOs)
(
Ferreira etal., 2013; Tyson etal., 2013; Kobayashi etal., 2016),
exoskeletons (
Nilsson etal., 2014; Bortole etal., 2015; Buesing
et al., 2015
), partial body weight support (Ng et al., 2008; Lee
etal., 2013
) and split-belt treadmill training systems (Reisman
etal., 2007
; Malone and Bastian, 2014), and robotic gait trainers
(
Pennycott et al., 2012; Mehrholz etal., 2013; Bae et al., 2014;
Hussain, 2014; Dundar et al., 2015). Each of these approaches
has shown varying levels of promise for improving post-stroke
walking function. However, a critical challenge is determin-
ing the treatment prescription – which approach to apply, how
much of the approach to apply, and how the approach should be
applied – that will maximize recovery of walking function for
any particular individual. Furthermore, there is currently no way
to identify whether a small amount of treatment provided by a
combination of approaches might be dramatically more eective
than a large amount of treatment provided by a single approach
(
Belda-Lois et al., 2011). Current treatment design methods
based on trial-and-error and subjective clinical judgment cannot
address these challenges, since they do not provide an objective
means for predicting a patient’s walking function following a
specied treatment or treatment combination.
One possible approach for overcoming this challenge is to
use patient-specic neuromusculoskeletal models to predict
post-treatment walking function for dierent neurorehabilitation
technologies (alone or combined) under consideration. Such
models should account for how the patient interacts with the
treatment approach (
Mooney and Herr, 2016) so that the optimal
prescription can be recommended based on objective predictions
of walking improvement. A number of studies have pursued
such modeling eorts by simulating the eects of FES (
Riener,
1999
; Heilman and Kirsch, 2003; Zhang and Zhu, 2007; Shao and
Buchanan, 2008
; Nekoukar and Erfanian, 2013; Sharma et al.,
2014
; Alibeji etal., 2015), exoskeletons (Fleischer and Hommel,
2008
; Afschri etal., 2014; Farris etal., 2014; Sawicki and Khan,
2015
), orthoses (Zmitrewicz etal., 2007; Crabtree and Higginson,
2009
; Silverman and Neptune, 2012), and strength training
(
Goldberg and Neptune, 2007; Knarr et al., 2014) on lower
extremity function and walking ability in the context of stroke,
spinal cord injury, amputee, and general rehabilitation. Few of
these studies focused on stroke (
Goldberg and Neptune, 2007;
Shao and Buchanan, 2008; Knarr etal., 2014), few used three-
dimensional models (
Fleischer and Hommel, 2008; Afschri
etal., 2014
; Farris etal., 2014; Knarr etal., 2014; Sawicki and Khan,
2015
), few used subject-specic models created by calibrating
critical neuromusculoskeletal model parameters to movement
data collected from an individual (Fleischer and Hommel, 2008;
Shao and Buchanan, 2008; Knarr et al., 2014), and only one
included modeling elements that accounted for subject-specic
neural control capabilities and limitations (Alibeji etal., 2015).
No study to date has predicted a stroke patient’s complete post-
treatment walking motion and speed resulting from application
of a specic neurorehabilitation intervention.
As a rst step toward optimizing patient interaction with
stroke neurorehabilitation technologies, this study describes
the development and evaluation of a subject-specic synergy-
controlled neuromusculoskeletal simulation framework that can
predict three-dimensional walking motions for an individual
post-stroke. e main question we address is whether actuating
a subject-specic neuromusculoskeletal model with muscle syn-
ergy controls (5 per leg) facilitates generation of accurate walking
predictions compared to actuating the model with muscle activa-
tion controls (35 per leg) or joint torque controls (5 per leg). We
hypothesize that synergy controls will work the best since they
combine a low number of control signals with a subject-specic
representation of the coupling between muscle activations within
each leg. We collect gait data from a stroke subject walking at
0.4, 0.5, 0.6, 0.7, and 0.8m/s on an instrumented treadmill and
use data from his self-selected speeds of 0.4–0.6m/s to develop
a subject-specic neuromusculoskeletal model. We incorporate
the subject-specic full-body model into a direct collocation
optimal control framework to predict new walking motions for
the subject. To evaluate the framework and the potential benets
of using synergy controls, we predict how the individual will walk
(including cadence and stride length) at 0.5 and 0.8m/s (condi-
tions for which experimental data are available for comparison)
using joint torque, muscle activation, or muscle synergy controls
and at 1.1m/s (a condition for which no experimental data are
available) using only synergy controls. With future simulation
of dierent neurorehabilitation approaches, our subject-specic
synergy-controlled neuromusculoskeletal simulation framework
may help identify optimal neurorehabilitation prescriptions that
3
Meyer et al. Muscle Synergies Facilitate Walking Predictions
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2016 | Volume 4 | Article 77
maximize recovery of walking function on an individual patient
basis.
METHODS
Experimental Data Collection
To assist with development and evaluation of our subject-
specic synergy-controlled neuromusculoskeletal simulation
framework, we collected experimental walking data from one
high-functioning hemiparetic male individual with chronic
stroke-related walking dysfunction (age 79years, LE Fugl-Meyer
Motor Assessment 32/34 pts, right-sided hemiparesis, height
1.7 m, mass 80.5 kg). All study procedures were approved by
the University of Florida Health Science Center Institutional
Review Board (IRB-01) and the Malcom Randall VA Medical
Center Research and Development Committee and included
approval to study individuals with stroke-related disability. Study
personnel obtained written informed consent prior to participant
enrollment and involvement in study activities. Study procedures
were conducted in accordance with the Declaration of Helsinki.
Motion capture (Vicon Corp., Oxford, UK), ground reaction
(Bertec Corp., Columbus, OH, USA), and electromyography
(EMG) data (Motion Lab Systems, Baton Rouge, LA, USA)
were collected simultaneously while the participant walked on
a split-belt instrumented treadmill (Bertec Corp., Columbus,
OH, USA) at speeds ranging from 0.4 to 0.8m/s in increments of
0.1m/s. 0.8m/s was the fastest speed at which the subject could
walk safely without assistance. is range of speeds included the
participant’s self-selected walking speed of 0.5m/s. More than 50
gait cycles were recorded at each walking speed. A static stand-
ing trial was collected for model scaling purposes. To facilitate
subsequent creation of subject-specic foot–ground contact
models, the participant wore Adidas Samba Classic sneakers,
which have a at sole and neutral midsole with no cushioning,
and we collected additional static trials where we used a marker
wand to trace the outline of each sneaker sole on the ground.
Motion capture data were obtained using a modied Cleveland
Clinic full-body marker set with additional markers added to the
feet (
Reinbolt etal., 2005). Marker motion and ground reaction
data were ltered at a variable cut-o frequency of 7/tfHz, where
tf is the period of the gait cycle being processed, using a fourth-
order zero phase lag Butterworth lter (
Hug, 2011). is variable
cut-o frequency would cause data collected at a normal walking
speed to be ltered at ~6Hz.
Electromyography data were collected from 16 muscles in each
leg and processed using standard methods (Lloyd and Besier,
2003
). A combination of surface and ne-wire EMG electrodes
was used. Surface EMG data were collected for gluteus maximus
and medius, semimembranosus, biceps femoris long head, rectus
femoris, vastus medialis and lateralis, medial gastrocnemius, tibi-
alis anterior, peroneus longus, and soleus. Fine-wire EMG data
were collected for adductor longus, iliopsoas, tibialis posterior,
extensor digitorum longus, and exor digitorum longus. All
EMG data were high-pass ltered at 40Hz (
Lloyd and Besier,
2003), demeaned, rectied, and then low-pass ltered at a vari-
able cut-o frequency 3.5/tfHz. Filtering was performed using a
fourth-order zero phase lag Butterworth lter. EMG data from
each muscle were normalized to the maximum value over all trials
and resampled to 101 time points per gait cycle (heel strike to heel
strike for the less involved le side) while keeping an additional
20 time points before the start of the cycle to permit modeling of
electromechanical delay. In addition, each processed EMG signal
was oset on a cycle-by-cycle basis so that the minimum value
was zero.
Neuromusculoskeletal Model
Development
e subject-specic neuromusculoskeletal model that served as
the foundation for our simulation framework incorporated four
modeling components to account for the unique neurophysio-
logical and musculoskeletal characteristics of the subject: (1) a
subject-specic lower-body kinematic model to simulate the
subject’s skeletal motion, (2) subject-specic foot–ground con-
tact models to simulate how the subject’s feet interact with the
ground, (3) subject-specic EMG-driven muscle moment models
to simulate the subject’s lower extremity joint moments, and (4)
a subject-specic muscle synergy control model to simulate the
subject’s neural control system. Below we describe each of these
four components in further detail. Unless otherwise noted, we
calibrated model parameters in each component using a single
representative walking trial collected at the subject’s self-selected
speed of 0.5m/s.
Subject-Specic Lower-Body Kinematic Model
Our neuromusculoskeletal model creation process started with
a generic full-body musculoskeletal model (
Arnold etal., 2010;
Hamner etal., 2010) developed in OpenSim (Delp etal., 2007).
e original model included 37 degrees of freedom (DOFs) and
44 Hill-type muscle-tendon actuators in each leg. We locked the
wrist and forearm pronation–supination angles to anatomically
reasonable values for walking, leaving the following 31 DOFs:
6 DOF ground-to-pelvis joint, 3 DOF hip joints, 1 DOF knee
joints, 1 DOF ankle joints, 1 DOF subtalar joints, 1 DOF
toejoints connecting rear foot and toe segments, 3 DOF back
joint, 3 DOF shoulder joints, and 1 DOF elbow joints. We also
eliminated nine muscle-tendon actuators without related EMG
data (extensor hallucis longus, exor hallucis longus, gemelli,
gracilis, pectineus, piriformis, quadratus femoris, sartorius,
tensor fascia latae), leaving 35 muscles per leg that actuated
hip exion-extension, hip adduction-abduction, knee exion-
extension, ankle exion-extension, and ankle inversion–ever-
sion on each leg. We then scaled the modied model using the
standing static trial marker data and the OpenSim “Scale Model”
tool, where distances between markers placed over identiable
landmarks were averaged between the two sides to maintain
bilateral symmetry following scaling.
Once the model was scaled, we calibrated joint and marker
positions in the torso, pelvis, and lower-body portion of the
OpenSim model using marker data from a representative walking
trial. e calibration approach was similar to one described pre-
viously (
Reinbolt etal., 2005, 2008) except that it was performed
on the scaled OpenSim model using the OpenSim-MATLAB
Application Programing Interface and included modications
to maintain correct bone geometry positions within the body
4
Meyer et al. Muscle Synergies Facilitate Walking Predictions
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org October 2016 | Volume 4 | Article 77
segments (Charlton et al., 2004). To facilitate the calibration
process, we created marker plates on the torso, pelvis, thighs,
and shanks to which all markers on the respective OpenSim
body segments were attached. To perform the actual calibration,
we used non-linear least squares optimization (lsqnonlin) in
MATLAB to adjust joint (knee, ankle, and subtalar in both legs)
and marker plate (torso, pelvis, thighs, and shanks) positions
and orientations in their respective body segments such that
marker errors from repeated OpenSim inverse kinematic analy-
ses were minimized. e optimization cost function included
penalty terms that prevented large changes in joint and marker
plate positions and orientations that would produce only small
improvements in marker tracking. Modication of the two hip
joint center locations was achieved by modifying the position
and orientation of the rigid marker plate on the pelvis. For joint
centers and orientations, symmetry between le and right sides
of the body was maintained during the kinematic calibration
process. Markers on the feet were not adjusted since their loca-
tions were well dened. e position and orientation of the toe
axis in each foot and of the back, shoulder, and elbow joints was
maintained from the scaled OpenSim model.
Subject-Specic Foot–Ground Contact Models
Following kinematic calibration, we created a subject-specic
foot–ground contact model for each foot of the OpenSim model
using recently developed methods (
Jackson et al., 2016). e
elastic foundation contact models were developed in MATLAB
and used a grid of contact elements that spanned the rear foot
and toes segments of each foot. To create the element grid, we
started with the shoe outlines obtained from the static trial
marker data and used principal component analysis to identify
the principal axes of each foot (rear foot and toes segments
together). Using these axes, we constructed an 11 (anterior-
posterior) × 8 (medial-lateral) grid of rectangular contact
elements for the le foot, where the edges of the grid extended
2.5mm beyond the edge of the shoe outline in both directions.
Forty-seven elements whose centers fell within the shoe outline
were retained in the contact model, while 41 elements whose
centers fell outside the shoe outline were removed. Given the
locations of the MATLAB contact element centers relative to the
foot markers from the static le shoe outline trial, we calculated
the locations of the element centers on the OpenSim rear foot
and toes segments. We then projected the le toes axis of the
OpenSim model onto the contact element grid. Elements whose
centers were posterior to the axis were assigned to the rear foot
segment, while elements whose centers were anterior to the axis
were assigned to the toes segment. e complete MATLAB/
OpenSim contact element grid for the le foot was mirrored to
the right foot by aligning the principal axes of the mirrored grid
with those of the right foot.
Each contact element in the foot–ground contact models
generated normal force using a linear spring with non-linear
damping and shear force using a continuous stick-slip friction
model. For any contact element i, the required time-varying
inputs for contact force calculations performed in MATLAB were
the penetration into the oor y
i
, the normal penetration rate
y
i
,
and the shear slip rate
v
s
i
of the element center in the laboratory
coordinate system as calculated by OpenSim. e normal contact
force F
i
for element i was calculated as (
Hunt and Crossley, 1975)
Fkyy cy
iiii
=−+()
()
0
1
(1)
where k
i
is the spring stiness unique to each spring, y
0
is the
spring resting length common to all springs on the same foot
and essentially adjusts the height of the oor, and c (=1e–2) is a
non-linear damping coecient common to all springs. e linear
spring also generates a small amount of force when the foot is o
the oor, and this negligible force transitions in a smooth and
continuous way to the large force produced when the spring is in
contact with the ground (
Anderson and Pandy, 2001; Ackermann
and van den Bogert, 2010
). e non-linear damping ensures
that the normal contact force does not exhibit damping-related
discontinuities when a spring enters or leaves contact. e shear
contact force f
i
for element i was calculated using a simple con-
tinuous and dierentiable friction model (
Ackermann and van
den Bogert, 2010
)
fF
v
v
ii
s
l
i
=
µtanh
(2)
where μ [=1 (
Ackermann and van den Bogert, 2010)] is a dynamic
friction coecient common to all springs and v
l
(=5cm/s) is a
latching speed common to all springs that denes the edge of a
linear transition region between zero slip rate and the start of
dynamic friction. Shear contact force f
i
was applied to the element
center in the direction opposite to the slip velocity vector. e
direction calculation included a small constant value of 1e−4 in
the denominator to avoid division by a small number when the
slip speed was near zero. Once F
i
and f
i
were calculated for each
contact element, the net contact force and torque due to all contact
elements in the rear foot segment were calculated with respect to
the rear foot origin, and similarly for the toes segment using the
toes origin (
Kane and Levinson, 1985). ese net contact forces
and torques were then applied to their corresponding segments
in the OpenSim model. is approach allowed rear foot and toes
contributions to total ground reaction force to be predicted by
the model.
We calibrated the spring stiness k
i
of each contact element
in both feet and the spring resting length y
0
for all contact ele-
ments in each foot using marker and ground reaction data from
a representative walking trial. We made two assumptions about
the spring stiness distribution across the bottom of the shoe to
simplify the calibration process. First, we assumed that the mir-
rored stiness distribution was the same for both feet. Second,
we assumed that the stiness distribution across the entire shoe
bottom could be approximated by a three-dimensional parabolic
surface, which possesses only six unknown coecients rather
than 47 unknown independent spring stiness values and pre-
vents neighboring springs from having dramatically dierent
stinesses. To calibrate these six coecients and the two resting
lengths, we formulated a direct collocation optimal control
problem that tracked experimental marker, ground reaction, and
inverse dynamic joint torque data for the entire body with higher
weight placed on matching marker position and ground reaction
data for the two feet. Tracked ground reaction quantities for each