Introduction

Modeling the brain response to motor tasks has been widely studied for different conditions such as movement control [1,2,3,4], motor tasks classification [5, 6], and sequential movements [7,8,9]. These works focus on the voluntary movement of the subject, however, brain activity can be unintentionally modified by some disturbances coming from the environment. In consequence, although several studies exist about brain activity related to motor control tasks, the analysis when limbs are externally manipulated has not been deeply explored.

Studying brain activity in healthy people when external forces disturb limbs is key for gaining an understanding about how the brain reacts and how brain activity relates to the stimulus received. Hence, understanding the relationship between brain activity and external wrist manipulation of healthy people is the first step to find potential biomarkers that could be used in the future to examine the variations when brain injury occurs in people with motor deficits, tracing the desired path in the rehabilitation process.

Therefore, the underlying aim of the present work is to model the general response of brain caused by wrist manipulation of young and healthy subjects to first comprehend the behavior of the healthy brain, and then find out a common relationship among individuals via the impulse response of the system. Understanding the link between the brain response and the external perturbation in the wrist of healthy people would allow the identification of potentially helpful biomarkers in the examination of motor impairment in the upper limbs [10].

In this sense, some proposals have been developed. For instance, Nozari et al. [11] built a model based on auto-regressive techniques and locally linear neuro-fuzzy (LLNF) networks, which achieved the highest performance. However, the model’s parameters are subject-specific, which makes it a non-general model. A NARMAX model structure is proposed by Gu [12], however, the different time delays of the model terms do not fully explain the physiological pathways of brain activity elicited by each realization.

It is important to realize that the previously mentioned models only consider information related to displacement and velocity of the movement. However, there is evidence that acceleration is involved in proprioceptive sensory processes, particularly in sensorimotor responses to external perturbations [13].

To fill these gaps and contribute to a better understanding of the described problem, the proposal is based on a nonlinear delay differential embedding to model the brain response caused by wrist manipulation. A delay-differential embedding transforms a time series into a geometric structure composed of successive derivatives and delayed versions of the signal [14]. Therefore, this model could estimate the non-linear terms of the sensorimotor system while approaching to a general framework of the brain response elicited by wrist manipulation when the same time delays are considered among the parameters of the induced movement. Moreover, the delay differential equations-based models have been already used to distinguish patients with Parkinson’s disease from healthy individuals, suggesting that time delays are important to detect motor impairments [15].

To achieve this aim, the 4TU Research Dataset [16] from the Delft University of Technology is used in this work. In the experiment, ten healthy volunteers are exposed to a robotic system that causes a controlled wrist movement, and the resulting brain response is measured via electroencephalogram (EEGFootnote 1). To validate the proposed model, the cross-correlation regarding the Variance Accounted For (VAFFootnote 2) and the Pearson Correlation Coefficient are computed. The best averaged performance of the model is 90.21% ± 4.46% VAF, and 95.14% ± 2.31% Correlation.

A remarkable fact is that the time delays (associated to the induced movement parameters, i.e. the angular position of the wrist, velocity, and acceleration) used in the model are the same among participants, indicating a widespread structure of the physiological pathways, which represents a new contribution that has not been reported yet. Then, the main contribution of this work is that, based on the model of the contralateral cortical response caused by wrist manipulation, the dot products (among the parameters of the induced movement) wrap up the nonlinear harmonics of the brain response at the same time delays, reflecting a general behavior among subjects.

Literature review

According to He and Yang [17], identification and modeling of neural systems involves both the linear and nonlinear approaches.

Yang et al. [18] found a correlation between the activity of the cortical region and motor tasks when measuring the cortical activity provoked by the flexion of the right-handed subjects wrist that modified the length of these muscles (isotonic contractions). In addition, Campfens et al. [19] found out that the addition of small continuous position perturbations to the wrist provoked significant coherence between the position perturbation signal and the EEG signal.

Specifically, to model the brain response due to wrist manipulation, Vlaar [20] and Tian [21] reported a truncated Volterra series-based model and a NARMAX-based model, explaining 46% and 69% of the brain activity variance, respectively. A comparison was performed by de Pont [22], arguing that those two methods did not include a parametric uncertainty level, which might improve the model evaluation method. To solve this, a Volterra series and Bayesian Inference-based model was proposed, however, the results were unstable due to the second-degree Volterra matrix being close to singular.

