Learning Cortico-Muscular Dependence through Orthonormal Decomposition of Density Ratios

We thank the reviewer for the comments. Please find below our replies to the concerns/questions.

1. Limited baselines.

We have added EEG-Conformer (available on GitHub) and Deep4 (from the braindecode repository on GitHub) in the attached pdf. All baselines were implemented following the official codes. Results indicate that FMCA-T outperforms the additional baselines on almost all tasks. Deep4 performs slightly better on cross-subject 11-movement classification task. We will update Table 1 in the revised manuscript.

2. Computational resources for training and inference.

Experiments on the simulated sinusoidal dataset and the additional simulated EEG and EMG dataset were conducted on an NVIDIA GeForce RTX 3090. Experiments on the experimental EEG and EMG dataset were conducted on an NVIDIA GeForce A5000. All training and inference can be conducted on a single GPU.

3. Limited dataset size.

We agree with the reviewer about needing a larger, diverse dataset for better validation. Although unavailable now, we hope our method will encourage large-scale datasets, which we plan to explore further.

4. Generalization to other modalities.

We thank the reviewer for proposing the possible extension to broader neural data modalities. This is an interesting question for further investigation. We expect our method can still apply to learning the dependence between other modality data, e.g., between EMG and kinematics or between EEG and visual data. We will continue studying such scenarios in future work.

5. Ethics review.

The public EEG and EMG dataset was reviewed and approved by the Institutional Review Board at Korea University (1040548-KU-IRB-17-181-A-2).

We thank the reviewer for the instructive comments. Please find below our replies to the concerns/questions.

1. Mathematical and algorithmic complexity. We thank the reviewer for the suggestion. We will add more explanations of our methodology to the supplementary.

2. Interpretation of results and physiological implications of the identified activations. In terms of interpretation, we found that the frontal central areas (FC) are most activated in the spatial-level dependence maps for most subjects. Since we used eigenfunctions that decompose the density ratio for classification, this implies that these FC areas contribute significantly to movement classification. It is reasonable for FC1 to show the most activation in the dependence map. A related previous study also found movement-related cortical potential changes on FC1-FC2 and C2 electrodes [1].

We have additionally added a frequency analysis using event-related desynchronization (ERD) in the beta band. In the general response, we addressed the differences and consistency between frequency topologies and our dependence ratio map.

3. Scalability. We appreciate the reviewer's suggestion. Given the limited time frame and the limited availability of paired EEG-EMG datasets, it is difficult to show scalability on new real-world datasets. Instead, we divided the dataset into smaller sub-datasets and ran experiments with increased dataset sizes to demonstrate how our framework scales up. Fig. 4(a)~4(c) in the attached letter show that: (a) the convergence speed of training errors is almost the same, with smaller datasets converging slightly faster; (b) downstream classification accuracies steadily increase when using more data; (c) dependence decreases when using more data.

We observed the expected trend in dependence scores - as more data is used, it implies greater uncertainty in the data.

The computational complexity was not particularly heavy, and we did not see it as an obstacle. We trained the model on an A5000 GPU and did not encounter significant difficulties.

4. Potential applications. Corticomuscular coupling is promising for quantitively measuring movement disorders. This application has been validated in stroke populations [2, 3] and in Parkinson’s disease [4]. Therefore, the proposed dependence measure can be used as an effective and clinically relevant neural marker to evaluate movement performance of stroke patients and to indicate Parkinson’s disease pathology. Moreover, as demonstrated in the manuscript, EEG eigenfunctions show excellent performance in movement classifications. Therefore, the proposed method also has great potential in improving brain-computer interfaces. We will add this discussion in the revised manuscript.

[1] Spring J N, Place N, Borrani F, et al. Movement-related cortical potential amplitude reduction after cycling exercise relates to the extent of neuromuscular fatigue[J]. Frontiers in human neuroscience, 2016, 10: 257.

[2] Chen X, Xie P, Zhang Y, et al. Abnormal functional corticomuscular coupling after stroke[J]. NeuroImage: Clinical, 2018, 19: 147-159.

