# Deep Feature Mining via Attention-based BiLSTM-GCN for Human Motor Imagery
Recognition

Yimin Hou, Shuyue Jia, , Xiangmin Lun, Shu Zhang, Tao Chen, Fang Wang, Jinglei
Lv  This work was supported by the National Natural Science Foundation of
China under Grant 31772059. Yimin Hou, Xiangmin Lun, Shuyue Jia, Tao Chen, and
Fang Wang are with the School of Automation Engineering, Northeast Electric
Power University, Jilin, China. Xiangmin Lun is also with the College of
Mechanical and Electric Engineering, Changchun University of Science and
Technology, Changchun, China. Shu Zhang is with the School of Computer
Science, Northwestern Poly-technical University, Xi’an, China. Jinglei Lv is
with School of Biomedical Engineering & Brain and Mind Center, University of
Sydney, Sydney, NSW, Australia.

###### Abstract

Recognition accuracy and response time are both critically essential ahead of
building practical electroencephalography (EEG) based brain-computer interface
(BCI). However, recent approaches have either compromised in the
classification accuracy or responding time. This paper presents a novel deep
learning approach designed towards both remarkably accurate and responsive
motor imagery (MI) recognition based on scalp EEG. Bidirectional Long Short-
term Memory (BiLSTM) with the Attention mechanism is employednvolutional
neural network (GCN) promotes the decoding performance by cooperating with the
topological structure of features, which are estimated from the overall data.
Particularly, this method is trained and tested on the short EEG recording
with only 0.4 second length, and the result has shown effective and efficient
prediction based on individual and group-wise training, with 98.81% and 94.64%
accuracy, respectively, which outperformed all the state-of-the-art studies.
The introduced deep feature mining approach can precisely recognize human
motion intents from raw and almost instant EEG signals, which paves the road
to translate the EEG based MI recognition to practical BCI systems.

###### Index Terms:

Brain-computer Interface (BCI), Electroencephalography (EEG), Motor Imagery
(MI), Bidirectional Long Short-term Memory (BiLSTM), Graph Convolutional
Neural Network (GCN), Attention Mechanism

##  I Introduction

Recently, brain-computer interface (BCI) plays a promising role in assisting
and rehabilitating patients from paralysis, epilepsy, and brain injuries via
interpreting neural activities to control the peripherals[1, 2]. Among the
non-invasive brain activity acquisition systems, electroencephalography
(EEG)-based BCI gains extensive attention recently given its higher temporal
resolution and portability. Hence, it has been popularly employed to assist
the recovery of patients from motor impairments, e.g., amyotrophic lateral
sclerosis (ALS), spinal cord injury (SCI), or stroke survivors[3, 4].
Specifically, researchers have focused on the recognition of motor imagery
(MI) based on EEG and translating the brain activities into specific motor
intentions. In such a way, users can further manipulate external devices or
exchange information with the surroundings[4]. Although researchers have
developed several MI-based prototype applications, there is still space of
improvement before the practical clinical translation could be promoted[2, 5].
De facto, to achieve effective and efficient control via only MI, both precise
EEG decoding and quick response are eagerly expected. However, few existing
works of literature are competent in both perspectives. In this study, we
explore the possibility of a deep learning framework to tackle the challenge.

Figure 1: The schematical overview consisted of 64-channel raw EEG signals
acquisition, the BiLSTM with Attention model for feature extraction, and the
GCN model for classification.

###  I-A Related Work

Lately, Deep Learning (DL) attracted increasing attention in many disciplines
because of its promising performance in classification tasks[6]. A growing
number of works have shown that DL will play a pivotal role in the precise
decoding of brain activities[2]. Especially, recent works have been carried
out on EEG motion intention detection. A primary current focus is to implement
the DL-based approach to decode EEG MI tasks, which have attained promising
results[7]. Due to the high temporal resolution of EEG signals, methods
related to the recurrent neural network (RNN)[8], which can analyze time-
series data, were extensively applied to filter and classify EEG sequences,
i.e., time points[9, 10, 11, 12, 13]. In reference[9], a novel RNN framework
with spatial and temporal filtering was put forward to classify EEG signals
for emotion recognition, and achieved 95.4% accuracy for three classes with a
9 seconds’ segment as a sample. Wang _et al._ and Luo _et al._ performed Long
Short-term Memory (LSTM)[14] to handle time slices’ signals, and achieved
77.30% and 82.75% accuracy, respectively[10, 11]. Reference [13] presented
Attention-based RNN for EEG-based person identification, which attained 99.89%
accuracy for 8 participants at the subject level with 4 seconds’ signals as a
sample. However, it can be found that in these studies, signals over
experimental duration were recognized as samples, which resulted in a slow
responsive prediction.

Apart from RNN, Convolutional Neural Network (CNN)[15, 16] has been performed
to decode EEG signals as well[17, 18]. Hou _et al._ proposed ESI and CNN, and
achieved competitive results, i.e., 94.50% and 96.00% accuracy at the group
and subject levels, respectively, for four-class classification. What is more,
by combining CNN with graph theory, the Graph Convolutional Neural Network
(GCN)[19, 20, 21, 22, 23] approach was presented lately, taking consideration
of the functional topological relationship of EEG electrodes[24, 25, 26, 27].
In references[24] and [25], a GCN with broad learning approach was proposed,
and attained 93.66% and 94.24% accuracy, separately, for EEG emotion
recognition. Song _et al._ and Wang _et al._ introduced dynamical GCN (90.40%
accuracy) and phase-locking value-based GCN (84.35% accuracy) to recognize
different emotions[26, 27]. Highly accurate prediction has been accomplished
via the GCN model. But few researchers have investigated the approach in the
area of EEG MI decoding.

Figure 2: Presented BiLSTM with the Attention mechanism for feature
extraction.

###  I-B Contribution of This Paper

Towards accurate and fast MI recognition, an Attention-based BiLSTM-GCN was
introduced to mine the effective features from raw EEG signals. The main
contributions were summarized as follows.
(i) To our best knowledge, this work was the first that combined BiLSTM with
the GCN to decode EEG tasks.
(ii) The Attention-based BiLSTM successfully derived relevant features from
raw EEG signals. Followed by the GCN model, it enhanced the decoding
performance by considering the internal topological structure of features.
(iii) The proposed feature mining approach managed to decode EEG MI signals
with stably reproducible results yielding remarkable robustness and
adaptability, even countering considerable inter-trial and inter-subject
variability.

###  I-C Organization of This Paper

The rest of this paper was organized in the following. The preliminary
knowledge of the BiLSTM, Attention mechanism and GCN was introduced in Section
II, which was the foundation of the presented approach. In Section III,
experimental details and numerical results were presented, followed by the
conclusion in Section IV.

##  II Methodology

###  II-A Pipeline Overview

The framework of the proposed method was demonstrated in Figure 1.
(i) 64-channel raw EEG signals were acquired via the BCI2000, and then the 4
seconds’ (experimental duration) signals were sliced into 0.4-second segments
over time, where the dimension of each segment was 64 channels $\times$ 64
time steps.
(ii) The Attention-based BiLSTM was put forward to filter 64-channel (spatial
information) and 0.4-second (temporal information) raw EEG data and derived
features from the fully-connected neurons.
(iii) The Pearson, Adjacency, and Laplacian Matrices of overall features were
introduced sequentially to represent the topological structure of features,
i.e., as a graph. Followed by the features and its corresponding graph
representation as the input, the GCN model was performed to classify four-
class MI tasks.

###  II-B BiLSTM with Attention

####  II-B1 BiLSTM Model

RNN-based approaches have been extensively applied to analyze EEG time-series
signals. An RNN cell, though alike a feedforward neural network, has
connections pointing backward, which sends its output back to itself. The
learned features of an RNN cell at time step $t$ are influenced by not only
input signals $\mathbf{x}_{(t)}$, but also the output (state) at time step
$t$-1. This design mechanism dictates that RNN-based methods can handle
sequential data, e.g., time points signals, by unrolling the network through
time. The LSTM and Gated Recurrent Unit (GRU)[28] are the most popular
variants of the RNN-based approaches. In Section II-D, the paper compared the
performance of the welcomed models experimentally, and the BiLSTM with
Attention displayed in Figure 2 outperformed others due to better detection of
the long-term dependencies of raw EEG signals.

$\mathbf{i}_{(t)}=\sigma\left(\mathbf{W}_{xi}^{T}\cdot\mathbf{x}_{(t)}+\mathbf{W}_{hi}^{T}\cdot\mathbf{h}_{(t-1)}+\mathbf{b}_{i}\right)$
(1)
$\mathbf{f}_{(t)}=\sigma\left(\mathbf{W}_{xf}^{T}\cdot\mathbf{x}_{(t)}+\mathbf{W}_{hf}^{T}\cdot\mathbf{h}_{(t-1)}+\mathbf{b}_{f}\right)$
(2)
$\mathbf{o}_{(t)}=\sigma\left(\mathbf{W}_{xo}^{T}\cdot\mathbf{x}_{(t)}+\mathbf{W}_{ho}^{T}\cdot\mathbf{h}_{(t-1)}+\mathbf{b}_{o}\right)$
(3)
$\mathbf{g}_{(t)}=\tanh\left(\mathbf{W}_{xg}^{T}\cdot\mathbf{x}_{(t)}+\mathbf{W}_{hg}^{T}\cdot\mathbf{h}_{(t-1)}+\mathbf{b}_{g}\right)$
(4)
$\mathbf{c}_{(t)}=\mathbf{f}_{(t)}\otimes\mathbf{c}_{(t-1)}+\mathbf{i}_{(t)}\otimes\mathbf{g}_{(t)}$
(5)
$\mathbf{y}_{(t)}=\mathbf{h}_{(t)}=\mathbf{o}_{(t)}\otimes\tanh\left(\mathbf{c}_{(t)}\right)$
(6)

As illustrated by Figure 2, three kinds of gates manipulate and control the
memories of EEG signals, namely, the input gate, forget gate, and output gate.
Demonstrated by the $\mathbf{i}_{(t)}$, the input gate partially stores the
information of $\mathbf{x}_{(t)}$, and controls which part of it should be
added to the long-term state $\mathbf{c}_{(t)}$. The forget gate that is
controlled by the $\mathbf{f}_{(t)}$ decides which piece of the
$\mathbf{c}_{(t)}$ should be overlooked. And the output gate, controlled by
$\mathbf{o}_{(t)}$, allows which part of the info from $\mathbf{c}_{(t)}$
should output, denoted as $\mathbf{y}_{(t)}$, as known as the short-term state
$\mathbf{h}_{(t)}$. Manipulated by the above gates, two kinds of states are
stored. The long-term state $\mathbf{c}_{(t)}$ travels through the cell from
left to right, dropping some memories at the forget gate and adding something
new from the input gate. After that, the info passes through a non-linear
activation function, tanh activation function usually, and then it is filtered
by the output gate. In such a way, the short-term state $\mathbf{h}_{(t)}$ is
produced.

