HAL Id: lirmm-00798599
https://hal-lirmm.ccsd.cnrs.fr/lirmm-00798599
Submitted on 8 Mar 2013
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
A closed loop musculoskeletal model of postural
coordination dynamics
Vincent Bonnet, Philippe Fraisse, Nacim Ramdani, Julien Lagarde, Soane
Ramdani, Benoit G. Bardy
To cite this version:
Vincent Bonnet, Philippe Fraisse, Nacim Ramdani, Julien Lagarde, Soane Ramdani, et al.. A closed
loop musculoskeletal model of postural coordination dynamics. CDC/CCC’09: 48th IEEE Conference
on Decision and Control, held jointly with the 28th Chinese Control Conference, Dec 2009, Shanghai,
China. pp.6207-6212, �10.1109/CDC.2009.5400945�. �lirmm-00798599�
A closed loop musculoskeletal model of postural coordination dynamics
Vincent Bonnet, Philippe Fraisse, Nacim Ramdani, Julien Lagarde, Sofiane Ramdani, Benoit G. Bardy
Abstract— A closed-loop model with actuator dynamics and
sensory feedback has been developed to capture the complex
postural behaviors observed in a human head tracking task. In
motor-control litterature, spindle feedback gains are scaled by
the central nervous system to adapt muscle stiffness depending
on the postural task. We propose to identify spindle reflex
equivalent feedback gains for several target’s frequency values.
Comparison with experimental results on human shows the
relevance of this modeling, since our musculoskeletal model
is able to exhibit reasonably well the behavioral invariants
observed in human postures.
Index Terms— Human Postural Coordination, Dynamical
Postural controller, Redundant Systems, Biological system mod-
eling, Parameter estimation.
I. INTRODUCTION
The coordination and control of human posture (the spatio-
temporal organization of body joints) and balance have
been modelled in several ways using many experimental
paradigms, in tasks involving postural responses to external
perturbations or to intentional tracking tasks.
During a visual tracking task, Bardy et al [1] analyzed,
in the framework of coordination dynamics [2], the joint
coordination in the sagittal plane. They proposed the use
of a collective variable to describe in a simple way the
complex dynamical biological-segmental-articular-muscular-
couplings. They had standing participants moving back and
forth in the sagittal plane in order to track the displacement
of a virtual target. This simple task allowed the observation
of several self-organized properties of the postural system,
such as phase transition, multistability, critical fluctuations,
hysteresis, and critical slowing down. The collective variable
able to capture both fully and in a very compact way these
dynamical properties, is the relative phase, i.e. the phase
difference between the ankle and the hip. Two coordination
modes were observed between ankle and hip depending on
the target’s frequency: An in-phase mode for low frequencies
and an anti-phase mode for high frequencies.
Martin et al. [3] used a constrained optimization process
to analyze Bardy et al.’s [1] results, and showed that the
location of the center of pressure (CoP) can drive the
selection of the coordination mode. In a previous work [4],
we studied in more details Bardy et al.’s [1] paradigm with
V. Bonnet, P. Fraisse and N. Ramdani are with the LIRMM UMR 5506
CNRS, Univ. Montpellier 2, 161 rue Ada, Montpellier, 34392 France.
N. Ramdani is with the INRIA Sophia-Antipolis, Nice FR-06902, France.
J. Lagarde, S. Ramdani and B. G. Bardy are with the EDM EA 2991,
Univ. Montpellier 1, EDM, 700 av. du pic Saint Loup, Montpellier, 34090
France.
B. G. Bardy is with Institut Universitaire de France, 103 Bd St Michel,
Paris, 75005 France
a method similar to the one used in [3]. Then we imple-
mented the obtained coordination modes onto the HOAP-3
and HRP2 humanoid robots. We showed that the in-phase
mode corresponds to the minimum energy mode for low
frequencies, and that the anti-phase mode is the only one
able to maintain balance for high frequencies.
However, the approach described above considers only
steady state behaviors and thus is not capable of capturing
the transient dynamics observed during human postural be-
haviors such as the hysteresis phenomenon for instance. Non-
linear coupled oscillators are classi cally used to model these
human dynamical coupling phenomena [5]. However, these
oscillators involve several unknown parameters which have
to be tuned and whose connection with the actual system is
difficult to delineate.
In tasks involving postural responses to external pertur-
bations, Nasher and Mc Collum [6] observed two postural
strategies at the muscular activation level. One of two
postural strategies are adopted according to perturbation
magnitude: either a large muscular activation level at the
ankles, or a coordinated activation of the hips and the
ankles. Hemami [7] and colleagues have found that rea-
sonable predictions behaviour can be made using linearized
dynamics to model small perturbations postural response.
Based on these postural strategies and assumptions, Kuo [8]
developed a double inverted pendulum (DIP) model with
an optimal control law weighting excursion of the center
of mass and deviations from the upright position, and with
feasible acceleration set (FAS) which is connected to the
physiological limitations and equilibrium constraint. Using
the FAS framework, Park et al. [9] optimized joint feedback
gains to fit human data for different kind of perturbations
and showed that the trajectories of joint angles and joint
torques scale with perturbation magnitude in agreement with
the postural strategy observations. In the same way, Jo et
al. [10] proposed a cerebrocerebello-spinomuscular model of
human posture with time delays and proportional-integrative-
derivative feedbacks ; their controller uses different set of
cerebellum gains depending on perturbation magnitude.
These modeling approaches have inspired our work in
the postural coordination framework. In this paper, we re-
port the non-linear closed loop model we are developing
in order to capture and predict postural sway movements
during head tracking tasks. Our model is composed of a
double inverted pendulum as biomechanical model, muscles
models at each joint and a classical proportional-derivative
controller i n the operational space. In this sequel, we will
report human experimentations first, in order to describe the
postural coordination concept, and then we will describe
the modelling and the identification of model parameters in
subsequent sections.
II. HUMAN EXPERIMENTATION
The aim of the experiment is to provide a database of
human behavior in order to identify spindle feedback gains
in section IV.
A. Methods
Following previous studies [11], [12], the experiment
consists in tracking a moving target with the head while
standing. Participants stood on a force platform in front of a
physical target moved by a linear motor in antero-posterior
direction, with the knees locked and the soles constantly in
contact with the ground (Fig. 1-2-3).
Fig. 1. Experimental device. Physical target moved by a linear motor, force
plate and motion capture d evice.
(a) t=0 sec (b) t=1 sec (c) t=2.5 sec
Fig. 2. Human typical experiments at 0.2Hz. Coordinative in-phase small
displacement of the ankle and the hip.
(a) t=0 sec (b) t=0.5 sec (c) t=1.3 sec
Fig. 3. Human typical experiments at 0.6Hz. Coordinative anti-phase
displacement of the ankle and the hip. The hip amplitude is larger than
the ankle one.
The experiment was performed on 11 healthy male sub-
jects, with mean age 25, mean weight 75kg and mean
size 1.79m. Target motion was sinusoidal with 10cm as
amplitude, the frequency increases from 0.1Hz to 0.65Hz
by 0.05Hz steps and during 10 periods. To capture the joint
positions, a motion capture s ystem (VICON NEXUS) was
used, with 8 cameras (MX13) tracking 15 makers on the
right side of the subject.
0.1 0.2 0.3 0.4 0.5 0.6
0
50
100
150
200
frequency (Hz)
a. (°)
Relative phase Ankle/Hip with Hilbert Transform
0.1 0.2 0.3 0.4 0.5 0.6
−0.2
−0.1
0
0.1
0.2
frequency (Hz)
b. (rad)
Joint positions mean amplitudes per step frequency
Ankle
Hip
0.1 0.2 0.3 0.4 0.5 0.6
20
40
60
80
100
c. (N.m)
frequency (Hz)
Torque mean amplitudes per step frequency
Ankle
Hip
Fig. 4. Typical human experimental results. (a) Ankle/hip relative phase,
showing a transition frequency around 0.4Hz (b) Peak to peak joint
positions. (c) Estimation of joint torque amplitudes.
B. Experimental results
Fig.4 shows typical results for a representative subject
(weight 75kg, size 1.80m). On Fig. 4a, the mean values of the
relative phase (Hilbert-transformed) between ankle and hip
positions are represented as a function of the frequency step.
The depicted error bars correspond to the standard deviations
during the 10 oscillations achieved at each frequency step.
A transition is observed from in-phase to anti-phase mode
around 0.4Hz. Joint positions are presented on Fig. 4b by
minima and maxima values. Each point is the mean value of
the maximum (or minimum) joint position reached during the
10 oscillation periods performed at each frequency step. For
the in-phase mode, i.e., at low frequencies, the joint positions
amplitude difference are small, with individual differences
in terms of joint amplitude. At the transition frequency,
the ankle amplitude become very small (Fig. 10), and the
relative phase between ankle and hip is difficult to estimate.
This strong reduction of the ankle amplitude is typical of
human phase transition [11], [12], [13]. The hip amplitude
is larger than the ankle amplitude for the anti-phase mode
as mentioned in [11], [12], [13]. Fig. 4(c) depicts mean
values for torque amplitude estimation at each fr equency
step. Torque values were estimated by using the inverse
dynamical model of the DIP. They indicate a larger ankle
torque amplitude for in-phase mode and a larger hip torque
amplitude for anti-phase in agreement with the ankle and hip
strategy reported in [6] and by Runge et al. [14].
These observations hold for all participants and are in
accordance with [11], [1], [13], even though the actual
transition frequency and joint amplitudes depend on the
specific subject body type.
III. MODELING POSTURAL BEHAVIOR
The model of postural coordination consists of 4 major
components: the rigid body mechanics, the muscle dynamics,
the reflex feedback loop and the closed-loop controller.
A. Biomechanical model
Barin [15] shows the relevance of an inverted pendulum
structure in the case of a human sagittal plane task. In
addition, the Bardy’s paradigm focused on the hip and ankle
joints, so a DIP in the sagittal plane is used in this paper as
a biomechanical model (Fig. 5).
Fig. 5. Double inverted pendulum used to model postural coordination
and a musculoskeletal model. The moment arms of muscular forces are the
pulleys radius.
Balance is described by the position of the CoP within the
BoS, which can be expressed as a function of the dynamic
parameters (eq.1).
X
CoP
= (Γ
1
− F
gx
d + m
0
k
0
g)/F
gy
(1)
where F
gx
is the horizontal ground reaction force, F
gy
the
vertical one, Γ
1
the ankle torque, and m
0
and k
0
foot
parameters. Euler’s equations were used for the calculation
of the ground reaction forces as proposed by Cahouet et
al. [16]. The muscle length evolutions, function of the
joint positions, are computed by using a double pulley
model at each joint (see Fig. 5). Muscle insertion points and
parameters are adapted from [17].
B. Joint actuator dynamics
Musculotendon contraction dynamics is often modeled
with Hill type models. For postural applications, Hill type
models are often linearized [18], but in our application, we
will keep the non-linear dynamic phenomenon. Furthermore
many musculoskeletal models as [19] consider the tendon
completely stiff. To provide a realistic model, the compliant
tendon in the model needs to be addressed because of the
long ankle tendons. Then a 3 element Hill-type model (Fig.
6) is used, based on [20]. The muscle model is composed
of an active contractile element (CE) a passive parallel ele-
ment (PE), which represents the intrinsic viscoelastic muscle
properties, and a serial li near tendon (SE). Normalized units
are used in this model to produce the output force F
Mt
,
which results of the interaction between the force produced
by muscle and the tendon stiffness, i.e. F
Mt
= F
ce
+F
pe
. The
contractile element force is the product of the force-length
F
l
and the force-velocity F
v
relationships and the activation
Fig. 6. Hill type musculotendon model
act, as follows:
F
ce
= act.F
l
.F
v
. (2)
Note that act
min
< act < 1 by structure of the model. The
retained F
l
relation is F
l
= exp(−((L
ce
− L
ce
0
)/(L
ce
sh
))
2
),
where L
ce
0
is the optimal muscle fiber length and L
ce
sh
a
constant shape parameter. The Hill force-velocity relation can
be expressed as follows:
F
v
=
0 V
ce
< −V
ce
max
b+aV
ce
b−V
ce
−V
ce
max
< V
ce
≤ 0
b+V
ce
(F
v
max
+(F
v
max
−1)a)
b+V
ce
V
ce
> 0
where a, b et F
v
max
are constant parameters, V
ce
is the
contraction velocity of muscle fiber. The maximum veloc-
ity contraction of the muscle fiber is commonly taken as
V
ce
max
= 10.L
ce
0
. The force generated by the PE is given
by : F
pe
= (e
−
pe
sh
pe
xm
(l
ce
−1)
−1
)/(e
pe
sh
− 1), where pe
sh
and pe
xm
are constant parameters. Regarding the tendon
stiffness, it can be approximated by a linear spring as shown
in Zajac [20].
C. Muscle redundancy and torque feedback
A human joint is actuated by a pair of agonist-antagonist
muscle group. Therefore muscle coactivation needs to be
addressed. Previous human experimentation on s inusoidal
standing perturbations [21] have shown a very small coac-
tivation between the agonist and antagonist muscle group.
Therefore, the coactivation will be neglected in our modeling.
The total joint net torque in our modeling is the following:
Γ = r
a
F
Mt
a
− r
p
F
Mt
p
(3)
where r
a
and r
p
are the constant lever arm (since a double
pulley model is used), F
Mt
a
and F
Mt
p
are the agonist and
antagonist force developed by the muscles.
In human, a muscle force feedback is given by the Golgi
Tendon Organs (GTO), since the GTO is in serial with the
tendon, the gain K
GT O
is known to be constant. Since there
is no coactivation in our model and since the lever arms are
constants, the GTO feedback is the product of each muscle
force with muscle lever arm.
Finally, in order to address muscle redundancy at the joint
level (for a muscle group), we use the control scheme based
on [22] depicted on Fig. 7. A torque feedback is given by
Golgi Tendon Organs. The pulling muscle selection is based
on the sign of the torque error. a
a
and a
p
are the agonist and
antagonist activation signals which appear in Eq. (2); they are
saturated to 1. (Q
T
)
+
is the transformation matrix fr om the
joint-space to muscle space and Γ
d
the desired joint torque
which is an image of an higher level supra-spinal control
variable U
ss
since K
H
is a homogeneity gain fixed to 1. The
muscle model activation must be positive and since (Q
T
)
+
is constant, the selection of the pulling muscle is based on
the sign of the torque error.
Fig. 7. Block diagram of the muscle redundancy management.
D. Spindle reflex feedback
In this study, the restraining reflex joint torques are con-
sidered in the model in order to limit t he joint motion and
to improve the fit with human joint positions. In human,
the contraction velocity and the length of the musculotendon
system are given respectively by II and Ia spindles sensors.
It’s well known that spindle feedback reflexes (SFR) are
important in postural muscles, especially in terms of stiffness
control and dis turbance rejection. Many spindle models exist,
at different level of physiology description, e.g. [23], but
generally they are modelled with a pure gain and a small
time delay [17] (35ms). Because SFR time delay is small
and since no other time delays are introduced in our model,
there will be neglected. Future works will introduce a non-
linear spindle model with time delay, in order to have a better
physiological description.
Spindle sensors are located in parallel with muscle fibers
(CE element), so their length change with L
ce
length. Then
CNS adapt the SFR gains, in order keep a good measure sen-
sitivity. In addition, the point-equilibrium-theory in motor-
control literature, proposed by Feldman [24], argues that the
CNS modulates the SFR gains to control the muscle stiffness
and hence the joint position. This theory is the target of a
vivid debate but there is converging evidence. Therefore, we
will tune the SFR gains in section IV-A.
E. Closed-loop modeling
The CNS needs to manage redundant sets of actuators and
sensors to perform the tracking task. Thus postural control
can be assumed as an optimal control problem for the CNS.
The Bardy et al.’s paradigm is composed at least of two main
tasks: a standing task and a tracking task.
Firstly, a proportional-derivative (P D) controller, in the
operational space, is chosen as neural controller model to
perform the tracking task. Since Masani et al. [25] argues
that P D controller is a good approximation of the control
strategy in human standing, specially for the robustness on
time delay (which will be the next step of our work). The
actual tracking task sets only the horizontal head position
(X axis in Fig. 5), so the actuated system is redundant
with respect to the task. Therefore we use a jacobian
pseudoinverse matrix to address the joint control vector. It
is well attested that the use of the pseudoinverse matrix
in kinematics redundant problems minimizes the norm of
the velocity vector ||
˙
θ||
2
at a given time. By analogy with
inverse kinematics, the pseudoinverse matrix used in the
control scheme depicted on Fig. 8 minimizes the norm of
the supraspinal vector ||U
ss
(t)||
2
. Then our model behaves
like an optimal controller.
Secondly, standing must be guaranteed, that is the DIP
must stay relatively close to the upright position and the
CoP inside the base of support. This equilibrium constraint
is managed by the SFR loops which constrain muscle lengths
(i.e. joint positions) close to rest positions. In other words, it
means that the vertical position of the head (Y axis in Fig.
5 ) is virtually bounded. Consequently, the proposed closed-
loop model scheme minimises an energetic criterion while
minimising the joint deviations from rest position.
Therefore our non-linear control scheme is composed
of a double inverted pendulum as biomechanical model, a
muscle group (see section III-C), a classical P D controller
in operational space, a torque feedback representing the GTO
and joint positions and velocities feedback representing the
SFR, as shown on Fig. 8.
On Fig. 8, DKM is the direct kinematics model, J
+
is
the pseudoinverse matrix and K
Ia
and K
II
respectively the
position and velocity feedback gains at spindles level.
IV. IDENTIFICATION RESULTS
A. Identification of equivalent spindle feedback gains
Based on motor-control literature and as argued in section
III-D, we propose to identify for each frequency step the
SFR gains, i.e. KII
1
, KIa
1
, KII
2
and KIa
2
, where index
1 stands for ankle and index 2 for hip. All other model
parameters are taken constant and known prior to identi-
fication, chosen according to physiological considerations.
They correspond to the same typical subject described in
section II. In addition, the operational space gains are taken
as P = 800 and D = 1000, in order to let the closed-loop
model follow the desired head position with a good accuracy
for all frequencies. As regarding physiological meanings, this
tuning process is equivalent to adjusting as close as possible
the dynamical reflex behaviors of each movement at a given
frequency, to the actual ones.
In this study, SFR gains identification is addressed in the
stochastic framework, where the maximum-likelihood ap-
proach makes it possible to derive a criterion to be optimized
to estimate these gains and an asymptotic uncertainty domain
associated with the est imated gains [26]. Under the usual
Gaussian assumptions for the probability density function of
the data and measurement errors, the maximum-likelihood
estimator boils down to the least-square estimator, which
minimizes the quadratic norm of the following output error:
J =
2
X
i=1
10T
X
t=0
(θ
i
hum
− θ
i
sim
)
2
(4)
