%%The option is available for: sn-basic.bst, sn-chicago.bst% % \documentclass[pdflatex,sn-nature]{sn-jnl}% Style for submissions to Nature Portfolio journals % \documentclass[pdflatex,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style \documentclass[pdflatex,sn-mathphys-num]{sn-jnl}% Math and Physical Sciences Numbered Reference Style %%\documentclass[pdflatex,sn-mathphys-ay]{sn-jnl}% Math and Physical Sciences Author Year Reference Style %%\documentclass[pdflatex,sn-aps]{sn-jnl}% American Physical Society (APS) Reference Style % \documentclass[pdflatex,sn-vancouver-num]{sn-jnl}% Vancouver Numbered Reference Style % \documentclass[pdflatex,sn-vancouver-ay]{sn-jnl}% Vancouver Author Year Reference Style % \documentclass[pdflatex,sn-apa]{sn-jnl}% APA Reference Style % \documentclass[pdflatex,sn-chicago]{sn-jnl}% Chicago-based Humanities Reference Style \usepackage{lineno} \usepackage{graphicx}% \usepackage{multirow}% \usepackage{amsmath,amssymb,amsfonts}% \usepackage{amsthm}% \usepackage{mathrsfs}% \usepackage[title]{appendix}% \usepackage{xcolor}% \usepackage{textcomp}% \usepackage{manyfoot}% \usepackage{booktabs}% \usepackage{algorithm}% \usepackage{algorithmicx}% \usepackage{algpseudocode}% \usepackage{listings}% \usepackage{afterpage} \usepackage{lipsum} \usepackage{placeins} \setcitestyle{super,comma,sort&compress,open={},close={}} % \bibliographystyle{naturemag} \usepackage{xr} \externaldocument{sn-article} \theoremstyle{thmstyleone}% \newtheorem{theorem}{Theorem}% meant for continuous numbers %%\newtheorem{theorem}{Theorem}[section]% meant for sectionwise numbers %% optional argument [theorem] produces theorem numbering sequence instead of independent numbers for Proposition \newtheorem{proposition}[theorem]{Proposition}% %%\newtheorem{proposition}{Proposition}% to get separate numbers for theorem and proposition etc. \theoremstyle{thmstyletwo}% \newtheorem{example}{Example}% \newtheorem{remark}{Remark}% \theoremstyle{thmstylethree}% \newtheorem{definition}{Definition}% \raggedbottom %%\unnumbered% uncomment this for unnumbered level heads \begin{document} \linenumbers \section*{Supplementary Notes}\label{SuppNote} \setcounter{figure}{0} \renewcommand{\thefigure}{S\arabic{figure}} \renewcommand{\theHfigure}{Supplement.\thefigure} % Create unique internal link anchor \begin{figure}[htbp!] \centering \includegraphics[width=1.0\textwidth]{SupFig/SupFigR1.png} \caption{\textbf{RoA and XCC are nonlinearly correlated, with dataset-dependent agreement that persists after duplicate removal (subset spike approach).} \textbf{a} Relationship between the spike train matching metric (RoA) and the waveform tracking metric (XCC), incorporating MU-specific random effects. A log-transformed mixed-effects model provided the best fit in both datasets, with RoA exerting a significant population-level effect on XCC ($p<.001$). Left: 1DoF dataset; Right: GRASP dataset. \textbf{b} Distribution of XCC values at each controlled RoA level, with horizontal lines indicating commonly used waveform-tracking thresholds. Data points above a given threshold indicate agreement between waveform tracking and spike train matching; data points below indicate disagreement. XCC variance is substantially higher at low RoA levels in both datasets. \textbf{c} Mean absolute deviation (MAD) of XCC at each RoA level, quantifying the decrease in XCC variability as RoA increases. A significant negative linear relationship between MAD and RoA was found in both datasets ($p<.001$). \textbf{d} Edge-level agreement between spike train matching and waveform tracking at each individual RoA level, computed as the proportion of MU pairs between the reference MU and the artificial MU for which both methods agree, given the applied waveform-tracking threshold. Agreement varies substantially across RoA levels and across threshold values in both datasets. Left: 1DoF dataset; Right: GRASP dataset. \textbf{e} Edge-level disagreement recomputed after excluding MU pairs with RoA $> 30\%$ (duplicate MUs \cite{holobar2010experimental}, shown as gray dots), retaining only pairs considered distinct by spike train matching. Disagreement between the two methods persists across all commonly used waveform-tracking thresholds, even after duplicate removal. ***: $p<.001$; **: $p<.01$; *: $p<.05$.} \label{supfigR1} \end{figure} \subsection*{Supplementary Note 1: Model selection and RoA-XCC correlation for both surrogate spike approach and subset spike approach}\label{SuppNote1} We examined four models (linear, quadratic, log-transformed and exponential, see Method \ref{Method-S5-SS1}) on two datasets (1DoF and GRASP) under two approaches of artificial MU generation (surrogate spike approach and subset spike approach; see Method \ref{Method-S5-SS1}). \par For the surrogate spike approach, the BIC value we found: $-1.90\times10^5$ (1DoF-linear), $-7.74\times10^4$ (GRASP-linear), $-3.12\times10^5$ (1DoF-quadratic), $-1.25\times10^5$ (GRASP-quadratic), $-3.13\times10^5$ (1DoF-log transformed), $-1.22\times10^5$ (GRASP-log transformed), $-3.75\times10^5$ (1DoF-exponential), $-1.52\times10^5$ (GRASP-exponential). \par We therefore select the NLME model with an exponential transformation as our final model for the correlation study. Moreover, we found a significant population-level effect of RoA on XCC in both datasets, described by $\mathrm{XCC}=1-\exp(-6.72\times\mathrm{RoA})$, $p<.001$ (1DoF) and $\mathrm{XCC}=1-\exp(-7.00\times\mathrm{RoA})$, $p<.001$(GRASP) \par For the subset spike approach, the BIC value we found: $-2.87\times10^5$ (1DoF-linear), $-1.21\times10^5$ (GRASP-linear), $-4.332\times10^5$ (1DoF-quadratic), $-1.78\times10^5$ (GRASP-quadratic), $-4.70\times10^5$ (1DoF-log transformed), $-1.87\times10^5$ (GRASP-log transformed), $-2.73\times10^5$ (1DoF-exponential), $-1.16\times10^5$ (GRASP-exponential). We therefore select the linear mixed-effects model with a log transformation as our final model for the correlation study. Moreover, we found a significant population-level effect of RoA on XCC in both datasets, described by $\mathrm{XCC}=1.00\ln(\mathrm{RoA})+0.58$, $p<.001$ (1DoF) and $\mathrm{XCC}=0.93\ln(\mathrm{RoA})+0.61$, $p<.001$(GRASP) % ----------------------------------------------------------------------- \subsection*{Supplementary Note 2: P values and effect sizes for dataset-dependent RoA-XCC differences across RoA levels}\label{SuppNote2} \hspace{15pt} In the main text (Section \ref{Result-S1-SS1}), we reported that the RoA-XCC correspondence is dataset-dependent, with significant differences between the 1DoF and GRASP datasets. Here we report the full $p$ values and Cohen's $r$ effect sizes at each of the nine RoA levels for both the surrogate spike and subset spike approaches. \par For the surrogate spike approach, significant differences between datasets were found at RoA $\geq$ 50\% (Wilcoxon rank-sum, $p < .05$). P Values at each level are as follows: $p_{10\%}=0.3444$, $p_{20\%}=0.7435$, $p_{30\%}=0.2552$, $p_{40\%}=0.0902$, $p_{50\%}=0.0311$, $p_{60\%}=0.0077$, $p_{70\%}=0.0020$, $p_{80\%}=0.0020$, $p_{90\%}=0.0189$. The effect sizes at each level are as follows: $r_{10\%} = 0.0078$, $r_{20\%} = 0.0027$, $r_{30\%} = -0.0094$, $r_{40\%} = -0.0140$, $r_{50\%} = -0.0178$, $r_{60\%} = -0.0220$, $r_{70\%} = -0.0254$, $r_{80\%} = -0.0255$, $r_{90\%} = -0.0193$. All effect sizes fall within the small effect size range ($|r| \leq 0.1$), indicating that while the differences are statistically significant given the large sample sizes ($n = 10{,}429$ for 1DoF; $n = 4{,}296$ for GRASP), they are not practically large. \par For the subset spike approach, significant differences were observed across all nine RoA levels (Wilcoxon rank-sum, $p < .05$). P Values at each level are as follows: $p_{10\%}=6.15\times10^{-19}$, $p_{20\%}=5.28\times10^{-19}$, $p_{30\%}=3.35\times10^{-18}$, $p_{40\%}=2.22\times10^{-18}$, $p_{50\%}=1.77\times10^{-17}$, $p_{60\%}=9.31\times10^{-19}$, $p_{70\%}=3.82\times10^{-19}$, $p_{80\%}=5.27\times10^{-19}$, $p_{90\%}=1.63\times10^{-19}$. The effect size is again consistently small: $r_{10\%} = 0.0078$, $r_{20\%} = 0.0027$, $r_{30\%} = -0.0094$, $r_{40\%} = -0.0140$, $r_{50\%} = -0.0178$, $r_{60\%} = -0.0220$, $r_{70\%} = -0.0254$, $r_{80\%} = -0.0255$, $r_{90\%} = -0.0193$. \par The negative direction of the effect sizes at higher RoA levels indicates that XCC values from the GRASP dataset tend to be marginally higher than those from the 1DoF dataset at high RoA, consistent with the visual pattern in Supplementary Fig. \ref{supfigR1}. \par These small but consistent effect sizes across both simulation approaches and across the full RoA range reinforce the conclusion that the RoA-XCC relationship is dataset-dependent, even if the magnitude of the difference is modest. This finding is practically significant because it implies that a waveform-tracking threshold calibrated on one dataset may not transfer optimally to another, even when the same muscles and task structure are involved. % ----------------------------------------------------------------------- \subsection*{Supplementary Note 3: Subset spike approach validation}\label{SuppNote3} \hspace{15pt} To verify that the findings from the surrogate spike approach are not an artifact of the spike generation procedure, we repeated the simulation using a subset spike approach, in which artificial MUs were generated using only a subset of the decomposed spikes rather than augmented surrogate spikes (see Methods \ref{Method-S5-SS1}). This approach introduces less noise at low RoA levels and provides independent validation of the surrogate spike findings. All four main findings from the surrogate spike approach were confirmed in the subset spike approach, as detailed below. \par \textbf{Nonlinear RoA-XCC correlation} The nonlinear relationship between RoA and XCC was confirmed in the subset spike approach for both datasets (Supplementary Fig. \ref{supfigR1}a). The best-fitting model differed from the surrogate spike approach: a log-transformed model provided the best fit under the subset spike procedure, as determined by BIC. The detailed non-linear relationships are reported in Supplementary Note 1. Despite this difference in functional form, the qualitative character of the relationship, with XCC increasing nonlinearly with RoA and exhibiting higher variability at low RoA levels, was preserved across both simulation approaches. \par \textbf{Dataset-dependent pattern} The dataset-dependent pattern in the RoA-XCC correspondence was confirmed in the subset spike approach. XCC values from the 1DoF and GRASP datasets differed significantly across all nine RoA levels (Wilcoxon rank-sum, $p < .05$ at all levels), though effect sizes remained small throughout (Cohen's $r < 0.1$ at all levels). The full effect size values for both simulation approaches across all nine RoA levels are reported in Supplementary Note 3. \par \textbf{Decreasing XCC variability} The finding that XCC variability decreases with increasing RoA was also confirmed in the subset spike approach. The regression statistics for both approaches are reported together in the main text (Section \ref{Result-S1-SS2}) for direct comparison. \par \textbf{Inevitable disagreement} The finding that disagreement between spike train matching and waveform tracking is inevitable under commonly used thresholds was confirmed in the subset spike approach. After restricting analysis to MU pairs with RoA $\leq$ 30\%, agreement under XCC thresholds of 0.7 and 0.9 followed the same pattern as in the surrogate spike approach (Supplementary Fig. \ref{supfigR1}). \par Taken together, these results confirm that the main findings are robust to the choice of spike generation procedure and are not an artifact of the surrogate spike simulation method. % ----------------------------------------------------------------------- \begin{figure}[htbp!] \centering \includegraphics[width=1.0\textwidth]{SupFig/SupFigR2_overlay.png} \caption{Schematic of the overlaid plot applied in the joint distribution evaluation and edge-level analysis. Each coordinate plane represents a pairwise comparison within a single decomposed MU set. To study the whole dataset's behavior, we overlay all comparisons.}\label{supfigR2_overlay} \end{figure} \begin{figure}[htbp!] \centering \includegraphics[width=1.0\textwidth]{SupFig/SupFigR2_example.png} \caption{Examples for MU pairs whose RoA $<$ 10\% and XCC $>=$ 0.9. 1\textsuperscript{st} and 2\textsuperscript{nd} MU: 1DoF Dataset, subject 9, day 1, finger 4, trial 2 extensor muscle, MU 5 and MU 23.}\label{supfigR2_example} \end{figure} % ----------------------------------------------------------------------- \subsection*{Supplementary Note 4: GMM fitting and BIC analysis for marginal distributions}\label{SuppNote4} \hspace{15pt} To characterize the modality of the RoA and XCC marginal distributions, we first confirmed non-unimodality using the Hartigan dip test (main text, Section \ref{Result-S2-SS2}), then estimated the number of modes by fitting GMMs with one to six components and selecting the optimal component count using the Bayesian Information Criterion (BIC; see Methods \ref{Method-S5-SS3}). \par The BIC decreased monotonically as the number of components increased for all four marginal distributions, namely 1DoF-RoA, 1DoF-XCC, GRASP-RoA, and GRASP-XCC, across both datasets (Fig. \ref{figR2}c). This pattern indicates that the GMM fitting procedure did not identify a clear optimal modal structure: each additional component continued to improve the BIC without converging to a stable solution within the range tested (one to six components). As a result, no definitive modal size can be concluded from this analysis. \par This finding has an important practical implication. Several prior studies have proposed determining spike train matching or waveform tracking thresholds analytically from the distribution of RoA or XCC values, based on the assumption that a bimodal distribution would reveal a natural separation between same and different MU pairs \cite{chen2022caution}. The absence of a clear bimodal or low-modal structure in the empirical distributions of both RoA and XCC, confirmed here across two independent datasets and two metrics, indicates that this approach is not reliably applicable in practice. The non-unimodal but indeterminate modal structure of these distributions further supports the conclusion that threshold determination requires dataset-specific calibration rather than distributional analysis alone. \par The marginal distributions for both datasets and both metrics are shown in Fig. \ref{figR2}a, b, with BIC curves in Fig. \ref{figR2}c. % ------------------------------------------- \subsection*{Supplementary Note 5: Full optimal waveform-tracking threshold trajectories across all spike train matching thresholds}\label{SuppNote5} \hspace{15pt} In the main text, we reported the waveform-tracking threshold that maximizes edge-, node-, and group-level agreement at the commonly used spike train matching threshold of 30\% RoA, as well as the globally optimal threshold combination for each dataset. Here, we provide the full trajectory of the optimal waveform-tracking threshold at every fixed spike train matching threshold, for all three agreement levels and both datasets. \subsubsection*{Edge-level agreement} \hspace{15pt} For the 1DoF dataset, the optimal waveform-tracking threshold increases monotonically as the spike train matching threshold varies from 1\% to 54\% RoA, then plateaus and remains stable once the spike train matching threshold exceeds 54\%. For the GRASP dataset, the optimal waveform-tracking threshold is largely stable: it remains at XCC = 0.81 as the spike train matching threshold varies from 4\% to 29\%, increases to XCC = 0.86 as the threshold varies from 30\% to 38\%, and stabilizes once the spike train matching threshold exceeds 38\%. The divergence in these trajectories between the two datasets is a further manifestation of the dataset-dependent pattern described in the main text (Section \ref{Result-S2-SS3}). \subsubsection*{Node-level agreement} \hspace{15pt} For the 1DoF dataset, the optimal waveform-tracking threshold is XCC = 0.50 at a spike train matching threshold of 5\%. It then remains at XCC = 0.80 as the spike train matching threshold varies from 10\% to 35\%, increases from XCC = 0.85 to 0.95 as the spike train matching threshold varies from 40\% to 55\%, and stabilizes at XCC = 0.95 as the spike train matching threshold varies from 60\% to 95\%. For the GRASP dataset, the trajectory is more complex: at a 5\% spike train matching threshold, the optimal waveform-tracking threshold is XCC = 0.75. Between 10\% and 15\%, it takes values between 0.85 and 0.90; between 15\% and 20\%, values between 0.80 and 0.90. It then stabilizes at XCC = 0.85 from 25\% to 45\%, and rises to XCC = 0.95 from 50\% to 95\%. These patterns are consistent with the node-level agreement results reported in the main text (Section \ref{Result-S2-SS4}). \subsubsection*{Group-level agreement} \hspace{15pt} For the 1DoF dataset, the optimal waveform-tracking threshold is XCC = 0.75 at a spike train matching threshold of 5\%, remains at XCC = 0.80 from 10\% to 30\%, increases from XCC = 0.80 to 0.95 as the spike train matching threshold varies from 30\% to 45\%, and stabilizes at XCC = 0.95 from 45\% to 95\%. For the GRASP dataset, the optimal waveform-tracking threshold is XCC = 0.75 at 5\%, takes values between 0.75 and 0.90 from 10\% to 20\%, remains at XCC = 0.90 from 20\% to 35\%, and stabilizes at XCC = 0.95 from 40\% to 95\%. These patterns are consistent with the group-level agreement results reported in the main text (Section \ref{Result-S2-SS5}). \par Across all three agreement levels and both datasets, two consistent patterns emerge. First, the optimal waveform-tracking threshold is non-decreasing as the spike train matching threshold increases, reflecting the positive RoA-XCC correlation established in the simulation. Second, the specific threshold values and the RoA level at which the optimal XCC stabilizes differ between datasets, reinforcing the conclusion that no single waveform-tracking threshold is universally optimal across datasets or spike train matching thresholds. These trajectories are provided to assist researchers in selecting threshold combinations appropriate for their specific dataset, pending a dedicated pilot study. \end{document}