Equation 1$-$Equation 6 describe the procedure of an LSTM cell, where
$\mathbf{W}$ and $\mathbf{b}$ are the weights and biases for different layers
to store the memory and learn a generalized model, and $\sigma$ is a non-
linear activation function, i.e., sigmoid function used in the experiments.
For bidirectional LSTM, BiLSTM for short, the signals $\mathbf{x}_{(t)}$ are
inputted from left to right for the forward LSTM cell. What is more, they are
reversed and inputted into another LSTM cell, the backward LSTM. Thus it
leaves us two output vectors, which store much more comprehensive information
than a single LSTM cell. Then, they are concatenated as the final output of
the cell.

####  II-B2 Attention Mechanism

The Attention mechanism, imitated from the human visual, has a vital part to
play in the field of Computer Vision (CV), Natural Language Processing (NLP),
and Automatic Speech Recognition (ASR)[29, 30, 31, 32]. Not all the signals
contribute equally towards the classification. Hence, an Attention mechanism
$\mathbf{s}_{(t)}$ is jointly trained as a weighted sum of the output of the
BiLSTM with Attention based on the weights.

$\mathbf{u}_{(t)}=\tanh\left(\mathbf{W}_{w}\mathbf{y}_{(t)}+\mathbf{b}_{w}\right)$
(7)
$\mathbf{\alpha}_{(t)}=\frac{\exp\left(\mathbf{u}_{(t)}^{\top}\mathbf{u}_{w}\right)}{\sum_{t}\exp\left(\mathbf{u}_{(t)}^{\top}\mathbf{u}_{w}\right)}$
(8) $\mathbf{s}_{(t)}=\sum_{t}\mathbf{\alpha}_{(t)}\mathbf{y}_{(t)}$ (9)

$\mathbf{u}_{(t)}$ is a Fully-connected (FC) layer for learning features of
$\mathbf{y}_{(t)}$, followed by a Softmax layer which outputs a probability
distribution $\mathbf{\alpha}_{(t)}$. $\mathbf{W}_{w}$, $\mathbf{u}_{w}$, and
$\mathbf{b}_{w}$ denote trainable weights and bias, respectively. It selects
and extracts the most significant temporal and spatial information from
$\mathbf{y}_{(t)}$ by multiplying $\mathbf{\alpha}_{(t)}$ with regard to the
contribution to the decoding tasks.

(a) GAA w.r.t. RNN-based Models

(b) GAA w.r.t. the BiLSTM Cell Size

(c) GAA w.r.t. the Attention Size of the BiLSTM

(d) GAA w.r.t. Num. of Extracted Features

(e) Loss w.r.t. RNN-based Models

(f) Loss w.r.t. the BiLSTM Cell Size

(g) Loss w.r.t. the Attention Size of the BiLSTM

(h) Loss w.r.t. Num. of Extracted Features

Figure 3: Models and Hyperparameters Comparison w.r.t. the RNN-based Methods
for Feature Extraction

###  II-C Graph Convolutional Neural Network

####  II-C1 Graph Convolution

(a) Pearson Matrix for Subject Nine

(b) Absolute Pearson Matrix for Subject Nine

(c) Adjacency Matrix for Subject Nine

(d) Laplacian Matrix for Subject Nine

Figure 4: The Pearson, Absolute Pearson, Adjacency, and Laplacian Matrices for
Subject Nine.

In the graph theory, a graph is presented by the graph Laplacian $L$. It is
computed by the Degree Matrix $D$ minus the Adjacency Matrix $A$, i.e.,
$L=D-A$. In this work, Pearson’s matrix $P$ was utilized to measure the inner
correlations among features.

$P_{X,Y}=\frac{E\left(\left(X-\mu_{X}\right)\left(Y-\mu_{Y}\right)\right)}{\sigma_{X}\sigma_{Y}}$
(10)

where $X$ and $Y$ are two variables regarding different features, $\rho_{X,Y}$
is their correlation, $\sigma_{X}$ and $\sigma_{Y}$ are the standard
deviations and $\mu_{X}$ and $\mu_{Y}$ are the expectations. Besides, the
Adjacency Matrix $A$ is recognised as:

$A=|P|-I$ (11)

Where $|P|$ is the absolute of the Pearson’s matrix $P$, and
$I\in\mathbb{R}\textsuperscript{\emph{N}$\times$\emph{N}}$ is an identity
matrix. In addition, the Degree Matrix $D$ of the graph is computed as
follows.

$D_{ii}=\sum_{j=1}^{N}A_{ij}$ (12)

Then, the normalized graph Laplacian is computed as:

$L=D-A=I_{N}-D^{-1/2}AD^{-1/2}$ (13)

It is then decomposed by the Fourier basis
$U=\left[u_{0},\ldots,u_{N-1}\right]\in\mathbb{R}^{N\times N}$. The graph
Laplacian is described as $L=U\Lambda U^{T}$, where
$\Lambda=\operatorname{diag}\left(\left[\lambda_{0},\ldots,\lambda_{N-1}\right]\right)\in\mathbb{R}^{N\times
N}$ are the eigenvalues of $L$. The convolutional operation for a graph is
defined as:

$y=g_{\theta}(L)x=g_{\theta}\left(U\Lambda
U^{T}\right)x=Ug_{\theta}(\Lambda)U^{T}x$ (14)

In which, $g_{\theta}$ is a non-parametric filter. Specifically, the operation
is as follows.

$y_{:,j}^{k+1}=\sigma\left(\sum_{i=1}^{f_{k-1}}Ug_{\theta}(\Lambda)U^{T}x_{:,i}^{k}\right)$
(15)

In which $y^{k}\in\mathbb{R}^{N\times f_{k-1}}$ denotes the signals, $N$ is
the number of vertices of the graph, $f_{k-1}$ and $f_{k}$ are the numbers of
input and output channels respectively, and the $\sigma$ denotes a non-
linearity activation function. What is more, $g_{\theta}$ is approximated by
the Chebyshev polynomials because it is not localized in space and very time-
consuming[33]. The Chebyshev recurrent polynomial approximation is described
as $T_{k}(x)=2xT_{k-1}(x)-T_{k-2}(x),\ T_{0}=1,\ T_{1}=x$. And the filter can
be presented as
$g_{\theta}(\Lambda)=\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{\Lambda})$, in
which $\theta\in\mathbb{R}^{K}$ is a set of coefficients, and
$T_{k}(\tilde{\Lambda})\in\mathbb{R}^{K}$ is the $k^{th}$ order polynomial at
$\tilde{\Lambda}=2\Lambda/\lambda_{max}-I_{n}$, and $I_{n}\in[-1,1]$ is a
diagonal matrix of the scaled eigenvalues. The Convolution can be rewritten
as:

$y=\sum_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{L})x$ (16)

####  II-C2 Graph Pooling

The graph pooling operation can be achieved via the Graclus multilevel
clustering algorithm, which consists of nodes clustering and one-dimensional
pooling[34]. A greedy algorithm is implemented to compute the successive
coarser of a graph and minimize the clustering objective, from which the
normalized cut is chosen[35]. Through such a way, meaningful neighborhoods on
graphs are acquired. Reference [23] proposed to carry out a balanced binary
tree to store the neighborhoods, and a one-dimensional pooling is then applied
for precise dimensionality reduction.

###  II-D Proposed Approach

The presented approach was a combination of the Attention-based BiLSTM and the
GCN as illustrated in Figure 1. The BiLSTM with the Attention mechanism was
presented to derive relevant features from raw EEG signals. During the
procedure, features were obtained from neurons at the FC layer. The GCN was
then applied to classify the extracted features. It was the combination of two
models that promoted and enhanced the decoding performance by a significant
margin compared with existing studies. Details were provided in the following.

First of all, an optimal RNN-based model was explored to obtain relevant
features from raw EEG signals. As detailed in Figure 3, in this work, the
BiLSTM with Attention model was best performed, which achieved 77.86% Global
Average Accuracy (GAA). The input size $\mathbf{x}_{(t)}$ of the model was 64
denoting 64 channels (electrodes) of raw EEG signals. The max time $t$ was
chosen as 64, which was 0.4 seconds’ segment. According to Figure 3, a higher
accuracy has obtained while increasing the number of cells of BiLSTM model. It
should, however, be noted in Figure 3 that when there were more than 256
cells, the loss showed an upward trend, which indicated the concern of
overfitting due to the increment of the model complexity. As a result, 256
LSTM cells (76.67% GAA) were chosen to generalize the model. Meanwhile, it was
apparent that, in Figure 3, as for the linear size of the Attention weights,
the majority of the choices did not make a difference. Thus, 8 neurons, with
79.40% GAA, were applied during the experiments empirically. Comparing Figure
3 and Figure 3, it showed that a compromise solution should be applied, which
took consideration of both performance and input size of the GCN. As a result,
the linear size of 64 (76.73% GAA) was utilized at the FC layer.

Besides, to prevent overfitting, a 25% dropout[36] for the BiLSTM and FC layer
was implemented. The model carried out batch normalization (BN)[37] for the FC
layer, which was activated by the softplus function[38]. And L2 norm with
$1\times 10^{-7}$ coefficient was applied to the Euclidean distance as the
loss function. 1,024 batch sizes were used to maximize the usage of GPU
resources. $1\times 10^{-4}$ learning rate was applied to Adam Optimizer[39].

Furthermore, 2nd order Chebyshev polynomial was applied to approximate
convolutional filters in the experiments. The GCN consisted of six graph
convolutional layers with 16, 32, 64, 128, 256, and 512 filters, respectively,
each followed by a graph max-pooling layer, and a softmax layer derived the
final prediction.

In addition, for the GCN model, before the non-linear softplus activation
function, BN was utilized at all of the layers except the final softmax.
$1\times 10^{-7}$ L2 norm was added to the loss function, which was a cross-
entropy loss. Stochastic Gradient Descent[40] with 16 batch sizes was
optimized by the Adam ($1\times 10^{-7}$ learning rate).

All the experiments above were performed and implemented by the Google
TensorFlow[41] 1.14.0 under NVIDIA RTX2080ti and CUDA10.0.

##  III Results and Discussion

###  III-A Description of Dataset