Another approach is the use of the Time-Delay Neural Network model (TDNN) [23]. It is noted that compared to the Volterra system with particle swarm optimization (PSO), the TDNN mathematical model has more adjustable parameters and takes less calculation time. As this approach is a neural network-based model, it has the black-box limitations, hence, limited physiological information can be extracted from this model.

Nozari et al. [11] used a strategy based on LLNF networks and auto-regressive approaches. In this model, a base is created with delayed inputs and outputs with time-lag limits. However, the number of parameters is not consistent throughout all the inputs and outputs, since it depends on the model structure and the time-lag limits, which is opposite to the present work where the number of parameters is fixed.

Gu et al. [12] identified a common model structure for healthy and young people by using the NARMAX method. They found that 13 model terms are similar for all the participants: five output linear terms, six input quadratic terms, one output quadratic term, and one constant. However, angular acceleration information is missing, and the different time lags do not provide specific explanation of the involved physiological pathways.

Materials and methods

Hardware and software

An agreement with the Laboratorio Nacional de Supercómputo del Sureste de México (LNSFootnote 3) was made to perform all stages of this work. The hardware used in this project has a maximal capacity of up to 50 GB for storage and 240 calculation cores simultaneously, up to 128 GB in RAM. Libraries of Python 3.8 were used. From Matlab 2022, the topoplot function of EEGLAB was used for plotting the cortical response topology.

Experiment and dataset

The Delft University of Technology made public a dataset of electroencephalographic recordings of 10 right-handed, young, and healthy subjects (six men and four women) [16] when applying force to keep position for wrist torque perturbations (active task) and doing nothing when a controlled movement is induced in their wrist (passive task). Only the data from the passive task is considered because it fits with the main goal of this work, that is to find the general brain response provoked by wrist manipulation. The experimental procedure was approved by the Human Research Ethics Committee of the TU Delft. All participants delivered written informed consent before the performance of experiments.

In the experiment (passive task), the subject sits in a chair, whose right hand is holding a joystick, and the arm is at rest and attached to a robotic joint (see left side of Fig. 1). Specifically, a wrist manipulation (realization) is the induced movement of a motor angular position. This movement profile is a multisine signal with frequencies 1, 3, 5, 7, 9, 11, 13, 15, 19, and 23 Hz and random phase. The first three frequencies have the highest amplitude, then it descends 20 dB at each frequency. A total of seven different realizations provoke movement of the right wrist. Three different realizations’ groups are shown at the top left of Fig. 1, where shadowed regions are removed to decrease transition effects. Finally, each realization is segmented into 210 periods of 1 s. These settings are established to allow similarity between realizations, but with a slight difference to avoid a learning process.

Fig. 1
figure 1

Image taken from benchmark “Cortical responses evoked by wrist joint manipulation” [20] and [24]. The permission from the main author was given

Full size image

Brain activity is recorded via EEG, with 126 channels and 2048 Hz sampling frequency. A fourth-order zero phase shift band-pass filter (cutoff frequencies at 1 Hz and 100 Hz) is applied to attenuate signals of no interest. A 50 Hz Notch filter is considered to attenuate the power line signal.

Methodology

From ten participants and seven realizations (inputs), preprocessing of the raw EEG recordings is performed before ICA. Then, the SNR values of all independent components are computed, and the signal with maximal SNR value and contralateral response (obtained from the topology maps) is chosen as the cortical response (a total of 70 cortical responses are computed). The found cortical responses (outputs) are then estimated via the nonlinear delay-differential embedding model. In this step, the best time-delays that mostly fit among all subjects and realizations are selected to compute the general response among all participants.

The suggested methodology is graphically abstracted in Fig. 2 and is divided into three sections: (A) wrist movement-related cortical response identification, (B) time delay terms estimation, and (C) model parameter approximation.

Fig. 2
figure 2

Methodology’s graphical abstract. A Wrist movement-related cortical response identification. B Best time delays terms estimation. C Model parameters approximation

Full size image

Cortical response identification

For each realization, a matrix \(\varvec{x}\) is built from the EEG recordings, and Independent Component Analysis (ICAFootnote 4) is applied as \(\varvec{x = As}\), where \(\varvec{A}\) is the weights matrix, and \(\varvec{s}\) the independent components matrix.