[3] Liu J, Wang J, Tan G, et al. Correlation evaluation of functional corticomuscular coupling with abnormal muscle synergy after stroke[J]. IEEE Transactions on Biomedical Engineering, 2021, 68(11): 3261-3272.

[4] Zokaei N, Quinn A J, Hu M T, et al. Reduced cortico-muscular beta coupling in Parkinson’s disease predicts motor impairment[J]. Brain communications, 2021, 3(3): fcab179.

We thank the reviewer for the detailed insightful feedback.

1.1 Novelty over CCA. We acknowledge the reviewer's observation that the final cost is the trace of a normalized canonical correlation matrix between two multivariate neural networks.

But we emphasize the link between cost optimization and joint density ratio $\frac{p(X,Y)}{p(X)p(Y)} = \sum_{k=1}^K \lambda_k \widehat{\phi_k} (X) \widehat{\psi_k} (Y)$, which the reviewer might have misunderstood.

Vanilla CCA uses linear models to maximize a correlation coefficient. After KICA was introduced (i.e., Kernelized CCA), the problem became minimizing matrix costs using universal approximators, with costs including log determinant, Frobenius norm, and matrix trace. Our trace cost is closer to HSIC. We have thoroughly compared our method with KICA and HSIC as baselines, both of which use kernels while we use neural nets.

Only recently, studies proved that optimizing matrix costs is equivalent to factorizing joint density ratios, including FMCA and Gaussian universal features. An informal proof: if $\widehat{\phi}$ and $\widehat{\psi}$ are two orthonormal sets, optimizing the norm of $\mathbb{E}[\mathbf{\widehat{\phi}}(X)\mathbf{\widehat{\psi}^\intercal}(Y)] = \int \mathbf{\widehat{\phi}}(X)\mathbf{\widehat{\psi}}^\intercal(Y) p(X,Y)dXdY$ decomposes the joint density ratio $\frac{p(X, Y)}{p(X)p(Y)}$ and finds its eigenfunctions.

This paper further shows that all such costs, including log determinant, Frobenius norm, and trace, are fundamentally equivalent, differing only in the convex functions applied to eigenvalues. We consider this a major novelty.

1.2 More discussions about novelties. Clarifying the question on $\mathbf{z}$, it is the variable making $\mathbf{X}$ and $\mathbf{Y}$ conditionally independent, such as movement types, not the network's latent features that approximate eigenfunctions. Lemma~1 explains why eigenfunction can be used for downstream tasks. Another novelty is analyzing dependence within the network for temporal and spatial activation.

1.3 Normalization. Eigenvalues of $\frac{p(X,Y)}{p(X)p(Y)}$ are well-regularized and don't need normalization: (1) All eigenvalues are positive and bounded by 1; (2) $\lambda_1=1$ regardless of dependence; (3) independence iif there exists only one non-zero eigenvalue. This is due to the implicit normalization dividing $p(X,Y)$ by $p(X)p(Y)$, corresponding to the cost's normalization $\mathbf{R}_F^{-\frac{1}{2}} \phi$ and $\mathbf{R}_G^{-\frac{1}{2}} \psi$.

The normalization is costs also enforces orthonormality constraints on features, so outputs won't be constants (trivial solution)

2.1 Simulated dataset. We have added additional experiments with simulated dataset, shown in Fig. 3 in the general response.

We simulated EEG and EMG signals for left/right motor and sensory activations in 20 subjects using EEGSourceSim. Motor sources were used to simulate the corresponding EMG signals. We trained the networks on paired EEG-EMG samples from 16 subjects' data as training samples, and visualized the spatial-level dependence maps for the remaining 4 subjects as test samples. Then we plotted ground truth brain activations calculated from motor ROI and forward matrices.

The spatial-level dependence is highly similar to ground truth activations, indicating the learned ratio captured real brain activations.

2.2 Dependence shows major activations around FC areas. We have added a frequency analysis using event-related desynchronization (ERD) in the beta band (Fig. 2) and complete spatial dependence maps of our method (Fig. 3), both for Subject 3.