The data collected from the EEG Motor Movement/Imagery Dataset[42] was
employed in this study. Numerous EEG trials were acquired from 109
participants performing four MI tasks, i.e., imagining left fist (L), right
fist (R), both fists (B), and both feet (F) (21 trials per task). Each trial
is a 4-second experiment’s duration (160 Hz sample rate) with one single
task[17]. In this work, a 0.4 seconds’ temporal segment of 64 channels’
signals, i.e., 64 channels $\times$ 64 time points, was regarded as a sample.
In Section III-B, we used a group of 20 subjects’ data ($S_{1}-S_{20}$) to
train and validate our method. 10-fold cross validation was carried out.
Further, 50 subjects ($S_{1}-S_{50}$) were selected to verify the
repeatability and stability of our approach. In Section III-C, the data set of
individual subjects ($S_{1}-S_{10}$), was utilized to perform subject-level
adaptation. For all the experiments, the dataset was randomly divided into 10
parts, where 90% was the training set, and the left 10% was regarded as the
test set. In Section III-B, the above procedure has carried out for 10 times.
Thus, it left us 10 results of 10-fold cross validation.

###  III-B Group-wise Prediction

It was suggested that the inter-subject variability remains one of the
concerns for interpreting EEG signals[43]. Firstly, a small group size (20
subjects) was adopted for group-wise prediction. In Figure 3, 63.57% GAA was
achieved by the BiLSTM model. After applying the attention mechanism, it
enhanced the decoding performance, which accomplished 77.86% GAA (14.29%
improvement). Further, we employed Attention-based BiLSTM-GCN model in this
work. It attained 94.64% maximum GAA[17] (31.07% improvement compared with the
BiLSTM model), and 93.04% median accuracy from 10-fold cross-validation. Our
method promoted the classification capability under subjects’ variability and
variations by taking the topological relationship of relevant features into
consideration. Meanwhile, as illustrated by Figure 5, the median values of
GAA, Kappa, precision, recall, and F1 Score were 93.04%, 90.71%, 93.02%,
93.01%, and 92.99%, respectively. To the knowledge of the authors, the
proposed method has achieved the best state-of-the-art performance in group-
level prediction. Besides, remarkable results of 10-fold cross-validation have
verified the repeatability and stability. Furthermore, the confusion matrix of
Test One (94.64% GAA) was given in Figure 5. 91.69%, 92.11%, 94.48%, and 100%
accuracy were obtained for each task. It can be observed that our method was
unprecedentedly effective and efficient in detecting human motion intents from
raw EEG signals.

(a) Box Plot for 10-fold cross validation

(b) Confusion Matrix of Test One

Figure 5: Box plot and confusion matrix for 10-fold cross validation.

By grouping signals from additional 30 subjects (in total 50 subjects), the
robustness of the method has been validated in Figure 6.

(a) GAA w.r.t. group-wise prediction

(b) ROC Curve w.r.t. group-wise prediction

Figure 6: GAA and ROC Curve for 20 and 50 subjects, separately.

Towards practical EEG-based BCI applications, it is essential to develop a
robust model to counter serious individual variability[43]. Figure 6
illustrated the GAA of our method through iterations. As listed in Figure 6,
we can see that 94.64% and 91.40% GAA were obtained with regard to the group
of 20 and 50 subjects, respectively. And the Area under the Curves (AUC) were
0.964 and 0.943. Indicated by the above results, the presented approach can
successfully filter the distinctions of signals, even though the data set was
extended. In other words, by increasing the inter-subject variability, the
robustness and effectiveness of the method were evaluated.

The comparison of group-wise evaluation was demonstrated measured by the
maximum of GAA[17] during experiments[44, 17]. Here, we compared the
performance of several state-of-the-art methods in Table I.

TABLE I: Comparison on group-wise evaluation

Related Work | Max. GAA | Approach |  | Num of
---
Subjects
Database
Ma _et al._ (2018) | 68.20% | RNNs | 12 | Physionet Database
Hou _et al._ (2019) | 94.50% | ESI + CNNs | 10
92.50% | 14
This work | 94.64% | Attention-based BiLSTM-GCN | 20

Table I listed the performance of related methods. Hou _et al._ achieved
competitive results. However, our method obtained higher performance (0.14%
accuracy improvement) even with doubling the number of subjects. It can be
found that our approach has outperformed those by giving the highest accuracy
of decoding EEG MI signals.

###  III-C Subject-specific Adaptation

The performance of individual adaptation has witnessed a flourishing
increment[45, 46, 47, 48, 49, 50, 18, 17]. The results of our method on
subject-level adaptation have been reviewed in Table II, and we compared the
results in Table III.

TABLE II: Subject-level Evaluation

No. of Subject | GAA | Kappa | Precision | Recall | F1 Score
---|---|---|---|---|---
1 | 94.05% | 92.06% | 94.20% | 94.32% | 94.16%
2 | 96.43% | 95.19% | 96.06% | 96.06% | 96.06%
3 | 97.62% | 96.79% | 97.33% | 97.08% | 97.18%
4 | 90.48% | 87.34% | 91.30% | 91.11% | 90.42%
5 | 95.24% | 93.61% | 95.96% | 95.06% | 95.38%
6 | 94.05% | 92.02% | 93.40% | 94.96% | 93.66%
7 | 98.81% | 98.40% | 98.81% | 99.07% | 98.92%
8 | 95.24% | 93.60% | 95.39% | 95.04% | 95.19%
9 | 98.81% | 98.39% | 99.11% | 98.68% | 98.87%
10 | 94.05% | 91.98% | 93.39% | 94.70% | 93.61%
Average | 95.48% | 93.94% | 95.50% | 95.61% | 95.35%

Results were given in Table II, from which the highest GAA was 98.81% achieved
by the Subject $S_{7}$ and $S_{9}$, and the lowest was 90.48% by $S_{4}$. On
average, the presented approach can handle the challenge of subject-specific
adaptation. It achieved competitive results, with an average accuracy of
95.48%. Moreover, Cohen’s Kappa coefficient (Kappa), precision, recall, and
f1-score were 93.94%, 95.50%, 95.61%, and 95.35%, respectively. The promising
results above indicated that the introduced method filtered raw EEG signals,
and succeed in classifying MI tasks.

(a) Loss w.r.t. Subject-level Validation

(b) ROC Curve w.r.t. Subject-level Validation

Figure 7: Loss and ROC Curve for subject-level evaluation.

As can be seen from Figure 7, the model has been shown to converge for the
subject-specific adaptation. Receiver Operating Characteristic Curve (ROC
Curve) with its corresponding AUC was visible in Figure 7.

TABLE III: Current studies comparison on subject-level prediction Related Work | Max. GAA | Approach | Database
---|---|---|---
Ortiz-Echeverri _et al._ (2019) | 94.66% | Sorted-fast ICA-CWT + CNNs | BCI Competition IV-a Dataset
Sadiq _et al._ (2019) | 95.20% | EWT + LS-SVM
Taran _et al._ (2018) | 96.89% | TQWT + LS-SVM
Zhang _et al._ (2019) | 83.00% | CNNs-LSTM | BCI Competition IV-2a Dataset
Ji _et al._ (2019) | 95.10% | SVM
Amin _et al._ (2019) | 95.40% | MCNNs
Dose _et al._ (2018) | 68.51% | CNNs | Physionet Database
Hou _et al._ (2019) | 96.00% | ESI + CNNs
This work | 98.81% | Attention-based BiLSTM-GCN

The comparison of subject-level prediction was put forward between the
presented approach and competitive models[45, 46, 47, 48, 49, 50, 18, 17]. The
Attention-based BiLSTM-GCN approach has achieved highly accurate results,
which suggested robustness and effectiveness towards EEG signal processing as
shown in Table III.

The presented approach has improved classification accuracy, and obtained
state-of-the-art results. The reason for the outstanding performance was that
the Attention-based BiLSTM model managed to extract relevant features from raw
EEG signals. The followed GCN model successfully classified features by
cooperating with the topological relationship of overall features.

##  IV Conclusion

To address the challenge of inter-trial and inter-subject variability in EEG
signals, an innovative approach of Attention-based BiLSTM-GCN was proposed to
accurately classify four-class EEG MI tasks, i.e., imagining left fist, right
fist, both fists, and both feet. First of all, the BiLSTM with Attention model
succeeded in extracting relevant features from raw EEG signals. The followed
GCN model intensified the decoding performance by cooperating with the
internal topological relationship of relevant features, which were estimated
from Pearson’s matrix of the overall features. Besides, results provided
compelling evidence that the method has converged to both subject-level and
group-wise predictions, and achieved the best state-of-the-art performance,
i.e., 98.81% and 94.64% accuracy, respectively, for handling individual
variability, which were far ahead of existing studies. The 0.4-second sample
size was proven effective and efficient in prediction compared with
traditional 4s trial length, which means that our proposed framework can
provide time-resolved solution towards fast response. Results on a group of 20
subjects were derived by 10-fold cross validation indicating repeatability and
stability. The proposed method is predicted to advance the clinical
translation of the EEG MI-based BCI technology to meet the diverse demands,
such as of paralyzed patients. In summary, the unprecedented performance with
the highest accuracy and time-resolved prediction were fulfilled via the
introduced feature mining approach.

## Acknowledgments

The authors would like to thank the Brain Team at Google for developing
TensorFlow. We further acknowledge PhysioNet for open source the EEG Motor
Movement/Imagery Dataset to promote the research.