The Signal-to-Noise Ratio (SNRFootnote 5) is calculated for all independent components as in Algorithm 1. The cortical response is identified by selecting the component with the maximal SNR and its topology showing a contralateral response. The resulting component is averaged by periods to obtain the expected signal for model purposes.

A similar method as Vlaar (see [20]) was applied to calculate the SNR of each component (see Algorithm 1). The main difference with the method of the current work is that ICA is performed for each realization, i.e. for each subject there are seven ICA processes, instead of just one as in Vlaar’s method. This new approach accelerates the process of cortical response identification, and reduces computational cost.

Algorithm 1
figure b

SNR algorithm

Full size image

From Algorithm 1, T is the total number of periods, \(\varvec{\bar{IC}}_{i}(n)\) is the average of the i-th independent component with \(i \in [1,125]\), and r is the index of the independent component with maximal SNR at each realization. The cortical response is represented as \(\varvec{y}_{M}(n)\), which is provoked by the corresponding multisine signal \(\varvec{u}_{M}(n)\). From here, the realizations M will be implicit, hence the cortical response is then represented as \(\varvec{y}(n)\), and the Mth realization that provokes this response is shown as \(\varvec{u}(n)\).

Proposed model

The general expression to estimate of the cortical response caused by wrist manipulation is given by Eq. 1

$$\begin{aligned} \begin{array}{c} \varvec{\hat{y}}(n) =\varvec{\Phi }^{T}\varvec{\Gamma }\\ \varvec{\Phi }=[a_{1}, a_{2}, a_{3}, ..., a_{s}]\\ \varvec{\Gamma } = [\varvec{u}(n-\tau _{1}), \varvec{\dot{u}}(n-\tau _{2}), \varvec{\ddot{u}}(n-\tau _{3}), ..., \\ \varvec{\dot{u}}^{p-1}(n-\tau _{2})\varvec{\ddot{u}}(n-\tau _{3})] \end{array} \end{aligned}$$
(1)

where \(\varvec{\hat{y}}(n)\) is the modeled signal of the cortical response, \(\varvec{\Phi }\) is the weight vector including the contribution of each unitary basis (regressor) in the matrix \(\varvec{\Gamma }\), and \(\varvec{u}(n)\) is the corresponding realization (input) that provokes the brain response \(\varvec{y}(n)\). As \(\varvec{\Phi }\) is a filter to obtain the cortical response from the input, it represents the system’s impulse response. The solution to find the weight vector \(\varvec{\Phi }\) is expressed in Eq. 2 (a full description of the model is found in section “Theory”).

$$\begin{aligned} \begin{array}{c} \varvec{\Phi } = (\sum _{n=1}^{N}\varvec{\Gamma }(n,\tau _{k})^{2}) ^{-1}\sum _{n=1}^{k}\varvec{\Gamma }(n,\tau _{k})\varvec{y}^{T}(n) \end{array} \end{aligned}$$
(2)

Finding the best time delays of the proposed model

Rather than making any assumptions about physiological information when setting time delays, the particle’s exploration and exploitation technique, explained in the next paragraphs, determines the optimal parameters that reflect the distinctive relationship between the physiological component and signal behavior.

From Eq. 1, the proposed model considers common delayed signals of the realizations and its derivatives, i.e., lag terms are regarded for the angular position, velocity, and acceleration (\(\tau _{1}\), \(\tau _{2}\) and \(\tau _{3}\), respectively). Consider \(\varvec{Q} = [\tau _{1}, \tau _{2}, \tau _{3}]\) as a 3D particle.

The best time delays are found as follows: the exploration of the entire space with a certain step gives the first approximation of the best solution. Then, the exploitation of this particle (with a minor step) gives the final best particle. A high value of step means a low resolution. To avoid local minima, 10 possible solutions were obtained (because of 10 subjects), and the final solution considered the particle Q that gave the maximal correlation among all subjects’ impulse responses (see Algorithm 2).

Algorithm 2
figure c

Algorithm to find the best time delays of the model

Full size image

In terms of time delays computation, the computational cost is not considered because of the high complexity of the task (see Algorithm 2). For that reason, the use of supercomputer resources allows to identify those optimal time delays parameters that reflect the widespread brain response.

Approximation of the model parameters

From \(\varvec{\Phi }\) (weight vector) and \(\varvec{\Gamma }\) (unitary bases matrix) of Eq. 11, a new matrix \(\varvec{Z}\) is built as in Eq. 3, which gathers the algebraic elements that reflect the cortical response when summing them. Each input and each subject has its corresponding matrix \(\varvec{Z}\), hence a total of 70 matrices are built from the 10 subjects and seven realizations.

