Experimental setup
Four Inertial Measurement Units (IMU) sensors (Three-Space Sensor™, Yost Labs Inc.) were strapped to the arms and forearms of the participant via Velcro® bands as in Fig. 1, panel a. Sensors were placed bilaterally, one proximally to the biceps, the other below the elbow on the ulnar bone. Pitch and Roll angles of each sensor were recorded using a custom written C++ code in Simulink®, yielding a total of 8 sensor inputs at a frequency of 50 Hz. The BoMI consisted of the mxq linear map H : Wm → Vn from the 8D vector space of arm motions, q ∈ Wm, into the 2D vector space of cursor position on a computer screen, p ∈ Vn, as defined by Eq. (1):
$$ \boldsymbol{p}=\mathrm{H}\left(\boldsymbol{s}-E\left[\boldsymbol{s}\right]\right) $$
(1)
The BoMI map H0 was defined via a simple calibration procedure. Each participant was instructed to move their arms spontaneously without task specifications in a comfortable range for 35 s. Subsequently, the dimensionality of the calibration dataset was reduced via Principal Component Analysis by extracting the first two components that expressed the greatest variance. H0 was thus defined to be the composition of the first two eigenvectors u1,u2 of the calibration covariance matrix. It is important to note here that H0 : W0m → V0n, where the vector space W0m is the particular subspace of possible arm motions that have been observed during calibration.
This transformation allowed to map the first and second eigenvector onto width (x-axis) and height (y-axis) coordinates on the screen respectively. The cursor displacement was then rescaled by the eigenvalues λ1, λ2 and the width w and height h of the desired workspace:
$$ {\boldsymbol{p}}_{\mathbf{1}}={\boldsymbol{u}}_{\mathbf{1}}\boldsymbol{s}\frac{w}{\sqrt{\lambda_1}},\kern0.5em {\boldsymbol{p}}_{\mathbf{2}}={\boldsymbol{u}}_{\mathbf{2}}\boldsymbol{s}\frac{h}{\sqrt{\lambda_2}}; $$
(2)
The cursor was visualized on the screen as a red circle of 1 cm diameter whose position was controlled by different arm postures according to H0.
Protocol
The study protocol was divided in three main phases. During the first phase, participants calibrated the interface using random arm motions, as described in the previous section. In the second phase, they practiced a reaching task to 5 target locations, shown in Fig. 1 – panel b, uniformly distributed on an ellipse with main axes 21.52 cm and 13.50 cm. The target location was displayed as a yellow circle of 1.5 cm diameter. The reaching task duration was set to 18 min. The practice was split into 4 blocks of respective duration 3′ – 5′ – 5′ – 5′. Participants were allowed to rest in between blocks as long as needed to reduce fatigue. During the third phase we evaluated participants’ ability to control the BoMI in an untrained novel task (tracking task). Participants were required to follow a target circle that moved along three complete laps of the infinity-shaped trajectory in Fig. 1 – panel b. The whole protocol lasted about 1 h and is summarized in Fig. 1 – panel c.
Reaching task
In the reaching task, participants were instructed to reach as many targets as possible in a given slot of time in order to maximize their score. For every target reached, the score would be incremented by 1 point. Furthermore, after the initial 12 trials, if the target was reached in less than 1.5 s, the score would be incremented by an extra point. Participants were given the score as a feedback at all times on the top-right edge of the screen. In order to reach a target, they had to hold the cursor within 1 cm of the target center for at least a second. One trial consisted of a complete reaching movement to a target.
Participants were also informed that on some trials the cursor position might not have been visible at the time of the target appearance. In this condition, which we called blind trials, participants were instructed to perform a quick motion that they thought would bring them to overlap with the target and hold the posture. After 2 s into the trial, the cursor would reappear to allow for correcting any offset. The absence of visual feedback in the first 2 s of movement allowed us to separate the effect of feedforward predictive control of the cursor movement from the corrective actions due to feedback mechanisms driven by visual error. Blind trials were interspersed pseudo-randomly in the target set.
Tracking task
During the tracking task, the target moved clockwise along an infinity-shaped trajectory at a rate of 60 s/lap. Participants were instructed to keep in the closest possible vicinity of the target, ideally overlapping with it at all times. Anytime the distance from the target fell below 5 cm, the target paused its motion.
In summary, participants engaged in three different tasks with the BoMI – a random movements task, a reaching task, and a tracking task. Random movements of the arms were used to calibrate the interface map H0. Participants only trained in the reaching task, and the tracking task was used to assess generalization of the reaching skills.
Subjects and conditions
Twenty volunteers (age 19–45, 8F) with no history of neuromotor impairment participated in the study procedures. All participants gave their informed consent prior to the test. All procedures were carried out in accordance with the ethical standards of the IRB and with the Declaration of Helsinki and Northwestern University IRB approved all human involvement in the study (IRB #STU00057856).
Participants were split in two groups and assigned to one of the following study conditions, as in Fig. 1 – panel c:
-
i.
Constant map (group C, N = 10): participants were assigned a constant BoMI map H0, defined by an individual calibration procedure (see section 2.1). Therefore, each individual had a personalized map which did not change throughout the session.
-
ii.
Adaptive Map (group A, N = 10): the initial BoMI map H0 was defined by a calibration procedure as section 2.1. During the reaching task, an adaptive procedure (see section 2.4) continuously updated the body-cursor mapping in time as follows. The BoMI map was initialized with H0 and kept constant for 3 trials. From the 4th trial until 1′ to the end of the reaching task the map was iteratively updated. The final minute of training was then preformed keeping the map constant to the latest available map at the final time step K, HK. The same map HK was also used during tracking.
Strategy for adapting the map
Let us assume that learning to control a BoMI entails a progressive reshaping of movement coordination towards a set of “synergies” – or covariation patterns – that lie on a manifold that has reduced dimension. This assumption is consistent with previous studies that showed that learning novel redundant visuomotor mappings is accompanied by a progressive reduction in complexity of movement coordination [25, 27,28,29]. Moreover, we assume that the movement manifold is a linear subspace of all the possible movement vectors (each vector is one sensor orientation axis). In this sense, the BoMI implements one particular subspace of movement vectors, that we call H.
According to our first hypothesis (HP1), once BoMI users have identified a suitable movement subspace for controlling the interface, learning in terms of improving movement performance would proceed rapidly. Therefore, we devised a strategy for facilitating the identification of the BoMI manifold through its online reorientation to better match the natural movement manifold of the user as this evolves through exposure to particular tasks.
By construction, the BoMI subspace H represents the directions of maximum movement covariance over a set of calibration data. In practice, we expect movement covariance to have a different structure from calibration while the user learns to control the interface [26]. One possible way to compensate for the difference in distributions and to account for non-stationarity introduced by learning is recomputing the full covariance matrix using a batch of new data or a moving window on the data [30]. This operation, despite being beneficial to user’s performance [21], is however expensive in terms of memory consumption and computation time, which grows with the square of the batch size and the dimensionality of the data. Moreover, having to recompute the full covariance matrix through time with data available on the fly not only is expensive, but also unnecessary, given that we are only interested in the first few eigenvectors of movement covariance.
A more efficient and practical solution is given by iterative methods that approximate the eigenvectors of the covariance matrix without the need of having all the data readily available [31]. Weng et al. [32] proposed an iterative formulation that produces an average of a series of stochastic estimates of the eigenvectors of the underlying process covariance, without directly approximating the full covariance matrix. We summarized the method hereafter.
Let us consider a stochastic process with zero mean (without loss of generality) that generates a sequence of N independent identically distributed d-dimensional data s(k),k = 1, . . N, with unknown dxd sample covariance matrix Γ(N). For a specific sample k in the sequence, Γ(k) = E[ s(k)s(k)T], hence \( {\Gamma}^{(N)}=\frac{1}{N+1}{\sum}_{k=1}^N{\Gamma}^{(k)} \). Knowing that the eigenvectors u and eigenvalues λ of the covariance matrix must satisfy the equality Γu=λu, we can further write an iterative expression for v(N) = λu(N) given Γ(N):
$$ {\mathbf{v}}^{(N)}=\frac{1}{N+1}{\sum}_{k=1}^N\ {\boldsymbol{s}}^{(k)}{\boldsymbol{s}}^{(k)T}{\mathbf{u}}^{(k)} $$
(3)
Where u(N) = v(N)/‖v(N)‖, ‖v(N)‖ = λ. Finally, we can re-write equation (3) recursively for obtaining an update rule for the eigenvectors of Γ to apply iteratively for every new sample s(k) as follows:
$$ {\mathbf{v}}^{(k)}=\frac{k-1}{k}{\mathbf{v}}^{\left(k-1\right)}+\frac{1}{k}\left(\ {\boldsymbol{s}}^{(k)}{\boldsymbol{s}}^{(k)T}\bullet \frac{{\mathbf{v}}^{\left(k-1\right)}}{\left\Vert {\mathbf{v}}^{\left(k-1\right)}\right\Vert}\right) $$
(4)
In order to recover orthogonality of the estimated eigenvectors, we can apply the Gram-Schmidt orthogonalization procedure on the vectors v(k) obtained by equation (4). In the case of non-centered data, the sample mean can be estimated as \( {\mathbf{m}}^{(k)}=\frac{k-1}{k}{\mathbf{m}}^{\left(k-1\right)}+\frac{1}{k}\ {\boldsymbol{s}}^{(k)} \).
From equation [7], we see that the estimate at the current time stamp k is composed of a memory element (the value of v at the previous step) and a novel element that depends on the mismatch between the sample covariance and the axes of the map and defines the direction of attraction of v. The effect of this computation is a progressive rotation of the original map towards the estimated directions of variance expressed by the user. The averaging component allows the formulation to be free of tuning parameters, as opposed to stochastic PCA approximations [31]. The drawback is that the estimate becomes progressively insensitive to new data as k increases, making the algorithm unsuitable for tracking non-stationary processes.
In order to cope with non-stationarity in the movement distribution, we have previously proposed a modified version of the iterative algorithm, that introduces the use of a fixed learning rate or amnesic parameter [32]. Considering a process with non-zero mean m, Equation (4) becomes as follows:
$$ {\mathbf{v}}^{(k)}=\left(1-\eta \right){\mathbf{v}}^{\left(k-1\right)}+\eta \left(\ {\boldsymbol{s}}^{(k)}-{\boldsymbol{m}}^{(k)}\right){\left({\boldsymbol{s}}^{(k)}-{\boldsymbol{m}}^{(k)}\right)}^{\boldsymbol{T}}\frac{{\mathbf{v}}^{\left(k-1\right)}}{\left\Vert {\mathbf{v}}^{\left(k-1\right)}\right\Vert } $$
(5)
$$ {\boldsymbol{m}}^{(k)}=\left(1-\eta \right){\boldsymbol{m}}^{\left(k-1\right)}+\eta\ {\boldsymbol{s}}^{(k)} $$
(6)
Equation (5) provides an updated estimate for the sample mean. The purpose of the amnesic parameter 0 < η < 1 in equation [5, 6] is to introduce an exponential downweighing effect over past samples so as to allow capturing new variation in the data. The estimate is subtracted to the data value at each iteration before computing the eigenvector update.
The amnesic parameter η can also be seen as the learning rate for map adaptation. The higher η, the faster new samples will steer the BoMI subspace towards recent experience. Schmitt et al. [33] suggested that in processes monitoring applications that make use of recursive PCA with exponential downweighing, the range of reasonable values for the amnesic parameter is commonly restricted to 0.9 < 1 − η < 0.9999, since values η > 0.01 might introduce instability in the model. On the other hand, values of η too close to 0 might lead to poor monitoring performance, which in our case translates in the algorithm having very low sensitivity to changes in covariance occurring at time scales faster than a few hours.
Convergence of non-adaptive formulation of the algorithm of equation (4) was provided in [34] under the assumption of Gaussian stationary processes. Even if convergence of equation (5) has not been mathematically proven yet, we have verified through simulations that the algorithm in equation (5) converges to the eigenvector of a series of template covariance matrixes in time for learning rates between 10− 4 and 10− 2. Learning rates in the range 0.001 to 0.003 yielded the best reconstruction in terms of estimated eigenvalues for a sample size of 500 to 3000 data points.
The problem of adopting learning rates that do not converge to zero with time has previously been addressed in other iterative formulations of the eigenvector problem since in practice the number will be nonzero due to numerical approximations and round-off errors [33]. In particular, [35, 36] provided proof of global convergence of the iterative PCA method of Oja and Karhunen [37], from which the algorithm of Weng et al. [32] was derived, in the case of constant learning rates. Other methods that employ recursive PCA showed that these algorithms are generally well behaved for a reasonable choice of the learning rate given the process dynamics. How to choose the optimal rate remains, however, an open problem.
In order to minimize sample-to-sample correlation when running the iterative PCA of equations (5, 6) online, we provided as input at each update step k the average of 5 samples of sensor data for which the sensor speed exceeded a threshold. This allowed suspending the update in conditions when no movement of the user’s arms was detected. After averaging, 500 and 3000 update steps correspond to roughly 1 and 5 min of continuous movement data when running at 50 Hz.
When selecting the learning rate, we assumed that the time scale of adaptation would differentially affect the development of a movement strategy. In particular, changing the map later into the training would not significantly affect the movement strategy, while a change early on would probably trigger a more significant co-adaptation. A previous study by Orsborne et al. [38, 39] suggested that the time scale of 1–2 min might be ideal when updating Kalman-filter based closed-loop decoders in brain machine interfaces in cases when the initial performance of the decoder is poor. Given the similarity of intents in our study, we followed a similar approach choosing η = 0.002 (time constant of approximately a minute).
Outcome measures
During the reaching practice, we quantified performance over trials in terms of two metrics:
Reaching Time – the time from the first appearance of a new target to when the target has been reached successfully. This metrics was only computed when the cursor was visible for the whole duration of the trial. A reduction of Reaching Time indicates improvement in performance.
Reaching Error – the Euclidean distance between the target position and the cursor position after 2 s from target appearance. Since the distance across target location was not constant in our design, we normalized the Euclidean distance by the relative distance from the start target position and final target to be reached in a trial. Therefore, the measure has a lower bound of 0, in case the target was reached at or before 2 s. This metric was computed also for blind trials. In this case, the metrics quantifies the distance between the target and the cursor position at the time of cursor re-appearance.
We quantified performance in the tracking task in terms of average Tracking Error, which is computed as the Euclidean distance between the position of the moving target on the screen and the cursor position at a given instant of time.
We then computed a series of metrics to characterize the movement distribution (distribution of sensor values) in absolute terms (#1), relative to another distribution (#2) or relative to the BoMI manifold (#3):
#1. Planarity Index – the variance accounted for by the first 2 principal components of the sensor values over a certain period of time. Previous studies have shown that when learning a novel sensorimotor transformation, individuals tend to progressively distribute movement variance along dimensions that are task-relevant [25]. Given that the dimensionality of the task is only 2, we expected our participants to increase planarity in time. After computing the covariance matrix of the 8 sensor values S over a certain interval of time or trials, and denoting λi the eigenvalues of S, the Planarity index is computed as follows:
$$ \mathrm{Planarity}=\frac{\lambda_1+{\lambda}_2}{\sum_{i=1}^8{\lambda}_i} $$
#2. Subspace Angle – expresses the angle between two hyperplanes embedded in a higher dimensional space. If the two hyperplanes are parallel, the angle is zero. If they are orthogonal, the angle is 90 deg. Given two rectangular matrices 8 × 2 H0 and H1 whose columns are vectors of unit length, the Subspace Angle is computed as:
$$ \mathrm{Subspace}\ \mathrm{Angle}={\cos}^{-1}\left(\left|{H_0}^T{H}_1\right|\right)\bullet \frac{180}{\pi } $$
#3. Variance accounted for by the BoMI (VAF) – is the percentage of overall sensors variance accounted for by the 2D BoMI manifold and quantifies how much the movement distribution of the user is aligned with the 2D distribution defined by the BoMI. For the control group the BoMI manifold is a hyperplane defined by the calibration. For the Adaptive group, the BoMI hyperplane changes iteratively throughout the whole session as a function of the user’s movement distribution. The ideal value of 100% indicates that the movement variance is completely captured by the 2D BoMI manifold. VAF < 100% indicate that some components of variance project in the BoMI null space. Given the 8 × 8 covariance matrix of the sensor values S and the 2 × 2 covariance matrix of the corresponding cursor positions on the screen C, we computed the VAF as the ratio between the trace of C and the trace of S:
$$ \mathrm{VAF}=\frac{tr\left(\mathrm{C}\right)}{tr\left(\mathrm{S}\right)}\bullet 100 $$
Data analysis and statistics
Sensors values and cursor positions were resampled in order to ensure a constant sample time and smoothed using a 3rd order Savitsky-Golay filter with 3 Hz cutting frequency. During reaching, trajectories in cursor and sensor spaces were segmented into trials excluding the holding phase on the target position. All analyses were carried out using MATLAB® R2017a.
Normality of the data was assessed using the Anderson-Darling test. When normality was met, we applied a t-test for comparing the performance of the two groups in the same block of trials or unit of time. Whenever we were interested in comparing the performance of one or both groups through time, we used a repeated measure analysis of variance (ANOVA) design with group as a between factor and repetitions as within factor. When sphericity was not met according to the Mauchly test, we reported the Greenhouse-Geisser correction for the p-values. If normality was not met, we applied a Mann Whitney U test for comparing the median of two populations. The repeated measures ANOVA was replaced by a Friedman’s ANOVA. The level of significance was set to 0.05. Post-hoc tests were carried out using the Tukey-Kramer procedure for multiple comparisons.
To compare the performance of participants in the beginning and the end of the reaching practice we computed the average Reaching Time and Reaching Error at 2 s during the first 20 and last 20 trials with complete vision of the cursor and Reaching Error at 2 s in the first and last 10 blind trials. To compare the evolution of movement distribution through time, we computed Planarity, Subspace Angle, and VAF metrics using the data from a 60 s time window either in the beginning or the end of the reaching practice. Whenever we were comparing performance metrics in trials with metrics qualifying movement distributions (i.e. Planarity), we used movement data within each trial that was included in the analysis.
In order to compute learning curves of Reaching Time and Planarity we followed two procedures.
#1: Computing learning curves across trials: here we were interested in obtaining a single model that could account for differences in learning dynamics within each experimental group. Hence, we used a non-linear mixed-effect model that accounts for inter-subject variability as a random factor affecting model parameters to fit a single exponential function to the Reaching Time and Planarity in first 90 trials with visual feedback (nlmefit function, Matlab 2017a Statistics and Machine Learning toolbox).
#2: Computing learning curves in time: in order to extract the dynamics of performance through time, we fitted an exponential model to the Reaching Time at each trial against the time elapsed from the beginning of the reaching practice to the end of the corresponding trial. Since we are here considering all the trials completed by a participant in the complete feedback condition, participants who completed greater amount of trials exhibited learning with different time-scales across trials of practice. For this reason, we fit two models to each individual’s data, a single and double exponential function using a non-linear least square regression (fit function, MATLAB 2017a) and we selected the model that yielded the higher adjusted R squared for each participant to extract the learning time constant(s).