## References

  * [1] C. E. Bouton, A. Shaikhouni, N. V. Annetta, M. A. Bockbrader, D. A. Friedenberg, D. M. Nielson, G. Sharma, P. B. Sederberg, B. C. Glenn, W. J. Mysiw, et al., “Restoring cortical control of functional movement in a human with quadriplegia,” Nature, vol. 533, no. 7602, p. 247, 2016.
  * [2] M. A. Schwemmer, N. D. Skomrock, P. B. Sederberg, J. E. Ting, G. Sharma, M. A. Bockbrader, and D. A. Friedenberg, “Meeting brain–computer interface user performance expectations using a deep neural network decoding framework,” Nature medicine, vol. 24, no. 11, p. 1669, 2018.
  * [3] J. J. Daly and J. R. Wolpaw, “Brain–computer interfaces in neurological rehabilitation,” The Lancet Neurology, vol. 7, no. 11, pp. 1032–1043, 2008\.
  * [4] J. Pereira, A. I. Sburlea, and G. R. Müller-Putz, “Eeg patterns of self-paced movement imaginations towards externally-cued and internally-selected targets,” Scientific reports, vol. 8, no. 1, p. 13394, 2018.
  * [5] M. Mahmood, D. Mzurikwao, Y.-S. Kim, Y. Lee, S. Mishra, R. Herbert, A. Duarte, C. S. Ang, and W.-H. Yeo, “Fully portable and wireless universal brain–machine interfaces enabled by flexible scalp electronics and deep learning algorithm,” Nature Machine Intelligence, vol. 1, no. 9, pp. 412–422, 2019.
  * [6] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” nature, vol. 521, no. 7553, pp. 436–444, 2015.
  * [7] F. Lotte, L. Bougrain, A. Cichocki, M. Clerc, M. Congedo, A. Rakotomamonjy, and F. Yger, “A review of classification algorithms for eeg-based brain–computer interfaces: a 10 year update,” Journal of neural engineering, vol. 15, no. 3, p. 031005, 2018.
  * [8] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” nature, vol. 323, no. 6088, pp. 533–536, 1986\.
  * [9] T. Zhang, W. Zheng, Z. Cui, Y. Zong, and Y. Li, “Spatial–temporal recurrent neural network for emotion recognition,” IEEE transactions on cybernetics, vol. 49, no. 3, pp. 839–847, 2018.
  * [10] P. Wang, A. Jiang, X. Liu, J. Shang, and L. Zhang, “Lstm-based eeg classification in motor imagery tasks,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 26, no. 11, pp. 2086–2095, 2018\.
  * [11] T.-j. Luo, F. Chao, et al., “Exploring spatial-frequency-sequential relationships for motor imagery classification with recurrent neural network,” BMC bioinformatics, vol. 19, no. 1, p. 344, 2018.
  * [12] N. F. Güler, E. D. Übeyli, and I. Güler, “Recurrent neural networks employing lyapunov exponents for eeg signals classification,” Expert systems with applications, vol. 29, no. 3, pp. 506–514, 2005.
  * [13] X. Zhang, L. Yao, S. S. Kanhere, Y. Liu, T. Gu, and K. Chen, “Mindid: Person identification from brain waves through attention-based recurrent neural network,” Proceedings of the ACM on Interactive, Mobile, Wearable and Ubiquitous Technologies, vol. 2, no. 3, p. 149, 2018.
  * [14] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [15] K. Fukushima, “Neocognitron: A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position,” Biological cybernetics, vol. 36, no. 4, pp. 193–202, 1980.
  * [16] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, et al., “Gradient-based learning applied to document recognition,” Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.
  * [17] Y. Hou, L. Zhou, S. Jia, and X. Lun, “A novel approach of decoding EEG four-class motor imagery tasks via scout ESI and CNN,” Journal of Neural Engineering, vol. 17, p. 016048, feb 2020.
  * [18] H. Dose, J. S. Møller, H. K. Iversen, and S. Puthusserypady, “An end-to-end deep learning approach to mi-eeg signal classification for bcis,” Expert Systems with Applications, vol. 114, pp. 532–542, 2018.
  * [19] M. Henaff, J. Bruna, and Y. LeCun, “Deep convolutional networks on graph-structured data,” arXiv preprint arXiv:1506.05163, 2015.
  * [20] J. Bruna, W. Zaremba, A. Szlam, and Y. Lecun, “Spectral networks and locally connected networks on graphs,” in International Conference on Learning Representations (ICLR2014), CBLS, April 2014, pp. http–openreview, 2014.
  * [21] D. K. Duvenaud, D. Maclaurin, J. Iparraguirre, R. Bombarell, T. Hirzel, A. Aspuru-Guzik, and R. P. Adams, “Convolutional networks on graphs for learning molecular fingerprints,” in Advances in neural information processing systems, pp. 2224–2232, 2015.
  * [22] M. Niepert, M. Ahmed, and K. Kutzkov, “Learning convolutional neural networks for graphs,” in International conference on machine learning, pp. 2014–2023, 2016.
  * [23] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in Advances in neural information processing systems, pp. 3844–3852, 2016.
  * [24] X.-h. Wang, T. Zhang, X.-m. Xu, L. Chen, X.-f. Xing, and C. P. Chen, “Eeg emotion recognition using dynamical graph convolutional neural networks and broad learning system,” in 2018 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pp. 1240–1244, IEEE, 2018.
  * [25] T. Zhang, X. Wang, X. Xu, and C. P. Chen, “Gcb-net: Graph convolutional broad network and its application in emotion recognition,” IEEE Transactions on Affective Computing, 2019.
  * [26] T. Song, W. Zheng, P. Song, and Z. Cui, “Eeg emotion recognition using dynamical graph convolutional neural networks,” IEEE Transactions on Affective Computing, 2018.
  * [27] Z. Wang, Y. Tong, and X. Heng, “Phase-locking value based graph convolutional neural networks for emotion recognition,” IEEE Access, vol. 7, pp. 93711–93722, 2019.
  * [28] K. Cho, B. van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder–decoder for statistical machine translation,” in Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), (Doha, Qatar), pp. 1724–1734, Association for Computational Linguistics, Oct. 2014.
  * [29] D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” arXiv preprint arXiv:1409.0473, 2014\.
  * [30] K. Xu, J. Ba, R. Kiros, K. Cho, A. Courville, R. Salakhudinov, R. Zemel, and Y. Bengio, “Show, attend and tell: Neural image caption generation with visual attention,” in International conference on machine learning, pp. 2048–2057, 2015.
  * [31] Z. Yang, D. Yang, C. Dyer, X. He, A. Smola, and E. Hovy, “Hierarchical attention networks for document classification,” in Proceedings of the 2016 conference of the North American chapter of the association for computational linguistics: human language technologies, pp. 1480–1489, 2016\.
  * [32] J. K. Chorowski, D. Bahdanau, D. Serdyuk, K. Cho, and Y. Bengio, “Attention-based models for speech recognition,” in Advances in neural information processing systems, pp. 577–585, 2015.
  * [33] D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” Applied and Computational Harmonic Analysis, vol. 30, no. 2, pp. 129–150, 2011.
  * [34] I. S. Dhillon, Y. Guan, and B. Kulis, “Weighted graph cuts without eigenvectors a multilevel approach,” IEEE transactions on pattern analysis and machine intelligence, vol. 29, no. 11, pp. 1944–1957, 2007.
  * [35] J. Shi and J. Malik, “Normalized cuts and image segmentation,” Departmental Papers (CIS), vol. 22, no. 8, pp. 888–905, 2000.
  * [36] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: a simple way to prevent neural networks from overfitting,” The journal of machine learning research, vol. 15, no. 1, pp. 1929–1958, 2014\.
  * [37] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” arXiv preprint arXiv:1502.03167, 2015.
  * [38] R. H. Hahnloser, R. Sarpeshkar, M. A. Mahowald, R. J. Douglas, and H. S. Seung, “Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit,” Nature, vol. 405, no. 6789, p. 947, 2000.
  * [39] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * [40] T. Zhang, “Solving large scale linear prediction problems using stochastic gradient descent algorithms,” in Proceedings of the twenty-first international conference on Machine learning, p. 116, ACM, 2004.
  * [41] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., “Tensorflow: A system for large-scale machine learning,” in 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), pp. 265–283, 2016\.
  * [42] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley, “Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals,” Circulation, vol. 101, no. 23, pp. e215–e220, 2000.
  * [43] H. Tanaka, “Group task-related component analysis (gtrca): a multivariate method for inter-trial reproducibility and inter-subject similarity maximization for eeg data analysis,” Scientific Reports, vol. 10, no. 1, pp. 1–17, 2020.
  * [44] X. Ma, S. Qiu, C. Du, J. Xing, and H. He, “Improving eeg-based motor imagery classification via spatial and temporal recurrent neural networks,” in 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 1903–1906, IEEE, 2018.
  * [45] C. J. Ortiz-Echeverri, S. Salazar-Colores, J. Rodríguez-Reséndiz, and R. A. Gómez-Loenzo, “A new approach for motor imagery classification based on sorted blind source separation, continuous wavelet transform, and convolutional neural network,” Sensors, vol. 19, no. 20, p. 4541, 2019\.
  * [46] M. T. Sadiq, X. Yu, Z. Yuan, Z. Fan, A. U. Rehman, G. Li, and G. Xiao, “Motor imagery eeg signals classification based on mode amplitude and frequency components using empirical wavelet transform,” IEEE Access, vol. 7, pp. 127678–127692, 2019.
  * [47] S. Taran and V. Bajaj, “Motor imagery tasks-based eeg signals classification using tunable-q wavelet transform,” Neural Computing and Applications, vol. 31, no. 11, pp. 6925–6932, 2019.
  * [48] R. Zhang, Q. Zong, L. Dou, and X. Zhao, “A novel hybrid deep learning scheme for four-class motor imagery classification,” Journal of neural engineering, vol. 16, no. 6, p. 066004, 2019.
  * [49] N. Ji, L. Ma, H. Dong, and X. Zhang, “Eeg signals feature extraction based on dwt and emd combined with approximate entropy,” Brain sciences, vol. 9, no. 8, p. 201, 2019.
  * [50] S. U. Amin, M. Alsulaiman, G. Muhammad, M. A. Bencherif, and M. S. Hossain, “Multilevel weighted feature fusion using convolutional neural networks for eeg motor imagery classification,” IEEE Access, vol. 7, pp. 18940–18950, 2019.


# Motor strategies and adiabatic invariants: The case of rhythmic motion in
parabolic flights

N. Boulanger1<EMAIL_ADDRESS>F. Buisseret2,3
<EMAIL_ADDRESS>V. Dehouck1,4<EMAIL_ADDRESS>F.
Dierick2,5,6<EMAIL_ADDRESS>O. White4<EMAIL_ADDRESS>1 Service de Physique de l’Univers, Champs et Gravitation, Université de Mons,
UMONS Research Institute for Complex Systems, Place du Parc 20, 7000 Mons,
Belgium  2 CeREF, Chaussée de Binche 159, 7000 Mons, Belgium  3 Service de
Physique Nucléaire et Subnucléaire, Université de Mons, UMONS Research
Institute for Complex Systems, 20 Place du Parc, 7000 Mons, Belgium  4
Université de Bourgogne INSERM-U1093 Cognition, Action, and Sensorimotor
Plasticity, Campus Universitaire, BP 27877, 21078 Dijon, France  5 Faculté des
Sciences de la Motricité, Université catholique de Louvain, 1 Place Pierre de
Coubertin, 1348 Louvain-la-Neuve, Belgium  6 Centre National de Rééducation
Fonctionnelle et de Réadaptation – Rehazenter, Laboratoire d’Analyse du
Mouvement et de la Posture (LAMP), Luxembourg, Grand-Duché de Luxembourg

###### Abstract