$$\begin{aligned} \varvec{Z}=\left[ \Phi [1] \left[ \begin{array}{c} \Gamma [1,1] \\ \Gamma [2,1] \\ ... \\ \Gamma [n,1] \\ \end{array} \right] ,..., \Phi [s] \left[ \begin{array}{c} \Gamma [1,s] \\ \Gamma [2,s] \\ ... \\ \Gamma [n,s] \\ \end{array} \right] \right] \end{aligned}$$
(3)

After building 70 matrices from \(\varvec{\Phi }\) (weight vector) and \(\varvec{\Gamma }\) (unitary bases matrix) and applying SVD (\(\varvec{Z=U\sum V}^{T}\)) to acquire its principal components, different patterns can be found by visually inspecting the matrix \(\varvec{V}\), which give information about the contribution of each unitary basis.

Model validation

To validate the proposed model, two setups are considered. In the first setup (labeled as 10sub), 10 subjects are considered to apply the leave-one-out cross-validation method. In the second setup (labeled as 11sub), all the EEG signals coming from 10 subjects are averaged and the resulting signals are established as belonging to an imaginary eleventh subject. This action aims to highlight the common activity among subjects caused by the movement induced in the wrist. The leave-one-out cross-validation method is applied as well. The performance metrics for this project are the VAF that calculates a variance ratio between the real and the modeled cortical responses elicited by wrist movement, the Pearson Correlation Coefficient (Corr) to see how similar signals are, and the RMSE that quantitatively calculate the estimation error (see Eqs. 4, 5, and 6, respectively).

$$\begin{aligned}{} & {} \begin{array}{c} VAF = \Big [1 - \frac{var(\varvec{y(t)} - \varvec{\hat{y}}(t))}{var(\varvec{y}(t))}\Big ]*100\\ \end{array} \end{aligned}$$
(4)
$$\begin{aligned}{} & {} \begin{array}{c} Corr = \frac{\sigma _{\varvec{y}(t),\varvec{\hat{y}}(t)}}{\sigma _{\varvec{y}(t)}\sigma _{\varvec{\hat{y}}(t)}}\\ \end{array} \end{aligned}$$
(5)
$$\begin{aligned}{} & {} \begin{array}{c} RMSE = \sqrt{\frac{1}{N}\sum _{1}^{N}[\varvec{y}(t) - \varvec{\hat{y}}(t)]^{2}}\\ \end{array} \end{aligned}$$
(6)

A polynomial filter was considered because there is intrinsic noise when averaging the estimated output signals. Therefore, the Savitzky–Golay filter (widely used to smooth biomedical signals [25]) is chosen, by considering a balance between smoothness and distortion of the signal, and effects of bias (over-smoothness) as well, when choosing the filter parameters [26, 27].

Theory

Consider \(\varvec{u}(n)\) as the input signal (realization), and \(\varvec{y}(n)\) as the output signal (component with maximal SNR and contralateral cortical response). Then, calculate the input’s numerical derivative and acceleration, filter and normalize them, as in Eq. 7.

$$\begin{aligned} \begin{array}{c} \varvec{\dot{u}}_{D}(n) = \frac{d}{dn}\varvec{u}(n), \ \ \ \ \ \ \ \ \varvec{\ddot{u}}_{D}(n) = \frac{d^{2}}{dn^{2}}\varvec{u}(n)\\ \varvec{\dot{u}}_{SG}(n) = SavGol(\varvec{\dot{u}}_{D}(n)) \\ \varvec{\ddot{u}}_{SG}(n) = SavGol(\varvec{\ddot{u}}_{D}(n))\\ \varvec{\dot{u}}(n) = \frac{\varvec{\dot{u}}_{SG}(n)}{std[\varvec{\dot{u}}_{SG} (n)]}, \ \ \ \varvec{\ddot{u}}(n) = \frac{\varvec{\ddot{u}}_{SG}(n)}{std[\varvec{\ddot{u}}_{SG}(n)]} \\ \end{array} \end{aligned}$$
(7)

where SavGol is the Savitzky–Golay filter used to reduce intrinsic noise caused by numerical processes [25, 26], and \(std[\varvec{x}(n)]\) is calculated as in Eq. 8.