The frequency topology shows that left-central areas (e.g., C3) are commonly activated across sessions and movements while frontal-central (FC) and central-parietal areas have unique patterns for each movement. In comparison, our dependence maps show the most activated areas around FC areas.

We argue that the difference in activation occurs because we are quantifying nonlinear dependence, not correlations. It has been experimentally proven that the eigenfunctions of this dependence yield high accuracies in downstream classification tasks, which means that the activation maps should show the areas that make each movement most distinguishable. Thus the frontal central and central parietal areas are activated in dependence maps. In frequency topologies, C3 is activated for all movements, but it's unable to be used as an identifier for individual movements. This causes the differences between our maps and frequency maps. More details can be found in the general response.

The public dataset used has been validated before publication. No activation or deactivation was found on peripheral or frontal electrodes in the frequency analysis-based topographies.

2.3 Electrode arrangement and normalization steps. The electrode arrangement is reliable and has been used in commercial equipment. In EEG analysis, if the target electrode is near the reference, a weighted average filter will be employed to enhance the signal. Our study instead had no clear target before analysis. Therefore, we applied an average reference across electrodes, which is commonly used in EEG analysis. During preprocessing, all EEG channels were normalized to the first sample's maximum amplitude. Trials exceeding an absolute amplitude of 5 were discarded. No channel-wise normalization was used that could affect the variance.

3.1 Choice of nonlinearity differs from EEGNet. In EEGNet paper, the authors found that nonlinear activations did not improve performance, which is why they chose linear functions, thus making it a matter of preference.

Although linearity may help prevent overfitting in cases of long signal durations, our experiments showed no significant overfitting. Using ReLU or Sigmoid may provide a more stable approximation when computing cross-layer ratios, as the ratios should be non-negative and bounded.

We thank the reviewer for the comments.

1. Classification baselines could be stronger. We have now added EEG-Conformer (available on GitHub) and Deep4 (from the braindecode repository on GitHub) as new baselines in the general response letter. The results show that FMCA-T surpasses the added baselines in almost all tasks, with Deep4 slightly ahead in cross-subject 11-movement classification.

2. Independence criterion with the eigenvalues. The reviewer correctly pointed out our mistake. A more precise argument would be: two variables are independent (i.e., $p(X, Y) = p(X)p(Y)$) if and only if the spectrum has a single positive eigenvalue $\lambda = 1$.

Define the linear operator $\psi := \mathcal{L}\phi = \int \frac{p(X,Y)}{p(X)p(Y)} \phi(X) \, d\mathbb{P}_X$ with $\phi$ and $\psi$ as functions of $X$ and $Y$. It can be confirmed that $\phi=1$ yields $\psi=1$, making it an eigenfunction of the operator with eigenvalue $\lambda=1$. Two variables are independent if and only if the only positive eigenvalue is this trivial one. All eigenvalues of the operator are non-negative and bounded by $1$. More than one positive eigenvalue implies dependence. We will correct the error.

3. Applicability and limitations of conditional independence assumption.

a. We demonstrated the system's ability to capture discrete and easily distinguishable factors like movement types, participants, and sessions. But it struggles with finer-grained factors, such as sub-movements that include arm-reaching in six directions, hand-grasping three objects, and wrist-twisting with two motions.

One possible explanation is that non-invasive techniques like EEG face more challenges to extract finer movement information compared to invasive techniques. Consequently, the ratio $\rho(X,z)$ becomes trivial, and finer information may not be easily extracted from the dependent components between EEG and EMG.

b. Future work is needed for scenarios with a continuous conditioned variable $\mathbf{z}$. In Lemma 1, we proposed the expression $\frac{p(X,Y)}{p(X)p(Y)} = \sum_{z\in \mathcal{Z}} p(z) \rho(X, z)\rho(Y, z) dz = \sum_{k}\lambda_k \phi_k(X)\psi_k(Y)$, which applies when $\mathbf{z}$ is discrete. This may not account for continuous variables such as the force of muscle contractions and kinematics (joint angle).

4. Density of texts. We appreciate the reviewer's feedback and will improve the formatting.