The role of gravity in human motor control is at the same time obvious and
difficult to isolate. It can be assessed by performing experiments in variable
gravity. We propose that adiabatic invariant theory may be used to reveal
nearly-conserved quantities in human voluntary rhythmic motion, an individual
being seen as a complex time-dependent dynamical system with bounded motion in
phase-space. We study an explicit realization of our proposal: An experiment
in which we asked participants to perform $\infty-$ shaped motion of their
right arm during a parabolic flight, either at self-selected pace or at a
metronome’s given pace. Gravity varied between $0$ and $1.8$ $g$ during a
parabola. We compute the adiabatic invariants in participant’s frontal plane
assuming a separable dynamics. It appears that the adiabatic invariant in
vertical direction increases linearly with $g$, in agreement with our model.
Differences between the free and metronome-driven conditions show that
participants’ adaptation to variable gravity is maximal without constraint.
Furthermore, motion in the participant’s transverse plane induces trajectories
that may be linked to higher-derivative dynamics. Our results show that
adiabatic invariants are relevant quantities to show the changes in motor
strategy in time-dependent environments.

##  I Introduction

Gravity obviously plays a role in human motor control, but it is not easy to
isolate. On one hand, the perception and internal representation of gravity in
the brain results from a widely distributed sensory process, including the
vestibular system, vision, and somatosensory information White et al. (2020).
On the other hand, the constant immersive nature of the body in the Earth
gravitational field calls for holistic experimental approaches in contrast to
focused interventions on isolated body parts. In these settings, participants
are exposed to variable gravito-inertial fields generated by human centrifuges
or parabolic flights. The rationale behind these experimental approaches rest
on Einstein’s equivalence principle (see e.g. Einstein (1916); Wald (1984)):
Physics in an accelerated spacecraft (1g) is undistinguishable from physics on
the ground. It does, however, not mean that the brain does not attempt to
identify the possible different sources that give rise to the same
consequences. Beyond the brain’s role, it is well-known that variable gravity
induces various changes in human physiology, typically at the cardiovascular
Aubert et al. (2016), neural White et al. (2016) and musculoskeletal levels
Lang et al. (2017). In the present work, we do not focus on a particular
physiological aspect but rather adopt a global (bio)mechanical point of view.
Our main goal is to show that tools derived from mechanics may lead to
conserved measurable parameters that the brain can exploit to plan and execute
actions instead of relying on an estimate of gravity acquired through a noisy
and distributed process.

Finding such conserved quantities demands to know the underlying dynamics,
which is not an easy task in human motion where a same – even simple – action
requires sometimes very different motor commands. For instance, consider
reaching for a cup of coffee on the breakfast table or in an aircraft subject
to turbulences, or drinking while seated or while walking. Planning efficient
actions is challenging for the brain. It is therefore natural that the central
nervous system relies on constants in this jungle of variability. Previous
research demonstrated that some classes of actions result from an optimization
process in which movements features are taken into account, such as minimizing
jerk Viviani and Flash (1995), metabolic cost Alexander (1997); Berniker et
al. (2013) or maintaining a nearly constant mechanical energy in level walking
Cavagna et al. (1976). Here we study rhythmic motion in time-changing gravity
as a peculiar case of time-dependent dynamical system with bounded motion. The
most powerful tools known so far to study such systems are adiabatic
invariants Landau and Lifchitz (1988); Henrard (1998); Jose and Saletan
(1998). They have been applied in a wide range of applications such as plasma
physics Notte et al. (1993) and cosmology Cotsakis et al. (1998). In
biomechanics, several studies have showed the invariance of the action
variable in time-independent conditions Kugler and Turvey (1987); Kugler et
al. (1990); Turvey et al. (1996); Kadar et al. (1993). In a previous work, we
went a step further and proposed to use adiabatic invariants to study the
changes in the motion of upper arm rhythmic movements about the elbow at free
pace and amplitude in a centrifuge where the perceived gravity’s intensity
changed stepwise from 1 to 3 $g$ and from 3 back to 1 $g$ Boulanger et al.
(2020). The direction of $\vec{g}$ was unchanged. It appeared that the
behaviour of the adiabatic invariant $I=\frac{1}{2\pi}\oint_{\Gamma}P\,dQ$
computed from the latter one-degree-of-freedom motion was compatible with a
theoretically predicted linear increase with $g$ Boulanger et al. (2019).

An obvious direction to generalize the framework of Ref. Boulanger et al.
(2020) is that of motions involving more than 1 degree of freedom. We first
show in Sec. II that a linear link between the adiabatic invariant and $g$ is
expected for any potential energy, only assuming a separable dynamics. Our
model is then applied to analyse the data of Ref. White et al. (2008) in Sec.
III. In this last work, participants were asked to continually perform an
$\infty-$shaped trajectory during a parabolic flight. The three-dimensional
kinematics of the hand has been recorded and adiabatic invariants can be
computed from it. Gravity during the parabolas varied between 0 and $1.8$ $g$.
We present our results in Sec. IV and discuss them in Sec. V.

##  II The model

###  II.1 Adiabatic invariants in variable gravity

We assume that human voluntary rhythmic motion in variable gravity $g(t)$ may
be seen as a dynamical system with bounded motion in phase space, described by
a separable Hamiltonian
$H(I_{\alpha},\theta^{\alpha},\lambda(t))=\sum^{D}_{\alpha=1}H_{j}(I_{\alpha},\theta^{\alpha},\lambda(t))$.
The Hamiltonian depends on action-angle variables $I_{\alpha}$ and
$\theta^{\alpha}$ respectively and on a time-dependent function $\lambda(t)$,
accounting for the modifications induced by $g(t)$. The various ingredients
underlying the latter assumption deserve further comments. First, Hamiltonian
dynamics being the most powerful formulation of classical Mechanics, it is
rather natural to adopt a Hamiltonian approach. This being said, not every
dynamics is Hamiltonian: A sufficient criterion for a Hamiltonian to exist is
that the total time-derivative applied to the Poisson bracket of two functions
defined on phase-space obeys the Leibniz rule (Jose and Saletan, 1998, Chapter
5). The latter criterion cannot a priori be checked from our experimental
setup: We will a posteriori confirm that the observed phase-space trajectories
are not incompatible with a Hamiltonian dynamics. Second, the dynamics has a
priori no reason to be separable. The main interest we have in assuming the
separability is that it allows a clear separation between vertical and
horizontal directions (with respect to $\vec{g}$), the dynamics in the
vertical direction being intuitively the most strongly impacted by the
variations of $g\,$. Again, the separability hypothesis will only be checked a
posteriori by observing the phase-space trajectories, see next Section.

It can then be shown that (Landau and Lifchitz, 1988, Eqs (50,10)-(50,11))

$\displaystyle\dot{I}_{\alpha}$ $\displaystyle=$ $\displaystyle-\frac{\partial
H}{\partial\theta^{\alpha}}=-\left(\frac{\partial\Lambda}{\partial\theta^{\alpha}}\right)_{I_{\alpha},\lambda}\dot{\lambda}\
,$ (1a) $\displaystyle\dot{\theta}^{\alpha}$ $\displaystyle=$
$\displaystyle\frac{\partial H}{\partial
I_{\alpha}}=\omega_{\alpha}+\left(\frac{\partial\Lambda}{\partial
I_{\alpha}}\right)_{\theta^{\alpha},\lambda}\dot{\lambda}\ ,$ (1b)

where the partial derivatives have to be computed while keeping constant the
indexed variables and where $\omega_{\alpha}$ are the motion’s frequencies.
The function $\Lambda$ is the action of the system; it is sufficient for our
purpose to state that it is a periodic function of the angle variables. Hence,
according to Landau and Lifchitz (1988),
$\Lambda=\sum^{+\infty}_{\ell_{1}=-\infty}\dots\sum^{+\infty}_{\ell_{D}=-\infty}{\rm
e}^{i\vec{\ell}\cdot\vec{\theta}}\Lambda_{\vec{\ell}}\ $ with
$\Lambda_{\vec{\ell}}\in\mathbb{C}$,
$\vec{\ell}=(\ell_{1},\dots,\ell_{D})\in\mathbb{Z}^{D}$ and

$\frac{\partial\Lambda}{\partial\theta^{\alpha}}=\sum^{+\infty}_{\ell_{1}=-\infty}\dots\sum^{+\infty}_{\ell_{D}=-\infty}i\ell_{\alpha}{\rm
e}^{i\vec{\ell}\cdot\vec{\theta}}\Lambda_{\vec{\ell}}\ .$ (2)

We moreover assume that $\lambda=\lambda_{0}+\epsilon g(t)$ with $\epsilon
g(t)\ll\lambda_{0}$, i.e. that the modifications induced by variable gravity
may be computed at first-order in $\epsilon$. Equation (1a) therefore leads to

$\dot{I}_{\alpha}=-\epsilon\,\dot{g}\sum^{+\infty}_{\ell_{1}=-\infty}\dots\sum^{+\infty}_{\ell_{D}=-\infty}i\ell_{\alpha}{\rm
e}^{i\vec{\ell}\cdot\vec{\theta}}\Lambda_{\vec{\ell}}\ ,$ (3)

or

$\frac{dI_{\alpha}}{dg}=-\epsilon\,\sum^{+\infty}_{\ell_{1}=-\infty}\dots\sum^{+\infty}_{\ell_{D}=-\infty}i\ell_{\alpha}{\rm
e}^{i\vec{\ell}\cdot\vec{\theta}}\Lambda_{\vec{\ell}}\ .$ (4)

Let us define the times $t_{n}$ such that the values $\theta^{\alpha}(t_{n})$
are all equal (modulo $2\pi$). The existence of $t_{n}$ is guaranteed for
periodic dynamics such as the one that we consider here. Then it can be said
that

$\left.\frac{dI_{\alpha}}{dg}\right|_{t=t_{n}}=I_{1}\Rightarrow
I_{\alpha}(t_{n})=I_{\alpha;0}+I_{\alpha;1}\,g(t_{n}),$ (5)

with $I_{\alpha;0}$ and $I_{\alpha;1}$ two real numbers such that
$|I_{\alpha;1}/I_{\alpha;0}|\ll 1$. Equation (5) defines our model: The
adiabatic invariant is expected to behave linearly in $g$ when computed at a
given position in the consecutive cycles performed.

###  II.2 Definition of $I_{\alpha}$

The action variables are defined from positions $(Q^{\alpha})$ and momenta
$(P_{\alpha})$ degrees of freedom as follows:

$I_{\alpha}=\frac{1}{2\pi}\oint_{\Gamma_{\alpha}}P_{\alpha}\,dQ^{\alpha}\ ,$
(6)

where $\Gamma_{\alpha}$ is the projection of the bounded trajectory in the
plane $(Q^{\alpha},P_{\alpha})$ for fixed $\alpha\,$. Note that, with a
kinetic energy of the standard form
$E_{c}\sim\sum^{D}_{\alpha=1}\dot{Q}^{\alpha\,2}$, one is led to a form for
the adiabatic invariant which is straightforward to compute:

$I_{\alpha}(t)\sim\int^{t+T}_{t}\dot{Q}^{\alpha\,2}(u)\,du\,,$ (7)

with $T$ the period of the phase-space cycle $\Gamma_{\alpha}$ starting at
$t$.

This last equation provides a way to compute the adiabatic invariant from
experimental data provided $Q^{\alpha}(t)$ is known, which is not so obvious
since a mathematical description of voluntary human motion may involve higher
derivative dynamics, see e.g. Nelson (1983); Hogan (1984); Hagler (2015). Two
cases should therefore be considered. First, the motion’s dynamics does not
involve higher derivative terms. In this case $Q^{\alpha}$ may directly be
identified to, say, one anatomical landmark’s trajectory $x^{\alpha}(t)$, and
$P_{\alpha}\sim\dot{x}^{\alpha}$. Second, the motion’s dynamics is a higher-
derivative one. Then $Q^{\alpha}$ and $P_{\alpha}$ can, in principle, be
computed from $x^{\alpha}(t)$ but their definition is more involved. We refer
the interested reader to the case of Pais-Uhlenbeck oscillator Pais and
Uhlenbeck (1950), that is a higher-derivative generalization of standard
harmonic oscillator for which adiabatic invariants can be analytically
computed Boulanger et al. (2019).

##  III The experiment

###  III.1 Parabolic flights

During parabolic flights, participants were asked to continually perform an
$\infty$-shaped trajectory oriented crosswise to the body around two virtual
obstacles situated 3 m in front of them. An optoelectronic device (OptoTrak
3020 system, Northern Digital, Waterloo, Ontario, Canada) recorded the
position of three infrared LEDs placed on the object with a resolution of 0.1
mm. A three-dimensional accelerometer fixed on the floor of the aircraft
recorded its acceleration. The two synchronized acquisition systems recorded
parameters at a sampling rate of 200 Hz. During a parabola, the aircraft
performs a series of manoeuvres to allow for changes of effective gravity.
This allows to run experiments at 0 (microgravity), 1 and approximately 1.8
$g$ (hypergravity), albeit for a short time. The micro and hyper gravity
phases last around 20 s with transition periods shorter than 5 s. Typical
plots of the motion performed and of the $g(t)$-profile are shown in Fig. 1.
The cartesian frame we use is also displayed.

Figure 1: (a) Typical plot of the $\infty-$shaped motion in frontal plane
$(x,y)$ during one parabola. A participant in FREE condition has been chosen.
(b) $g(t)$ profile during the same parabola. (c) The Cartesian frame is
displayed.

Participants executed $\infty-$shaped movements in two conditions. In the free
condition (FREE), the motion was self-paced. Four participants, totalling 24
parabolas, performed the motion in FREE condition. In the metronome condition
(METRO), participants had to adopt 1.5-second constant pace prompted by a
metronome. Seven participants, totalling 42 parabolas, performed the motion in
METRO condition. Before starting the parabolic flights, participant’s health
was assessed by their individual National Centres for Aerospace Medicine as
meeting the requirement ”Jar Class II” for parabolic flight. No participant
reported sensory or motor deficits and they all had normal or corrected-to-
normal vision. All participants gave their informed consent to participate in
this study and the procedures were approved by the European Space Agency (ESA)
Safety Committee and by the local ethics committee. Their motion was recorded
during 6 consecutive parabolas. We refer the reader to Ref. White et al.
(2008) for a more detailed presentation of the experiment.

Figure 2: (a) Typical speed-position plot of the motion in $(x,v_{x})$ plane
during several cycles, same participant as Fig. 1. (b) and (c) Same data in
the $(y,v_{y})$ and $(z,v_{z})$ planes respectively. The crosses show the
trajectory obtained from Eq. (8) with $A_{0}=0.092$ m, $A_{1}=0.065$ m,
$A_{2}=0.060$ m, $\omega=2\pi$ rad/s, $\phi_{1}=0$ rad and $\phi_{2}=-1$ rad.

###  III.2 Phase-space trajectories and action variables

The speeds $v_{\alpha}=\dot{x}^{\alpha}$ are first computed from the positions
$x^{\alpha}$ recorded by the optoelectronic device through a finite
differentiation. Typical speed-position plots are shown in Fig. 2.

The $x$ and $y$ directions show (quasi)-periodic trajectories of elliptic
type, compatible with a standard Hamiltonian dynamics. Hence we proceed as
follows to compute the action variables. First, we identify $Q^{\alpha}$ to
$x^{\alpha}$ and $P_{\alpha}$ to $\dot{Q}^{\alpha}$ – up to an arbitrary mass
scale that is set equal to 1 kg. Second, the beginning and end of each cycle
$\Gamma_{\alpha}$ in phase-space plane $(Q^{\alpha},P_{\alpha})$ are computed.
The end of the cycle starting at $t$ is chosen as the time $t^{\star}$ which
is the smallest time after $t$ at which the euclidean distance between
$(Q^{\alpha}(t^{\star}),P_{\alpha}(t^{\star}))$ and
$(Q^{\alpha}(t),P_{\alpha}(t))$ is minimal. Once $t^{\star}$ is identified,
the adiabatic invariant is computed by quadrature from Eq. (4):
$I_{\alpha}(t)=\int^{t+t^{\star}}_{t}\dot{Q}^{\alpha\,2}(u)\,du\,$. Then, to
apply our model, only adiabatic invariants corresponding to a given value of
the angle variable have to be collected. We only consider the instants at
which $P_{\alpha}=0$ and $Q^{\alpha}$ was maximal since they are easily
identified.

As can be seen in Fig. 2, the trajectory in $(z,v_{z})$ plane intersects
itself during one cycle. The underlying dynamic is therefore called _non-
autonomous_ , in the theory of dynamical systems.Participant’s motion in the
$z-$direction actually contains two distinct frequencies. One is the whole
$\infty$-shaped movement’s pulsation, say $\omega$, and the other is a
forward-backward oscillation at $2\omega$. A trajectory of the form

$z(t)=A_{0}+\sum_{j=1,2}A_{j}\sin(j\omega t+\phi_{j})$ (8)

has the qualitative features of what is observed in the $(z,v_{z})$ for
appropriate values of the real constants $\omega$, $A_{i}$ and $\phi_{i}$.

Several effective models may produce trajectories such as (8). (1) An
oscillator with pulsation $\omega$ plus an external, time-dependent, periodic
force with pulsation $2\omega$. A textbook example is the Duffing oscillator.
The complete set of solutions of a forced, non-harmonic oscillator is unknown
a priori but some special solutions are known that perfectly match the
observed motion. (2) A system of coupled harmonic oscillators oscillating near
its equilibrium position provided that the frequencies of some normal modes
are equal to $\omega$ and $2\omega$. The dynamics of the various joints of the
arm may be approximated by such a system. (3) A second-order Pais-Uhlenbeck
harmonic oscillator. As shown in Appendix A, an autonomous higher-derivative
oscillator may actually mimic the dynamics of one peculiar degree of freedom
in a system of coupled oscillators. An obvious advantage in resorting to Pais-
Uhlenbeck harmonic oscillators is that the effective dynamics in the
$z-$direction would be autonomous and that adiabatic invariants are well-
defined once phase-space is properly built from the position degree of freedom
and its time derivatives Boulanger et al. (2019). Hence, our model can in
principle be adapted to the $z-$direction. However, such higher-derivative
adiabatic invariants involve not only $\dot{z}$ but at least $\ddot{z}$ and
$\dddot{z}$ Boulanger et al. (2019). The experimental precision reached in the
measurement of $z$ does not allow for a reliable computation of those higher
derivatives from our experimental data and we chose not to make further
computations as far as the $(z,v_{z})$ plane is concerned.

##  IV Action variables in terms of gravity: The results

We have linearly fitted $I_{\alpha}$ versus $g$ for each available parabola in
order to check whether Model (5) is observed at an individual level or not.
The variable $g$ refers to the average value of $g(t)$ within the considered
cycle $\Gamma_{\alpha}$. The parameters of the fit are $I_{1;\alpha}$ (slope),
$r_{\alpha}$ (Pearson’s correlation coefficient) and $I_{0;\alpha}$
(intercept). A two-way ANOVA may be performed on the parameters of the fit
with factors condition (FREE or METRO) and parabola number (1 to 6). The
latter factor is introduced to check whether a learning effect is present or
not during the consecutive parabolas a given participant has experienced. The
ANOVA was performed using SigmaPlot software (v.11.0, Systat Software, San
Jose, CA, United States of America), with significance level 0.05. It appears
that no significant effect of the parabola number can be found in the fit
parameters which means that the model parameters are stable over time, with
values set from the outset. Interactions between condition and parabola number
are not significant either. However, the condition has a significant impact on
$I_{0;\alpha}$, $I_{1;y}$ and $r_{y}$, as shown in Table 1.

From Table 1, the following global features of participant’s motion can be
deduced. (a) The action variable $I_{x}$ does not show a well-defined linear
behaviour versus $g$: both the slopes and the Pearson’s correlation
coefficients are comparable with 0. The action variable $I_{x}$ can be
considered as constant with $g$, although it is significantly lower in METRO
condition than in the FREE condition. Let us note that the $x-$direction is
orthogonal to gravity while the $y-$direction is aligned with gravity (see
frame in Fig. 1). This observation may intuitively explain why there is no
trend versus $g$ for the $x-$dynamics. (b) However, the trend of $I_{y}$ vs
$g$ is compatible with model (5) for positive $I_{1;y}$ well smaller than the
intercept $I_{0;y}$. It is coherent with our initial assumption to work at
first-order in $g$. Furthermore, the slope and the intercept are significantly
lower in METRO condition than in the FREE one. (c) The clearest linear trend
is observed for the action variable $I_{y}$ in the FREE condition. An example
of linear fit is displayed in Fig. 3.