$$\begin{aligned} \begin{array}{c} std[\varvec{x}(n)] = \sqrt{\frac{1}{N-1}\sum _{n=1}^{N} \Big (\varvec{x}(n) - \frac{1}{N}\sum _{n=1}^{N}\varvec{x}(n)\Big )^{2}}.\\ \end{array} \end{aligned}$$
(8)

The unitary bases are built as in Eq. 9.

$$\begin{aligned} \begin{array}{c} \varvec{\Gamma }(n, \tau _{k}) = [\varvec{u}(n-\tau _{1}),\varvec{\dot{u}}(n-\tau _{2}), \\ \varvec{\ddot{u}}(n-\tau _{3}),\varvec{u}^{2}(n-\tau _{1}), \varvec{u}^{3}(n-\tau _{1}) ,...,\\ \varvec{u}^{p-1}(n-\tau _{1})\varvec{\dot{u}}(n-\tau _{2}), \varvec{u}^{p-1}(n-\tau _{1})\varvec{\ddot{u}}(n-\tau _{3}), \\ \varvec{\dot{u}}^{p-1}(n-\tau _{2})\varvec{\ddot{u}}(n-\tau _{3})]\\ \end{array} \end{aligned}$$
(9)

\(\tau _{1}\), \(\tau _{2}\) and \(\tau _{3}\) are the time delays of the angular position, velocity, and acceleration, respectively, and p is the model degree, which is established via a trade-off between the performance of 15 different model’s degree and its computational cost (number of parameters to calculate).

The estimated output signal is calculated as in Eqs. 10 and 11.

$$\begin{aligned}{} & {} \begin{array}{c} \varvec{\hat{y}}(n,\tau _{k}) = a_{1}\varvec{u}(n-\tau _{1}) + a_{2}\varvec{\dot{u}}(n-\tau _{2}) +... \\ + a_{s-1}\varvec{u}^{p-1}(n-\tau _{1})\varvec{\ddot{u}}(n-\tau _{3}) + \\ a_{s}\varvec{\dot{u}}^{p-1}(n-\tau _{2})\varvec{\ddot{u}}(n-\tau _{3}) \end{array} \end{aligned}$$
(10)
$$\begin{aligned}{} & {} \begin{array}{c} \varvec{\hat{y}}(n, \tau _{k}) =\varvec{\Phi }^{T}\varvec{\Gamma }(n, \tau _{k}) \end{array} \end{aligned}$$
(11)

where \(\varvec{\Phi }=[a_{1}, a_{2},..., a_{s-1}, a_{s}]\) is the weight vector that contains the contribution of each unitary basis located in the column space of the matrix \(\varvec{\Gamma }(n, \tau _{k})\). Finally, calculate \(\varvec{\Phi }\) as in Eq. 14.

$$\begin{aligned}{} & {} \begin{array}{ll} \varvec{\tilde{y}}(n,\tau _{k}) = \varvec{y}(n) - \varvec{\hat{y}}(n, \tau _{k}) \\ \end{array} \end{aligned}$$
(12)
$$\begin{aligned}{} & {} \begin{array}{ll} \varvec{f}(n,\tau _{k},\Phi ) = \frac{1}{2N}\sum _{n=1}^{N}\big (\varvec{y}(n) - \varvec{\Phi }^{T}\varvec{\Gamma }(n,\tau _{k})\big )^{2}\\ \end{array} \end{aligned}$$
(13)

Performing \(\partial {\varvec{f}(n, \tau _{k}, \varvec{\Phi })}/\partial {\varvec{\Phi }} = \textbf{0}\), the solution of the parameters vector \(\varvec{\Phi }\) is found in Eq. 14

$$\begin{aligned} \begin{array}{ll} \varvec{\Gamma }(n,\tau _{k})^{2} = \varvec{\Gamma }(n,\tau _{k})\varvec{\Gamma }^{T}(n,\tau _{k})\\ \varvec{\Phi } = \big (\sum _{n=1}^{N}\varvec{\Gamma }(n,\tau _{k})^{2}\big ) ^{-1}\sum _{n=1}^{k}\varvec{\Gamma }(n,\tau _{k})\varvec{y}^{T}(n) \end{array} \end{aligned}$$
(14)

Results

Cortical response topology