Figure 3: Adiabatic invariant $I_{y}$ versus $g$ computed from experimental data in the FREE condition for the same participant as in Fig. 1 (points), compared to the best linear fit of the form (5). Pearson’s correlation coefficients is also indicated. $I_{y}$ has been normalized so that its average value is $1$ at $1g$.  Direction | Condition |  $I_{1;\alpha}$ (kg.m) | $r_{\alpha}$ |  $I_{0;\alpha}$ (J.s)
---|---|---|---|---
$x$ | FREE | $[-0.014,0.005]$ | $[-0.150,0.082]$ | $[0.327,0.371]$
| METRO | $[-0.010,0.004]$ | $[-0.182,0.109]$ | $[0.177,0.205]$
| $p$ | 0.884 | 0.969 | $<0.001$
$y$ | FREE | $[0.017,0.033]$ | $[0.291,0.518]$ | $[0.187,0.296]$
| METRO | $[0.005,0.14]$ | $[0.117,0.332]$ | $[0.106,0.120]$
| $p$ | $<0.001$ | $0.014$ | $<0.001$
Table 1: 95 % confidence intervals for the slopes $I_{1;\alpha}$, Person’s
correlation coefficients $r_{\alpha}$ and intercepts $I_{0;\alpha}$ obtained
through the fit (5) of the computed $I_{\alpha}$ vs $g$ in each parabola for
all conditions. The $p-$ values of the ANOVA for the effect of condition are
also given.