After applying ICA, the cortical responses caused by the seven realizations were identified by choosing the component with maximal SNR and plotting its topology showing a contralateral activity. Figure 3 shows the cortical responses for ten subjects when the first realization M0 was applied.

Fig. 3
figure 3

Topology of the cortical response for ten subjects and first realization M0

Full size image

From Fig. 3, the brain activity is similar for most of the subjects and is located on the contralateral motor cortex (central left area). It must be remembered that the movement is induced in the right hand. These results are similar to those obtained by Vlaar in [20], with the advantage of reduced computational complexity and, therefore, a faster method to obtain the cortical response.

Establishing the model degree

In terms of the number of parameters, whether the 13th-degree model is chosen (considering 273 parameters), some information of the algebraic bases might be missed caused by the truncation of the eigenvalues matrix \(\varvec{\sum }\) (whose dimension is 256, the same as the sampling frequency). Therefore, the closest model is selected, that is, the 12th-degree model with 234 parameters to compute (see right side of Fig. 4).

Fig. 4
figure 4

Left: VAF Vs Model degree. Right: Number of parameters according to the model degree, realization M=1

Full size image

Moreover, in terms of performance (see left side of Fig. 4), when the model degree is 10, the VAF is between 70% and 80% with 165 parameters to find. Considering a 12th-degree model, there are 234 parameters to calculate, achieving VAF around of 90%. When the model degree is 13, there are 273 parameters to calculate, and VAF increases to almost 100%, however, as previously indicated, some information about the algebraic bases is missed.

In this sense, regarding the 12th-degree model, a lower degree model will have lower performance, and a higher one, reaching VAF 100%, could possibly have overfitting effects, which is not desired when computing a general model because it limits the addition of more subjects to the experiment.

Therefore, for mathematical modeling purposes in the present work, the 12th-degree model allows to build an algebraic basis to explain the cortical response elicited by wrist manipulation.

Best time delays of the model

Exploration This process begins with a low resolution (\(step = 6\) units, equivalent to 23.4 ms, considering a 256 Hz sampling frequency) through the entire space delimited from 1 to 128 (\(\tau _{b1} = \tau _{b2} = \tau _{b3} = 1\), and \(\tau _{e1} = \tau _{e2} = \tau _{e3} = 128\)), where 128 samples represent a time delay of 500 ms. The best solution was found as \(\varvec{Q}=[72,12,18]\) samples.

Exploitation After exploration, the exploitation of the particle is performed considering the highest resolution (\(step = 1\)) around the latter best particle. The new limits for each axis are determined by adding/subtracting the previous step value (\(step = 6\)) from the previous best delay particle values.

The best time delays of the model for the angular position, velocity, and acceleration were 281 ms, 47 ms, and 70 ms, respectively, corresponding to \(\tau _{1}=72\), \(\tau _{2}=12\), and \(\tau _{3}=18\) samples, respectively, with a 256 Hz sampling frequency.

Model validation

VAF and Correlation metrics are important for this research because they reflect the model performance by focusing on different aspects. While the VAF metric highlights the variance ratio between the error signal and the cortical response (100% VAF indicates identical signals, therefore, high VAF values reflect a better model performance), the Pearson Correlation shows the similarity between the identified cortical response and its estimation. Results of the cross-validation process for both setups (10sub and 11sub) are shown in Table 1.

Table 1 Model validation results
Full size table

RMSE measure displayed unexpected results when comparing the cross-validation processes in M6 and M0 with respect to the 10sub setup (see Table 1). While VAF and Correlation metrics indicate an improvement in model performance (90.2359% to 91.0885% VAF and 0.9507 to 0.9571 Correlation), the RMSE metric indicates a worse performance (0.0926 to 0.1096). Regarding the 11sub setup, the same is true. This can be explained by the fact that even little changes in the signal have a big impact on the model’s performance when using the RMSE metric. Hence, it is advised to use VAF and correlation to evaluate model performance of related works.

Finally, from Table 1, the best averaged performance of the model is 90.21% ± 4.46% VAF, and 95.14% ± 2.31% Correlation when considering the 11sub setup, indicating that highlighting the average behavior of the cortical response elicited by wrist manipulation emphasizes the intrinsic relationships among subjects and improves the model performance.

Discussion

In this section, the modeled cortical response, the widespread and predominant operation mode, and the corresponding algebraic basis are explained. A physiological interpretation of the best time delays of the model based on previous works is given at the end of this section.

Analysis of the modeled cortical response

The upper plot of Fig. 5 shows in blue the expected cortical response for subject S_01 when perturbation M1 is applied, and the orange signal represents the estimated cortical response calculated from the first cross-validation setup 10sub. Intrinsic noise is represented by some pikes centered at samples 10, 50, and 120, approximately (shadowed in red circles).

Fig. 5
figure 5

Expected (blue) and estimated cortical response (orange) correspond to the first cross-validation test without the SavGol filter (upper plot) and after the filter (lower plot) when realization M1 is applied

Full size image

After the SavGol filter (best parameters at \(window lenght=7\) and \(polyorder=3\)), these pikes are attenuated and smoothly decreased, hence, the model performance improved by increasing the correlation in percentage (from 93.74% to 94.81%), the VAF values (from 87.8336 to 89.7890%), and decreasing the RMSE (from 0.1259 to 0.1153).

General predominant operation mode and algebraic basis

The left side of Fig. 6 shows the average of the impulse response signals \(\varvec{\Phi }\) (weight vector) for all subjects in the cross-validation setup 10sub for realization M0. Visually comparing the graphs, it is possible to see that all subjects present a similar impulse response for each realization, which means that the intrinsic behavior is the same among all of them, but slightly different among realizations.

Fig. 6
figure 6

Predominant operation mode or impulse response (left), and Correlation Profile (right) for realization M0

Full size image

Table 2 shows that, of the 234 parameters in the proposed 12th-degree model, 10 only include even powers and contribute to a quarter of the general brain response. Additionally, nine of those model terms (dot product between various powers of the input’s angular velocity, \(\varvec{u}_{\tau _{2}}\), and acceleration, \(\varvec{u}_{\tau _{3}}\)) pike at 10 Hz and 12 Hz, and then descend in amplitude at 22 Hz, 24 Hz, 32 Hz, and 34 Hz. This suggests that those model terms are physiologically related to the intermodulation of the excited frequencies in the input, since only odd frequencies are contained in the perturbation, highlighting the dominance of the perturbation’s even powers (which is consistent with earlier findings by Vlaar [28]) and the impact of movement acceleration on the common cortical response due to \(\varvec{u}_{\tau _{3}}\).

Table 2 Cumulative contribution of the induced movement terms with more energy (unitary bases) to the widespread brain response
Full size table

In terms of correlation profiles, they illustrate the model’s performance as the model terms are added one at a time, peaking after all terms have been included. The correlation profiles for the cross-validation setup 10sub for realization M0 (see the right side of Fig. 6) reflect the algebraic bases that produce each subject’s cerebral response, indicating a common tendency of the model performance among individuals, and, therefore, highlighting a common behavior of brain responses among participants.

Sometimes the sign is inverted, but at the end, the correlation between the expected and estimated output has a positive sign and maximal value. A similar case appears in [12], but no explanation is given. According to our results, the change of sign might be provoked by the realization’s odd and even powers, i.e., since having only input’s odd frequencies, the sign of the estimated signal might change when summing up odd or even components.

Finally, a broader approach is presented in the current paper, which considers not only the angular position [20] and velocity [21], but the acceleration as well, allowing a richer framework in which the respective time delays reflect specific physiological pathways.

Contribution of the realization’ even and odds powers to the cortical response

After building 70 matrices from \(\varvec{\Phi }\) (weight vector) and \(\varvec{\Gamma }\) (unitary bases matrix) (see subsection 2.3.4) and applying SVD (\(\varvec{Z=U\sum V}^{T}\)) to acquire its principal components, two patterns appeared corresponding to the even and odd powers of the realization. The average values of these components in the 70 matrices (considering all subjects and all realizations) were calculated, so the same behavior appeared as shown in Fig. 7.

Fig. 7
figure 7

Averaged \(\varvec{V}^{T}\) for input’s even powers (left). Averaged \(\varvec{V}^{T}\) for input’s odd powers (right)

Full size image

From the SVD process, the matrix \(\varvec{V}^{T}\) gave information about the most relevant unitary bases that relatively contributed to the cortical response (see Table 3), since that information is found in different principal components.

Table 3 Average values of the 70 matrices (all subjects and all realizations) of the most relevant even and odd unitary bases that relatively contribute to the cortical response
Full size table