The global trend of $I_{y}$ vs $g$ can be observed by averaging $I_{y}$ over
participants by condition (FREE and METRO) and by gravity condition, i.e. by
gathering computed adiabatic invariants into bins of 0.1 $g$, ranging from 0
to 1.8 $g$. Only the bins containing more than 10 points were finally kept.
This threshold is arbitrary but avoids almost empty bins in the fast
transition regions between 0 and 1 $g$ and between 1 and 1.6 $g$. The results
of this analysis are displayed in Fig. 4. The observed trends are compatible
with $I_{y}^{FREE}=0.210+0.015\,g$ and $I_{y}^{METRO}=0.106+0.015\,g$; an
ANCOVA further shows that the slopes are not significantly different with
condition ($p=0.994$), while the intercepts significantly depend on condition
($p<0.001$). Finally, it is also worth highlighting the fact that some bins
(e.g. $[1.3;1.4[$) capture values of gravity in the ascending but also
descending parts of the parabolic profile.

Figure 4: Adiabatic invariant $I_{y}$ versus $g$ computed from experimental
data in the FREE (grey points) and METRO (yellow points) conditions. A linear
fit is given (solid lines) with its 95% confidence interval (colored bands) in
each condition. Pearson’s correlation coefficients are also indicated. The
1-$g$ bin is marked with a vertical dashed line.

##  V Discussion of the results

The trajectory in the $z-$direction contains two distinct frequencies leading
to intersecting trajectories in the $(z,v_{z})$ plane, see Fig. 2. Effective
dynamics in this plane therefore cannot be described by a time-independent
standard dynamics: Either higher-derivatives or time-dependent forces have to
be included. We believe that the appearance of such features could be
explained by shoulder biomechanical constraints. The main shoulder movements
required to execute the $\infty-$shaped movement are abduction and adduction.
During shoulder abduction-adduction movements, rotations are usually observed
Assi et al. (2016). It is well known that external rotation of the shoulder is
adopted during abduction to clear the major tubercule of humerus from beneath
acromion for preventing impingement Peat (1986); Hurov (1986). The movement
strategy spontaneously chosen by participants is therefore not located on a
single plane with constant $z$, leading to the observed nontrivial pattern.As
previously said, the current experimental accuracy along the $z-$axis does not
allow for a more detailed study of a potential higher-derivative effective
dynamics. Note that it has already been successfully conjectured that a
higher-derivative action principle such as $S=\int\dddot{x}^{2}\,dt$ – a jerk-
based cost function – may constrain non-rhythmic voluntary human motion Hogan
(1984). However, such an action principle does not lead to periodic solutions,
that is why a Pais-Uhlenbeck oscillator seems more relevant to us. In the
motion we observe, the two frequencies have an integer ratio, therefore
stability of the motion is not guaranteed in such a resonant case Boulanger et
al. (2019). It has recently been understood that extra interaction terms may
stabilize periodic solutions of resonant higher-derivative oscillators
Kaparulin et al. (2020): We hope to investigate the applicability of such
models to human rhythmic motion in future works.

By definition, the adiabatic invariant in $x$ and $y-$directions is
proportional to

$I_{\alpha}\sim T\left\langle E_{c,\alpha}\right\rangle,$ (9)

where $T$ is a given cycle duration in the vertical direction – the duration
of the cycle starting at the same time is twice that value in the
$x-$direction –, and where $\left\langle E_{c,\alpha}\right\rangle$ is the
averaged kinetic energy on the considered cycle. The protocol of White et al.
(2008) is such that $T^{FREE}<T^{METRO}$: The pace imposed by the metronome
was chosen to be slower than participants’ spontaneously chosen paces. Since,
at given $g$, the adiabatic invariant in FREE condition is always larger than
in METRO condition, it can be concluded that $\left\langle
E_{c,\alpha}^{FREE}\right\rangle>\left\langle
E_{c,\alpha}^{METRO}\right\rangle$. The smaller kinetic energy in METRO
condition thus follows from the fact that participants have to move slower
than in FREE condition in order to follow metronome’s pace. Note that
$I_{1;y}^{METRO}<I_{1;y}^{FREE}$: The extra constraint imposed by the
metronome actually prevents the participant from optimally adapting his/her
motion when $g$ is changing, assuming that the optimal motor strategy is
reached in FREE condition. It is also known that $T$ is a decreasing function
of $g$ in either METRO or FREE conditions White et al. (2008). Since
$I_{1;y}>0$, $\left\langle E_{c,y}\right\rangle$ has to be an increasing
function of $g$: The participant’s arm moves with higher typical vertical
speed at higher values of $g$.

Definition (9) explicitly makes appear the links between adiabatic invariant
and kinetic energy. Another interpretation of the adiabatic invariant,
focusing on external forces, is relevant to clarify the influence of gravity
on it. To this aim, the virial theorem may be used to state that

$I_{\alpha}\sim T\left\langle F_{\alpha}x^{\alpha}\right\rangle,$ (10)

where $F_{\alpha}$ is an external force acting on the point-like object whose
trajectory is $x^{\alpha}$. The latter force should involve muscular forces as
well as gravity. One can reasonably assume that
$F_{\alpha}=F_{0\alpha}+F_{1\alpha}g$, in coherence with the previously found
linear trend of $I_{\alpha}$ vs $g$. A priori, $F_{1y}\gg F_{1x}$ since
gravity’s influence should mostly concern the vertical direction. Changes
induced by $F_{1x}$ are probably unnoticeable up to our current experimental
precision.

Around $1g$, $I_{y}$ is lower than expected from the 95% confidence interval
of the linear fit. In that familiar environment, participants “know” the most
economic strategy when they are allowed to move freely in Earth’s gravity. In
METRO condition, that drop in $I_{y}$ is not observable: The extra constraint
imposed by the metronome does not allow participants to follow that optimal
strategy. Furthermore, our results also reveal that microgravity is a special
case. While the linear fit holds true for the whole explored gravitational
values ($[0g;1.8g]$), there is a significant gap between 0.3 $g$ and 0.7 $g$
in our data. Hypogravity values are not explored. Our study again reveals that
0 $g$ acts as a singular value for the brain White et al. (2008). Finally,
adiabatic invariants do not behave like parameters measured in most motor
control investigations. Indeed, while errors in reaching movements perturbed
by force fields require tens of trials to vanish Shadmehr and Mussa-Ivaldi
(1994), safety margins in object manipulation in altered gravity need an
exposure to 6 parabolas to decrease to normal values Augurelle et al. (2003)
and some behaviours in conflicting force-fields or visuomotor rotations do
even not adapt at all Cothros et al. (2008), adiabatic invariants seem to be
set to their nominal values from the outset.

##  VI Concluding comments

To conclude this work, it is worth linking them to well-known frameworks in
motor control.

Participants have many more kinematic degrees of freedom than necessary to
fulfill the demanded task, i.e. the $\infty-$shaped movement. The coordination
of kinematically redundant systems was formulated by Bernstein as the degrees
of freedom problem Bernstein (1967). The main difficulty of Bernstein’s
problem is that the nervous system must conciliate two apparently conflicting
abilities: (1) the realization of a movement from the choice of one among an
infinite number of motor patterns; (2) the absence of univocal relationship
between the movement realized and motor patterns used, known as motor
equivalence. Although it remains unclear as to whether and how the brain can
estimate adiabatic invariants, such quantities puts constraints on the allowed
strategies, i.e. strategies keeping $I_{\alpha}$ invariant at constant $g$. In
this picture, an increase (decrease) in $I_{\alpha}(g)$ may be related to an
increase (decrease) of the allowed motor patterns.

When subjects are free to point to a target, they automatically scale movement
duration with movement amplitude and choose a trade off between movement speed
and accuracy to touch the target. It is known as Fitts’s law Fitts (1954). The
adiabatic invariant is the area of a closed trajectory in phase space:
$I_{\alpha}\sim A_{\alpha}\,v^{{\rm max}}_{\alpha}$, with $A_{\alpha}$ and
$v^{{\rm max}}_{\alpha}$ the amplitude and maximal speed of the movement in
the direction $\alpha$ respectively. Its invariance at given $g$ implies that,
if maximal speed increases (decreases), amplitude decreases (increases). In
our experiment, the maximal speed is an obvious measure of movement’s speed,
and the amplitude can be seen as an index of precision. Indeed,the instruction
given to the participant is to avoid 2 targets by turning around. Thus, if
amplitude decreases (increases), the participant increases (decreases) the
chances of hitting the target, and he/she is less (more) precise. The
adiabatic invariant can then be seen as an explicit realization of the speed-
accuracy-trade off scenario. The modification of its value with $g$ actually
changes the acceptable values of maximal speed and amplitude involved in this
trade off.

In summary, our results indicate that adiabatic invariants deserve a
particular attention in biomechanical approaches of human motion. They are
indeed able to capture one individual’s reaction to time-dependent external
conditions, even in extreme cases such as variable gravity. Adiabatic
invariants seem very robust in this context. Further studies are now needed to
clarify the links between adiabatic invariant theory and celebrated motor
control paradigms such as speed-accuracy-trade off.

Acknowledgements This research was supported by grants from Prodex and IAP;
Belgian Federal Office for Scientific, Technical, and Cultural Affairs; Fonds
Spécial de Recherche; and Canadian Space Agency Contract 9F007-033026.

*

##  Appendix A Higher derivative dynamics and rhythmic motion

Let us consider a system with $N$ degrees of freedom $x^{\alpha}$ described by
the Lagrangian

$\displaystyle
L=\frac{1}{2}g_{\alpha\beta}(x)\dot{x}^{\alpha}\dot{x}^{\beta}-U(x^{\gamma})$
(11)

with $g_{\alpha\beta}$ the components of a real, symmetric and positive-
definite matrix $G$ that we call the kinetic matrix. Note that it is not
necessarily constant and may depend on the dynamical variables. If necessary
after a translation of the origin of the coordinates, we may assume that
$x^{\alpha}=0$ ($\forall~{}\alpha$) is an equilibrium position:
$\left.\frac{\partial U}{\partial x^{\alpha}}\right|_{x^{\gamma}=0}=0\,$. Such
a Lagrangian may model the motion of several joints, the potential energy $U$
being an a priori complicated function of the degrees of freedom. If only
small oscillations around equilibrium position are considered, the equations
of motion read
$\ddot{x}^{\alpha}+\gamma^{\alpha\delta}U_{\delta\beta}x^{\beta}=0\,$, with
$\gamma^{\alpha\beta}$ the components of the inverse of the matrix $G_{0}$ of
components $\gamma_{\alpha\beta}:=g_{\alpha\beta}(0)\,$. In other words, one
has $\gamma^{\alpha\delta}\gamma_{\delta\beta}=\delta^{\alpha}_{\beta}\,$. One
defines the potential matrix $U$ with components

$U_{\alpha\beta}=\left.\frac{\partial^{2}U}{\partial x^{\alpha}\partial
x^{\beta}}\right|_{x^{\delta}=0}.$ (12)

Solving the eigenequation
$V^{\alpha}{}_{\beta}\,\xi^{\beta}_{a}=\lambda_{a}\,\xi^{\alpha}_{a}$ for the
matrix $V=G_{0}^{-1}U$ with components
$V^{\alpha}{}_{\beta}=\gamma^{\alpha\delta}\,U_{\delta\beta}\,$, with
$a=1,\dots,N\,$, allows to solve the equations of motion in terms of the
normal coordinates $Q^{a}(t)$:

$x^{\alpha}=\xi^{\alpha}_{a}\,Q^{a}(t)\
\qquad\mbox{with}\qquad\ddot{Q}^{a}=-\lambda_{a}\,Q^{a}$ (13)

and without any summation over the index $a\,$. Therefore, for the given
dynamical system described by the variables $x^{\alpha}$, any small
oscillatory motion about the minimum of the potential in configuration space
can therefore be decomposed as a linear combination of elementary oscillations
along the normal modes, each one at the frequency $\nu_{a}=\omega_{a}/2\pi\,$,
where $\lambda_{a}=\omega_{a}^{2}\,$. In particular, for appropriate initial
conditions it is possible to only excite the normal mode $Q^{a}(t)$ for a
given value of the index $a\,$. The dynamical system as a whole will then
oscillate at the single frequency $\nu_{a}\,$, without exciting the modes
$Q^{b}(t)$ with $b\neq a\,$. The interested reader may find a detailed
discussion about small oscillations around an equilibrium position in (Landau
and Lifchitz, 1988, Chapter 5), or in (Arnold, 1989, Part 2, Chapter 5) for a
more precise mathematical formulation.

Note that, from the datum of the normal modes $Q$’s with their frequencies
$\nu$’s, one can go back and access the information contained in the kinetic
and potential matrices $G\,$ and $U\,$. This is because the normal modes are
orthogonal with respect to the metric $G\,$, and using the latter metric
together with the eigenvalues of $V:=G^{-1}U\,$ gives $U\,$ up to a reordering
of the dynamical variables $x^{\alpha}\,$.

In a human rhythmic motion, if the participant is asked to perform a periodic
motion, say with the forearm, one observes that the projection of the motion
of the hand along the three spatial directions gives rise to a very small set
of frequencies that are all integer multiples of a fundamental one. In this
description, we neglect the quasi-periodic motion of the forearm due to
physiological noise. Of course, the forearm is a very complicated system with
dozens of components linked in a complicated fashion, giving rise to a
configuration space $\mathbb{Q}$ of very large dimension $N\,$. In principle
(but not in practice), it is possible to describe it by a Lagrangian of the
form (11) and there will $N$ normal modes $Q^{a}$’s with possible degeneracies
in the frequencies. The observed motion of the forearm of the participants is,
instead, very simple and degenerate.

Instead of trying to find the realistic Lagrangian description (11) of the
forearm, from the sole observation of a very limited set of forced periodic
motions with distinct frequencies $\omega_{a}$’s, motions that we view as
analogous to the distinct normal modes of a dynamical system, we propose an
effective model whose purpose it to reproduce those “normal modes” without any
diagonalisation of any potential matrix $V\,$. The operator

$F=\Pi^{n}_{a=1}\left(1+\frac{1}{\omega^{2}_{a}}\frac{d^{2}}{dt^{2}}\right)$
(14)

is such that $Fx^{\alpha}=0\,$ _for all_ $\,\alpha$ because $FQ^{a}=0$ for all
$a\,$, by construction. Here, by an abuse of notation we have denoted by
$Q^{a}(t)$ the $n\ll N$ simple and pure harmonic modes observed in the
participant’s motion. The latter describes the motion in a configuration space
of very large dimension $N\,$ whereas we effectively reduce the dynamics to a
configuration space of dimension $n$ way smaller than $N\,$. Therefore, in our
effective description of the motion based on a specific set of harmonic
oscillations observed in the forearm’s motion, every single dynamical variable
$x^{\alpha}$ for $\alpha$ fixed obeys the equation of motion of a Pais-
Uhlenbeck oscillator whose Lagrangian reads
$L_{P-U}=-\frac{1}{2}\,x^{\alpha}\,F\,x^{\alpha}$ Pais and Uhlenbeck (1950).
If at least two frequencies $\omega_{a}$ are different, the effective dynamics
of a given degree of freedom can be mimicked by a particular solution of the
equations of motion of a higher-derivative harmonic oscillator.

## References

  * White et al. (2020) O. White, . J. Gaveau, L. Bringoux, and F. Crevecoeur, J. Neurophysiol. 124, 4 (2020).
  * Einstein (1916) A. Einstein, Annalen Phys. 49, 769 (1916).
  * Wald (1984) R. M. Wald, _General relativity_ (Chicago Univ. Press, Chicago, IL, 1984).
  * Aubert et al. (2016) A. Aubert et al., npj Microgravity 2, 16031 (2016).
  * White et al. (2016) O. White et al., npj Microgravity 2, 16023 (2016).
  * Lang et al. (2017) T. Lang, J. Van Loon, S. Bloomfield, L. Vico, A. Chopard, J. Rittweger, A. Kyparos, D. Blottner, I. Vuori, R. Gerzer, et al., npj Microgravity 3, 8 (2017).
  * Viviani and Flash (1995) P. Viviani and T. Flash, Journal of Experimental Psychology: Human Perception and Performance 21, 32 (1995).
  * Alexander (1997) R. Alexander, Biol Cybern. 76, 97 (1997).
  * Berniker et al. (2013) M. Berniker, M. K. O’Brien, K. P. Kording, and A. A. Ahmed, PLOS ONE 8, 1 (2013).
  * Cavagna et al. (1976) G. Cavagna, H. Thys, and A. Zamboni, J Physiol. 262, 639 (1976).
  * Landau and Lifchitz (1988) L. Landau and E. Lifchitz, _Physique théorique Tome 1 : Mécanique_ (E. MIR, Moscow, 1988).
  * Henrard (1998) J. Henrard, _The Adiabatic Invariant in Classical Mechanics_ (Dessy, 1998), pp. 60–73.
  * Jose and Saletan (1998) J. Jose and E. Saletan, _Classical dynamics: a contemporary approach_ (Cambridge Univ. Press, Cambridge, 1998).
  * Notte et al. (1993) J. Notte, J. Fajans, R. Chu, and J. S. Wurtele, Phys. Rev. Lett. 70, 3900 (1993).
  * Cotsakis et al. (1998) S. Cotsakis, R. L. Lemmer, and P. G. L. Leach, Phys. Rev. D 57, 4691 (1998).
  * Kugler and Turvey (1987) P. Kugler and M. Turvey, _Information, Natural Law, and the Self-Assembly of Rhythmic Movement_ (London: Routledge, London, 1987).
  * Kugler et al. (1990) P. Kugler, M. Turvey, R. Schmidt, and L. Rosenblum, Ecological Psychology 2, 151 (1990).
  * Turvey et al. (1996) M. Turvey, K. Holt, J. Obusek, et al., Biol. Cybern. 74, 107 (1996).
  * Kadar et al. (1993) E. Kadar, R. Schmidt, and M. Turvey, Biol. Cybern. 68, 421 (1993).
  * Boulanger et al. (2020) N. Boulanger, F. Buisseret, V. Dehouck, F. Dierick, and O. White, Phys. Rev. E 102, 062403 (2020), URL https://link.aps.org/doi/10.1103/PhysRevE.102.062403.
  * Boulanger et al. (2019) N. Boulanger, F. Buisseret, F. Dierick, and O. White, Eur. Phys. J. C79, 60 (2019), eprint 1811.07733.
  * White et al. (2008) O. White et al., J Neurophysiol 100 (2008).
  * Nelson (1983) W. Nelson, Biol. Cybern. 46, 135 (1983).
  * Hogan (1984) N. Hogan, J. Neurosci. 4, 2745 (1984).
  * Hagler (2015) S. Hagler (2015), URL arXiv:1509.06981.
  * Pais and Uhlenbeck (1950) A. Pais and G. E. Uhlenbeck, Phys. Rev. 79, 145 (1950).
  * Assi et al. (2016) A. Assi, Z. Bakouny, M. Karam, A. Massaad, W. Skalli, and I. Ghanem, Human movement science 50, 10 (2016).
  * Peat (1986) M. Peat, Physical therapy 66, 1855 (1986).
  * Hurov (1986) J. Hurov, Journal of hand therapy 22, 328 (1986).
  * Kaparulin et al. (2020) D. S. Kaparulin, S. L. Lyakhovich, and O. D. Nosyrev, Phys. Rev. D 101, 125004 (2020), eprint 2003.10860.
  * Shadmehr and Mussa-Ivaldi (1994) R. Shadmehr and F. Mussa-Ivaldi, J Neurosci. 14, 3208 (1994).
  * Augurelle et al. (2003) A. Augurelle, M. Penta, O. White, and J. Thonnard, Exp Brain Res. 148, 533 (2003).
  * Cothros et al. (2008) N. Cothros, J. Wong, and P. L. Gribble, PLOS ONE 3, 1 (2008).
  * Bernstein (1967) N. Bernstein, _The Coordination and Regulation of Movements_ (Oxford: Pergamon Press, Oxford, 1967).
  * Fitts (1954) P. M. Fitts, Journal of Experimental Psychology 47, 381 (1954).
  * Arnold (1989) V. I. Arnold, _Mathematical methods of classical mechanics_ , Graduate Text in Mathematics (Springer-Verlag, 1989), second edition ed.