Making a comparison between Tables 2 and 3, the 11.3819% relative contribution of \(\varvec{u}^{10}_{\tau _{1}}\) can be interpreted as 19.0871 minus 16.654 (see Table 2), resulting a 2.4331% contribution to the general cortical response. In this sense, contribution of the even and odd powers of the input (\(\varvec{u}_{\tau _{1}}\)) to the cortical response are 6.2052% and 5.1411%, respectively. These results show that the even behavior is dominant, which is consistent with the previous work of Vlaar [20].

Yang et al. [18] found out that there exists a nonlinear coupling for both integer and non-integer harmonics of the realization. Additionally, closed-loop connections like the NARX model can estimate more wrist dynamics information [12] than open-loop connections such the Volterra model [20]. However, the brain response has more harmonics which are not explained and its source of generation is not well understood. Our results indicates that the input’s even and odd powers produce such harmonics, i.e. the different powers reproduce the feedback effect of the closed-loop connection.

These findings suggest the patterns order in terms of the general contribution percentage to the cortical response in both the even terms (\(\varvec{u}^{10}_{\tau _{1}}\), \(\varvec{u}^{8}_{\tau _{1}}\), and \(\varvec{u}^{12}_{\tau _{1}}\)) and the odd terms (\(\varvec{u}^{9}_{\tau _{1}}\), \(\varvec{u}^{7}_{\tau _{1}}\), and \(\varvec{u}^{11}_{\tau _{1}}\)) could be used as a potential footprint to characterize the common behavior of the healthy and young brain response provoked by specific wrist manipulation. Therefore, the arrangement of each pattern discovered in the even and odd powers can also be suggested as an additional potential biomarker.

Physiological interpretation

The proprioceptive system is a key element to convert mechanical events (sensed via mechanoreceptors, such as muscle spindles or Golgi tendon organs [28]) into neural signals [29], where EEG recordings of brain activity reflect the contribution of many neurons at a macroscopic level. After applying ICA to the filtered EEG recordings, the component with maximal SNR (founded at the contralateral side of the perturbation) represents the activation source of the brain response caused by the external wrist manipulation.

Furthermore, regarding both \(\tau _{2}\) and \(\tau _{3}\) time delays parameters of the proposed model (47 ms and 70 ms, respectively), information of the perturbation’s angular velocity and acceleration is sent to the brain via the groups Ia-afferent and 1b-afferent, respectively, and processed by a partially modulated left lateral model [30], reflecting the change rate of muscle stretching and explosive behaviors provoked by the force applied to the wrist [18], respectively. Hideaki Onishi et al. [31] reported similar time courses (36.2 ± 8.2 ms, and 86.1 ± 12.1 ms) for passive finger movements.

The wrist angular position is the slowest signal sent to the brain (281ms time delay) via the group II-afferent, suggesting a long latency response involving more brain areas and processed by a fully modulated left lateral model [30].

In summary, the obtained results suggest that the intermodulation products arising from delayed signals of the velocity and acceleration of the wrist perturbation (and controlled by the groups Ia-afferent and 1b-afferent, respectively) play a significant role in the nonlinear behavior of the general cortical response, highlighting the dominance of the unexcited even frequencies associated with the even powers of the model terms.

Therefore, those model terms that exhibit the largest amplitudes of the general impulse response when considering all subjects and all realizations are suggested as possible biomarkers to be used in future research on the brain activity of individuals with upper limbs motor impairments.

Conclusions

In this work, a better understanding of the relationship between a mechanical wrist perturbation and brain response in healthy people was achieved, allowing the discovery of potential biomarkers that reflected this link. In the future, with new EEG recordings of patients with motor impairments, the methodology proposed in this work can be replicated to evaluate the behavior of those potential biomarkers for people with some motor deficiencies.

In terms of complexity reduction, a simpler model (from the 5th-degree model and higher where specific details are not important, but the general brain response to specific wrist manipulation), and the cortical response’s reconstruction from even frequencies would help to reduce the number of model terms to compute.

Furthermore, the identification of the cortical response might be simplified by using fewer EEG electrodes. Moreover, in terms of time delays computation, reducing time consuming using optimization methods, for example, metaheuristics, can be another research topic.

Therefore, regarding those elements could help to develop useful interfaces as medical applications to identify potential biomarkers that reflect the patient’s current state in real time. This would make it easier for therapist experts to evaluate and compare the effectiveness of rehabilitation sessions for people with upper limbs motor impairments.