Temporally stable beta sensorimotor oscillations and corticomuscular coupling underlie force steadiness

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY-NC-ND 4.0 https://creativecommons.org/licenses/by-nc-nd/4.0/ Temporally stable beta sensorimotor oscillations and corticomuscular coupling underlie force steadiness © 2022 The Authors. Published by Elsevier Inc. Published version Mongold, Scott J.; Piitulainen, Harri; Legrand, Thomas; Ghinst, Marc Vander; Naeije, Gilles; Jousmäki, Veikko; Bourguignon, Mathieu Mongold, S. J., Piitulainen, H., Legrand, T., Ghinst, M. V., Naeije, G., Jousmäki, V., & Bourguignon, M. (2022). Temporally stable beta sensorimotor oscillations and corticomuscular coupling underlie force steadiness. Neuroimage, 261, Article 119491. https://doi.org/10.1016/j.neuroimage.2022.119491 2022 NeuroImage 261 (2022) 119491 Contents lists available at ScienceDirect NeuroImage journal homepage: www.elsevier.com/locate/neuroimage Temporally stable beta sensorimotor oscillations and corticomuscular coupling underlie force steadiness Scott J. Mongold a , ∗ , Harri Piitulainen b , c , Thomas Legrand a , Marc Vander Ghinst d , e , Gilles Naeije d , f , Veikko Jousmäki c , g , Mathieu Bourguignon a , d , h a Laboratory of Neurophysiology and Movement Biomechanics, UNI –ULB Neuroscience Institute, Université libre de Bruxelles (ULB), Brussels, Belgium b Faculty of Sport and Health Sciences, University of Jyväskylä, Jyväskylä, Finland c Department of Neuroscience and Biomedical Engineering, Aalto University School of Science, Espoo, Finland d Laboratoire de Cartographie fonctionnelle du Cerveau, UNI –ULB Neuroscience Institute, Université libre de Bruxelles (ULB), Brussels, Belgium e Service d’ORL et de chirurgie cervico-faciale, CUB Hôpital Erasme, Université libre de Bruxelles (ULB), Brussels, Belgium f Centre de Référence Neuromusculaire, Department of Neurology, CUB Hôpital Erasme, Université libre de Bruxelles (ULB), Brussels, Belgium g Aalto NeuroImaging, Aalto University School of Science, Espoo, Finland h BCBL, Basque Center on Cognition, Brain and Language, 20009 San Sebastian, Spain a r t i c l e i n f o Keywords: Mu rhythm Beta sensorimotor oscillations Corticomuscular coherence Isometric contraction Magnetoencephalography Motor controling Primary sensorimotor cortex Muscle electromechanical coupling a b s t r a c t As humans, we seamlessly hold objects in our hands, and may even lose consciousness of these objects. This phenomenon raises the unsettled question of the involvement of the cerebral cortex, the core area for voluntary motor control, in dynamically maintaining steady muscle force. To address this issue, we measured magnetoencephalographic brain activity from healthy adults who maintained a steady pinch grip. Using a novel analysis approach, we uncovered fine-grained temporal modulations in the beta sensorimotor brain rhythm and its coupling with muscle activity, with respect to several aspects of muscle force (rate of increase/decrease or plateauing high/low). These modulations preceded changes in force features by ∼40 ms and possessed behavioral relevance, as less salient or absent modulation predicted a more stable force output. These findings have consequences for the existing theories regarding the functional role of cortico-muscular coupling, and suggest that steady muscle contractions are characterized by a stable rather than fluttering involvement of the sensorimotor cortex. 1. Introduction As humans, we rely on our hands to interact with the environment, using them to communicate, touch, and importantly, hold items, i.e., phone, coffee, and keys. Remarkably, we are able to lose awareness of the very object in our hand, even as we maintain grip, begging the question: is voluntary control necessary for sustained, low intensity steady contractions? The sensorimotor cortex is unequivocally implicated in voluntary muscle contraction, yet, its role in maintaining steady contractions is unsettled. Does it play a sustained role in this highly dynamic process, or a phasic role in correcting when the goal is no longer matched (the phone is slipping off the hand)? At least in animals, stereotyped motor actions such as walking do not require corticomuscular communication after initiation ( Purves, 1999 ). When attempting to sustain an isometric contraction, the applied force is never constant but rather fluctuates around an average value (as reviewed in Enoka and Farina, 2021 ). Typically, force variability is quantified over an entire isometric contraction, usually main- ∗ Corresponding author. E-mail address: [email protected] (S.J. Mongold) . tained over several seconds ( Jones et al., 2002 ; Laidlaw et al., 2000 ; Tracy and Enoka, 2002 ; Ushiyama et al., 2017 ). Surprisingly, there has been little investigation into non-global, ‘dynamic’ measures of force variability. A previous report indicated that force fluctuations during isometric contractions and brain activity are phase-coupled at frequencies below 3 Hz (Bourguignon et al., 2017) , a coupling akin to that occuring with submovements at 1–4 Hz during slow tracking movements ( Dipietro et al., 2011 ; Hall et al., 2014 ) and to that termed the corticokinematic coherence that occurs during fast repetitive finger movements at movement frequency and harmonics ( Bourguignon et al., 2012 , 2011 ; Piitulainen et al., 2013a ). It was suggested to reflect the processing of proprioceptive afferents’ signaling based on the results of an effective coupling analysis ( Bourguignon et al., 2015 , 2017 , 2019 ; Piitulainen et al., 2013a , 2013b ) Still, the cortical mechanisms underlying regulation (as opposed to monitoring) of force fluctuations are unclear. One means to assess corticomuscular communication is corticomuscular coherence (CMC), a coupling clearly distinct from corticokinematic coherence in terms of frequency channel implicated and reactivhttps://doi.org/10.1016/j.neuroimage.2022.119491 . Received 22 February 2022; Received in revised form 12 July 2022; Accepted 15 July 2022 Available online 28 July 2022. 1053-8119/© 2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 ity to movements ( Bourguignon et al., 2019 ). CMC captures the phase coupling that occurs between brain and muscle activities, mainly at beta ( Baker et al., 1997 ; Conway et al., 1995 ; Salenius et al., 1997 ) and alpha frequencies ( Piitulainen et al., 2015a ), i.e., the main components of the sensorimotor rhythm that reflect the state of activation of sensorimotor cortices ( Pineda, 2005 ). CMC usually peaks during sustained isometric contractions, decreases during dynamic contractions, and possesses somatotopic representation in the primary motor cortex and primary somatosensory cortex contralateral to the contracted muscle ( Hari and Salenius, 1999 ; Kilner et al., 1999 ; Salenius et al., 1997 ; Salenius and Hari, 2003 ). A host of studies suggests it builds on the descending motor command (Bourguignon et al., 2019) , but is modulated by (re)afferent information ( Fisher et al., 2002 ; Kilner et al., 2004 ; Liu et al., 2019 ; Riddle and Baker, 2005 ). Previous studies have focused on the association between global measures of force stability and CMC magnitude, based on minute-long recordings. Studies assessing motor precision where CMC levels are compared between conditions suggest that increased CMC is associated with smaller errors between target and exerted forces ( Kristeva et al., 2007 ; Mendez-Balbuena et al., 2012) . Somewhat contrastingly, in studies assessing correlation across participants, CMC appears positively associated with the amplitude of force fluctuations ( Ushiyama et al., 2017 , 2011a ). However, no matter the reasons for the discrepancy, these studies did not look at the temporal dynamics of CMC in relation to force fluctuations, which is key to clarify the cortical involvement in force regulation. Existing CMC analysis methods do not allow for the study of force regulation. CMC is typically estimated based on second-long epochs ( Mendez-Balbuena et al., 2012 ; Ushiyama et al., 2017 , 2010 ). This has allowed for assessment of CMC modulation in response to wellcontrolled isolated events, such as force ramps ( Kilner et al., 2003 , 2000 ) or sensory stimulations ( Hari et al., 2014 ; Piitulainen et al., 2015 b), but not continuously throughout a contraction. Here, we introduce a novel analysis method to identify the temporal dynamics of CMC and brain rhythms in relation to continuous signals. We use this method to determine whether and how beta CMC and beta brain rhythm modulate in relation to force fluctuations during volitional contraction. We expected to observe significant modulation whereby increased CMC and beta brain rhythm lead to increased force. Indeed, CMC and beta rhythm are generally correlated (van Wijk et al., 2009 ) and the latter is known to be transmitted down the corticospinal tract to spinal motoneurons ( Baker et al., 2003 ). Since motor unit pools possess a close-to-linear transfer function ( Stegeman et al., 2010 ), these beta oscillations should modulated EMG power, resulting in corresponding changes in force (McAuley et al., 1997) . The main aim of the present study is to provide direct evidence that such modulations exist, arguing for a sustained role of the cortex in regulating steady contraction force. If so, we also aim to determine (i) the temporal dynamics of these modulations with respect to force fluctuations, (ii) their relevance for force steadiness, and (iii) how they relate to global cortical involvement in the task (global CMC and beta power depression or enhancement). 2. Materials and Methods This is a reanalysis of previously published data (Bourguignon et al., 2017) . 2.1. Participants Seventeen healthy human volunteers (7 females, 10 males; mean ± SD, age 34 ± 7 years, range 20–47 years) with no history of neuropsychiatric diseases or movement disorders participated in our study. All participants were right-handed (mean ± SD, score 90 ± 12, range 65– 100 on the scale from –100 to 100; Edinburgh handedness inventory; Oldfield, 1971 ). The study had prior approval by the ethics committee of the Helsinki and Uusimaa hospital district. The participants gave informed consent before participation, and they were compensated monetarily for travel expenses and lost working hours. 2.2. Experimental protocol Figure 1 illustrates the experimental paradigm. During MEG recordings, the participants were sitting with their left hand on the thigh and their right hand on a table in front of them. Participants’ vision was optically corrected with nonmagnetic goggles when needed. Participants were asked to maintain a steady isometric pinch grip of 2–4 N against a custom-made handgrip (connected to a rigid load cell; rigidity 15.4 N/mm; Model 1004, Vishay Precision Group, Malvern, PA, USA) with the right thumb and the index finger (see Fig. 1 A), and to fixate at a black cross displayed on the center of a screen placed 1 m in front of them. When the force stepped out of the prescribed limits, a triangle (pointing up or down) appeared on top of the black cross, indicating in which direction to adjust the force, and it disappeared as soon as the force was correctly returned within the limits (see Fig. 1 B and C). After a ∼1-min practice session, two 5-min blocks were recorded, with a minimum of 2-min rest between the blocks. Each block started with ∼10 s without contraction, after which participants were prompted to begin the contraction task. A 5-min task-free block was also recorded. Of note, a wide range of force (2–4 N) was chosen to ensure participants could maintain the contraction with only infrequent corrective feedback. Moreover, in practice, the force remained relatively steady, with fluctuations at frequencies above 1 Hz lower than 100 nN during periods when visual feedback did not prompt for corrections (Bourguignon et al., 2017) . 2.3. Measurements MEG. The MEG measurements were carried out in a magnetically shielded room (Imedco AG, Hägendorf, Switzerland) at the MEG Core of Aalto NeuroImaging, Aalto University, with a 306-channel wholescalp neuromagnetometer (Elekta Neuromag TM , Elekta Oy, Helsinki, Finland). The recording passband was 0.1–330 Hz and the signals were sampled at 1 kHz. Participants’ head position inside the MEG helmet was continuously monitored by feeding current into four head-tracking coils located on the scalp; the locations of the coils and at least 200 headsurface points (scalp and nose) with respect to anatomical fiducials were determined with an electromagnetic tracker (Fastrak, Polhemus, Colchester, VT, USA). EMG and force. Surface EMG was measured from the first dorsal interosseous . Active EMG electrodes were placed on the muscle bulk and signals were measured with respect to a passive electrode placed over the distal radial bone. Recording passband was 10–330 Hz for EMG signals and DC–330 Hz for the force signal. EMG and force signals were then sampled at 1 kHz and recorded time-locked to MEG signals. MRI. 3D-T1 magnetic resonance images (MRIs) were acquired with General Electric Signa® 3.0 T whole-body MRI scanner (Signa VH/i, General Electric, Milwaukee, WI, USA) or with 3T MAGNETOM Skyra whole-body MRI scanner (Siemens Healthcare, Erlangen, Germany) at the AMI Centre, Aalto NeuroImaging, Aalto University School of Science. 2.4. MEG preprocessing Continuous MEG data were preprocessed off-line with MaxFilter 2.2.10 (Elekta Oy, Helsinki, Finland), including head movement compensation. The tSSS preprocessing was applied with a correlation limit of 0.9 and segment length equal to the recording length (Taulu and Kajola, 2005; Taulu and Simola, 2006). Raw data files were opened with MNE Matlab ( Gramfort et al., 2014 ), and unless stated otherwise, all further analyzes were conducted with custom-made Matlab (Mathworks, Natick, MA) scripts. Independent 2 S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 Fig. 1. Experimental setup. A, Illustration of the isometric contraction task. A steady contraction is maintained on a custommade handgrip with the right thumb and the index finger. Surface EMG is measured from the first dorsal interosseous (top right electrode) and the flexor carpi ulnaris (not visible here) of the right hand, with reference electrode over the distal radial bone (bottom left electrode). B, Sixty seconds of raw force signal from a representative participant. The participant was prompted to start contracting 10 s after the beginning of the recording. Two horizontal dashed lines indicate the force limits (2– 4 N). Gray shaded areas represent periods wherein contraction force was out of the prescribed bounds for at least 1 time-bin 2 s around. Corresponding data were not analyzed. C, Visual feedback presented to the participants to help them regulate their contraction force. Cross on the screen informed them that the force was within the prescribed limits. Arrow pointing up (respectively down) prompted them to increase (respectively decrease) the force. Reproduced from Bourguignon et al. (2017) . component analysis was then applied to MEG signals filtered through 1– 25 Hz using FastICA algorithm (dimension reduction, 30; non-linearity, tanh) ( Hyvärinen et al., 2001 ; Vigario et al., 2000 ). Between 1 and 3 components corresponding to eye-blink and heartbeat artifacts were visually identified based on their topography and time-series. The corresponding components were subsequently subtracted from full-band and full-rank MEG signals. 2.5. Conventional CMC estimation We used coherence analysis to estimate CMC for all MEG sensors, and to identify the optimal MEG sensor for further analyses. Time points for which visual feedback was presented (force not properly kept between 2 and 4 N) were marked as bad to exclude periods during which contraction was intentionally corrected. Time points for which MEG signals exceeded 5 pT (magnetometers) or 1 pT/cm (gradiometers) were also marked as bad to avoid contamination of the data by any artifact not removed by the pre-processing. Continuous data from the recording blocks were split into 1000-ms epochs with 800-ms epoch overlap (Bortel and Sovka, 2007), leading to a frequency resolution of 1 Hz. Epochs less than 2-s away from time points marked as bad were discarded from further analyses (mean ± SD artifact-free epochs 2875 ± 130, range 2573–3031). Coherence spectra were computed between all MEG sensors and the non-rectified EMG signal following the formulation of Halliday et al. (1995) , and by using the multitaper approach (5 orthogonal Slepian tapers, yielding a spectral smoothing of ± 2.5 Hz) to estimate power- and cross-spectra (Thomson, 1982 ). Whether EMG data should be rectified or not is still debated ( Halliday and Farmer, 2010 ; McClelland et al., 2014 , 2012 ; Ward et al., 2013 ). In any case, reasonable CMC estimates were obtained here without rectification, as reflected by the large proportion of participants in whom it was significant (16 in 17; see results section). Data from gradiometer pairs were combined in the direction of maximum coherence as done in ( Bourguignon et al., 2015 ). Then, for each participant, the gradiometer pair with the highest coherence value in the 10–30-Hz band was selected among a predefined subset of 9 gradiometer pairs covering the primary sensorimotor (SM1) cortex. Gradiometer selection was based on the fact that CMC generally peaks topographically in the sensors overlying the primary sensorimotor region, contralateral to the side of the contracting muscle, as done previously ( Beck et al., 2021 ). Further confirming the validity of this approach, the gradiometer pair with maximum CMC was among this preselection in 16 out of 17 participants. Further analyses were performed with the selected gradiometer signal in the orientation yielding the maximum coherence (MEG SM1 ). 2.6. Force features We extracted four features from the contraction force signal to represent the most relevant aspects of force fluctuations: rate of change (increase and decrease separately) and plateauing (at a low or high level). First, the force signal was low-pass filtered at 30 Hz, as previously done in studies investigating force fluctuations during isometric contractions ( Danion and Galléa, 2004 ; Missenard et al., 2009 ; Witte et al., 2007 ) . Then, time points less than 2 s away from periods where the force changed over 50 ms by more than 5 SD (range for this threshold across participants, 0.13–0.35 N; mean ± SD, 0.25 ± 0.06) were marked as bad (and hence excluded from the analysis), completing those for which MEG amplitude was too high or the force was not in the prescribed range (2–4 N). From the low-pass filtered force signal F ( t ), we estimated the rate of force change as  𝐹 ( 𝑡 ) = ( 𝐹 ( 𝑡 ) − 𝐹 ( 𝑡 − 𝛿𝑡 ) )∕ 𝛿𝑡 where 𝛿𝑡 is the time interval between adjacent samples (1 ms here). The two first force features considered were derived from this rate of force change (  𝐹 ): 1) The rate of force increase signal was the rate of force change half-wave rectified. That is  𝐹 when  𝐹 is positive and 0 otherwise. 2) The rate of force decrease signal was the opposite of the rate of force change half-wave rectified. That is −  𝐹 when  𝐹 is negative and 0 otherwise. The two additional force features considered were derived from the force plateauing signal, computed as the inverse of the rate of force change full-wave rectified ( 1∕ | 𝐹 |). To render the inversion computationally stable, values of the force steadiness above the 95 th percentile were set to that percentile value. Of note, the choice of the clipping cutoff (percentile 90, 95 or 99) had virtually no impact on the results. These features were: 3 S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 3) The force plateauing high signal was the force plateauing signal at time points where the force trace was concave (second derivative of F,  𝐹 ( 𝑡 ) = ( 𝐹 ( 𝑡 − 𝛿𝑡 ) − 2 𝐹 ( 𝑡 ) + 𝐹 ( 𝑡 − 𝛿𝑡 ) )∕ 𝛿𝑡 2 , positive) and 0 elsewhere. 4) The force plateauing low signal was the force plateauing signal at time points where the force trace was convex (second derivative of F negative) and 0 elsewhere. 2.7. CMC and power modulation in relation to force features We next estimated how CMC and the power of MEG and EMG signals modulate in relation to the four force features. Key to this analysis is the use of Hilbert transformation to estimate coherence. In what comes next, we first present how the Hilbert transformation can be used to estimate coherence and power. We then generalize this approach to introduce a weighting by positive time-series, such as each of the 4 force features presented in the previous section, to estimate how CMC and power modulate in relation to such time-series. The rationale is that, if such weighting enhances CMC (or power), this means time periods when the time-series is high concur with time periods when CMC (or power) is high. Finally, we present the inclusion of time delays in force features to get at the temporal dynamics of CMC and power modulations. To estimate CMC with the Hilbert transformation, MEG SM1 and EMG signals were filtered through a 5-Hz-wide frequency band centered on individual peak CMC frequency. The band-pass filter used in that effect was designed in the frequency domain with zero-phase and 1-Hz-wide squared-sine transitions from 0 to 1 and 1 to 0 (e.g. a filter centered on 20 Hz rose from 0 at 17 Hz to 1 at 18 Hz and ebbed from 1 at 22 Hz to 0 at 23 Hz). To remove the impact of artifacts, further analyses were based on time-points at least 2 s away from MEG artifacts (time points of MEG signals exceeding 5 pT for magnetometers or 1 pT/cm for gradiometers) and appearance of visual feedback. The analytical signals for MEG SM1 and EMG were then created by means of the Hilbert transform. From these analytical signals denoted s 1 ( t ) for MEG SM1 and s 2 ( t ) for EMG, the following formula 𝑃 𝑖𝑗 = ⟨𝑠 𝑖 ( ⋅) 𝑠 𝑗 ( ⋅) ∗ ⟩(1) yielded MEG SM1 power ( i = j = 1), EMG power ( i = j = 2) and their cross-power ( i = 1, j = 2). Here, the mean denoted ⟨⋅⟩was taken across all artifact-free samples, and ∗ denotes complex conjugation. From these quantities, the —unweighted —coherence is estimated as Coh = |𝑃 12 |2 𝑃 11 𝑃 22 (2) We used these formulas as a basis to introduce weighting by any positive time-series, which in our case will be each of the four force features. Denoting w ( t ) such a positive time-series time-locked to MEG and EMG signals, weighted auto- and cross-power can be estimated as 𝑃 𝑤 𝑖𝑗 = ⟨𝑠 𝑖 ( ⋅) 𝑠 𝑗 ( ⋅) ∗ 𝑤 ( ⋅) ⟩∕ ⟨𝑤 ( ⋅) ⟩(3) and weighted coherence as 𝐶𝑜ℎ 𝑤 = |𝑃 𝑤 12 |2 𝑃 𝑤 11 𝑃 𝑤 22 (4) Notice that when w is constant, formulas (1) and (3) are equivalent. The interpretation of these weighted measures for non-constant w is straightforward. For example, an increase in Coh w compared with Coh indicates that CMC tends to be higher when the time-series considered ( w ) is high. The last key element to our derivation of weighted coherence is the introduction of a time delay 𝜏in w ( t ), so that 𝑤 𝜏( 𝑡 ) = 𝑤 ( 𝑡 − 𝜏) . From there, the temporal evolution of MEG SM1 and EMG power and their coherence with respect to w can be assessed as 𝑃 𝑤 𝑖𝑗 ( 𝜏) = ⟨𝑠 𝑖 ( ⋅) 𝑠 𝑗 ( ⋅) ∗ 𝑤 𝜏( ⋅) ⟩∕ ⟨𝑤 𝜏( ⋅) ⟩(5) and 𝐶𝑜ℎ 𝑤 ( 𝜏) = |𝑃 𝑤 12 ( 𝜏) |2 𝑃 𝑤 11 ( 𝜏) 𝑃 𝑤 22 ( 𝜏) (6) With these formulas, a peak increase in Coh w for a given force feature w at time delay, e.g. , 𝜏= 40 ms, indicates that the force feature considered leads to an increase in CMC 40 ms later. Conversely, a peak in Coh w at a negative time delay 𝜏indicates that high CMC leads to an increase in the force feature considered with time delay 𝜏. The analysis approach described above was used to estimate how beta CMC, EMG power, and MEG SM1 power modulate in relation to the four selected force features. In that analysis, modulations were computed for time delays ( 𝜏) ranging from –2 s to 2 s by steps of 10 ms. Power modulations were normalized by the mean across the baseline defined by 1 s < |𝜏|< 2 s. We also used this approach to estimate how full-band EMG power modulates in relation to the four force features. Full-band EMG power was obtained from the EMG signal high-pass filtered at 30 Hz, rectified, low-pass filtered at 50 Hz, and then squared. In that analysis, the temporal resolution ( 𝜏) was set to 1 ms to resolve the fast transient dynamics of EMG-force coupling. 2.8. Link between identified modulations and behavioral relevance We used cross-participant Pearson correlation to determine the degree of (in)dependence between and among the different modulations we observed, parameters quantifying force steadiness, and global MEG SM1 power suppression. These different measures are described below. From CMC and power modulations, we extracted the maximum (for increases) or minimum (for decreases) value for each individual within ± 50 ms around group-level peak values. These values were then converted to percentage of change relative to baseline (1 s < | 𝜏| < 2), giving rise to CMC and power modulation values. For regularization purposes, peak and baseline CMC values were incremented by 0.01 before the division. Force steadiness was estimated as the standard deviation (SD) of the 10-min of force signal filtered in three relevant frequency ranges: 0.5–3 Hz, 3–15 Hz and 15–30 Hz. Force fluctuations at 0.5–3 Hz were previously argued to be the most relevant for the maintenance of a steady contraction and to be monitored by the brain through the proprioceptive system (Bourguigno n et al., 2017) . Force fluctuations in the intermediate range (3–15 Hz) should relate more to intermittent motor control ( Gross et al., 2002 ) and capture the physiological tremor at ∼10 Hz ( Gilbertson, 2005 ; McAuley et al., 1997 ). Force fluctuations at 15–30 Hz were shown to be tightly linked to the presence of beta CMC (Bourguignon et al., 2017) ; indeed, because the beta sensorimotor rhythm shapes the corticospinal drive and the periphery receives this drive, it follows that EMG and force fluctuations at beta frequencies emulate those of the cortex, as has been suggested previously ( Gross et al., 2000 ). Global MEG SM1 power suppression was estimated as the MEG SM1 power at the frequency of peak CMC in the isometric contraction task divided by this same quantity estimated from the task-free recording (preprocessed in the same way as isometric contraction MEG data). In cross-participant Pearson correlation, data points were considered as outliers if they departed by over 2 SDs from the mean (z-score above 2). A maximum of one outlier was left out (the one with the highest z-score) in each correlation analysis. 2.9. Statistical analyses Significance of unweighted coherence. A threshold for statistical significance of the coherence ( p < 0.05 corrected for multiple comparisons) was obtained as the 95 th percentile of the distribution of the maximum coherence —across 10–30 Hz, and across the selection of 9 gradiometers —evaluated between MEG and Fourier-transform surrogate reference signals (1000 repetitions; Faes et al., 2004) . The Fourier-transform 4 S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 surrogate of a signal is obtained by computing its Fourier transform, replacing the phase of the Fourier coefficients by random numbers in the range [ − 𝜋; 𝜋], and then computing the inverse Fourier-transform ( Faes et al., 2004 ; Theiler et al., 1992 ). Of note, the proportion of statistically significant CMC values was already reported in our previous study (Bourguignon et al., 2017) , reaching 16 in 17 participants, which compares favorably to other findings ( Johnson and Shinohara, 2012 ) and largely surpasses the ∼50% in other studies ( Steeg et al., 2014 ). This is likely explained by several elements, such as the use of MEG that typically features ∼3 times higher signal to noise ratio than EEG ( Destoky et al., 2019 ), the use of 10-minute-long recordings (vs. 5- minute-long in many studies) and the use of 5-Hz spectral smoothing (vs. 1 Hz in many studies), with the latter two factors reducing the significance threshold by tenfold. Significance of CMC and power modulation by force features. For each participant and each force feature, we used surrogate-data-based statistics to estimate the significance of the modulation in the 3 following signals: beta CMC, beta MEG SM1 power, and beta EMG power. Amplitude increase and decrease were defined as the maximum and minimum (respectively) difference between modulation in the range -0.5 s to 0.5 s (| 𝜏| < 0.5 s) and mean baseline value (0.5 s < | 𝜏| < 1 s). Thresholds for statistical significance of amplitude increase and decrease ( p < 0.05) were obtained as the 97.5 th percentile of their distribution for surrogate force features (500 repetitions). Surrogate force features were obtained by applying a circular shift of multiples of 1 s to the genuine force features. Non-parametric permutation statistics were further used to estimate the significance at the group level of modulations by force features ( Nichols and Holmes, 2002 ). This analysis aimed at determining if the modulation of each of the 3 modulated signals by each of the four force features was consistent across participants. First, we estimated the modulation amplitude for the modulation signals —with baseline level set to zero —averaged across participants. A permutation distribution (10,000 permutations) for this modulation amplitude was then derived from modulation traces in which the sign was randomly changed for each participant. An exact p -value for each sample of the non-permuted data was obtained as the proportion of the values in the permutation distribution exceeding the value at this sample (Nichols and Holmes, 2002) . Time points corresponding to significant modulations ( p < 0.05) were identified. Confidence interval for increases and decreases in CMC and power modulations. We used bootstrap statistics ( Efron, 1979 ; Efron and Tibshirani, 1994 ) to determine a 95% confidence interval for the timing and amplitude of peaks and troughs of modulations of beta CMC, beta MEG SM1 power, beta EMG power, and full-band EMG power in relation to force features. The bootstrap applied was bias-corrected and accelerated with 1000 resamplings (Efron, 1979; Efron and Tibshirani, 1994) . It was applied to the time courses upsampled by a factor 10 using spline interpolation (i.e., 0.1 ms resolution for EMG full band, and 1 ms resolution for the 3 other modulation signals). Global significance of associations with multiple modulation values. Modulations in beta CMC and MEG SM1 power were revealed to be highly correlated. Hence, we jointly estimated the significance of their correlation with other parameters (e.g., force SD or global MEG SM1 power attenuation) using non-parametric permutation statistics. Individual correlation values were averaged after having rendered their signs consistent (i.e., modulations were multiplied by -1 when they showed a negative correlation with an arbitrary reference modulation, the choice of which had no impact on the results). A permutation distribution was built for the averaged correlation value by computing 10,000 times its value after having randomly shuffled the values for the other parameter. The significance level was estimated as the proportion of values in the permutation distribution that are higher than the value for non-permuted data. Fig. 2. Static CMC. Left side —CMC spectra for each participant (gray traces) and group-average (black trace). For each participant, the trace corresponds to the CMC measured at the sensor (among the 9 sensors overlying the left SM1 cortex) for which CMC was maximum in the 10–30-Hz range. Right side —Spatial distribution of the CMC averaged across participants and visualized using FieldTrip toolbox ( Oostenveld et al., 2011 ). CMC peaked above the left SM1 cortex. 3. Results 3.1. Static CMC Figure 2 (left side) presents the spectra of CMC in all participants. CMC maximal amplitude in the range 10–30 Hz was 0.058 ± 0.048; it was statistically significant in 16/17 participants ( p < 0.001) and marginally significant in the remaining participant ( p = 0.060). It should be noted that the peak frequency and to a larger extent the magnitude of CMC were variable across participants, in line with previous studies ( Ushiyama et al., 2011b ). Figure 2 (right side) shows the distribution of beta CMC in a representative participant. As expected, CMC was maximal in sensors above the left SM1 cortex. 3.2. Validity of force features Figure 3 presents a 500-ms excerpt of force signal and related force features. The force trace presents clear fluctuations of the order of 0.1 N. These fluctuations are well characterized by the four force features we selected. Clearly, the force features we extracted identify the short-lived events they were designed to capture. Figure 4 presents the auto- ( Fig. 4 A) and cross-correlation ( Fig. 4 B) of some of the force features. Force features were characterized by oscillations in their level of auto- and cross-correlations that essentially vanished within 200 ms of delays (| 𝜏|). The period of the oscillations in auto-correlation was about 70 ms for force increase and decrease, and 20 ms for force plateauing high and low. This indicates that CMC and power modulations should not display side peaks attributable to auto-correlation in force features. This is because the intrinsic temporal granularity (or temporal resolution) of CMC and power (about 200 ms since their computation was based on signals restricted to a 5-Hz band) is lower (i.e., higher duration) than that of the auto-correlation (20–70 ms). Beyond delays of 200 ms, correlation reached a non-zero baseline level 0.2–0.3. The fact that baseline correlation was not 0 is attributable to the non-normality of the force feature signals (Bishara and Hittner, 2015 ). Trivially, when the delay was 0 ( 𝜏= 0), auto-correlations were 1 for all force features. For such null delay, the cross-correlation was 0 between force increase and decrease, and between force plateauing high and low. Such a value was expected because the first signal of each of these two pairs was non-zero only when the other was at zero and vice versa. Beyond the initial peaks at 𝜏= 0, the level of auto- and crosscorrelation correlation was limited, with maximum values of about 0.1 (up to 0.3 for some participants) above baseline level. This indicates 5 S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 Fig. 3. Excerpt (500 ms) of force signal (A) and related force features (B–E) for a representative participant. Vertical arrows indicate the correspondence between force features seen in the raw force signal and peaks in the force feature signals. that force features were well suited to identify transient events that are not too interdependent. In summary, the limited interdependence of force features and their fast-decaying auto-correlation ensure adequacy for modulatory analysis of CMC and power in relation to transient events during continuous isometric contraction. 3.3. Physiological relevance of force features To determine the physiological relevance of the identified force features and validate the novel analysis approach, we first focused on the modulation of wide-band EMG power in relation to the force features. Figure 5 presents these modulations for each participant and averaged across participants. All participants displayed a significant modulation in wide-band EMG power in relation to all four force features ( p < 0.003). In relation to the force increase signal, EMG power showed a peak at –13 ms followed by a trough at 8 ms (see Fig. 5 for confidence intervals and power modulation values). The opposite response pattern was observed for the force decrease signal. In relation to force plateauing high, EMG power showed a trough at –2 ms; the reverse pattern was seen in response to force plateauing low. 3.4. Modulation of beta CMC, MEG SM1 power and EMG power in relation to force features Figure 6 presents the modulation of beta CMC, beta MEG SM1 power and beta EMG power in relation to the four investigated force features. These quantities were estimated based on 544 ± 31 s (mean ± SD across participants) of artifact-free, steady muscle contraction (no visual feedback) data. Overall, these three beta modulations displayed significant modulations in relation to the four force features. These modulations were statistically significant ( p < 0.05) in at least half of the participants (see Fig. 6 for exact numbers and Supplementary Fig. 1 for an illustration of the surrogate distributions used to assess statistical significance), except for the modulation of beta CMC in relation to force plateauing high , which was significant in 6 participants. The three measures showed a peak in relation to force increase and force decrease , and a trough in relation to force plateauing high and force plateauing low . All peaks or troughs were consistent across participants (i.e., in the same direction and with similar latency), so that they were significant at the group level ( p < 0.05 in all 12 instances; see Fig. 6 for exact p -values). In all instances, the latency of the modulation tended to precede the force features by ∼40 ms (see Fig. 6 for bootstrap confidence intervals). Some force features induced stronger modulations than others. Peaks of beta CMC were higher in relation to force increase than decrease (27.6% vs. 17.5%, t 16 = 2.37, p = 0.012); the same occurred also for beta MEG SM1 power (12.0% vs. 5.1%, t 16 = 3.66, p = 0.0021) and beta EMG power (10.7% vs. 5.6%, t 16 = 2.62, p = 0.018). Also, the suppression of beta CMC was more pronounced in relation to force plateauing low compared to high (–19.5% vs. –10.5%; t 16 = 3.46, p = 0.0067); the same occurred also for beta MEG SM1 power (–5.3% vs. –3.5%, t 16 = 3.69, p = 0.0020) and beta EMG power (–6.0% vs. –3.8%, t 16 = 3.87, p = 0.0014). Importantly, force plateauing low gave rise to stronger CMC depression in all participants showing significant modulation in relation to either of the force plateauing signals. Moreover, these significant ef- 6 S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 Fig. 4. Auto-correlation (A) and cross-correlation (B) of force features. Traces are in gray for every participant and in black for the group-average. Auto-correlation for force decrease and plateauing high were highly similar to those for force increase and plateauing low (respectively). Cross-correlation between force increase and plateauing low were highly similar to those for the 3 other pairs involving increase/decrease on the one hand and plateauing high/low on the other hand. Fig. 5. Modulation of wide-band (30–330 Hz) EMG power in relation to force features. From left to right, plots are for force increase, force decrease, force plateauing high, and force plateauing low. Traces are in gray for every participant and in black for the group-average. Power modulation is expressed in percentage of change relative to baseline (mean across 1 s < | 𝜏| < 2 s), as function of the time delay ( 𝜏) with respect to the considered force feature. Arrows indicate the location of peaks and troughs as well as a bootstrap confidence interval for their timing and modulation amplitude. fects evidenced by group-level comparisons were well mirrored by individual data. Among the participants showing significant modulations in response to either force increase or decrease , the modulation in relation to force increase was higher in 10/13 participants for CMC, 13/15 for beta MEG SM1 power, and 12/15 for beta EMG power. Among the participants showing significant modulations in relation to either force plateauing low or high , the modulation was stronger in relation to force plateauing low in 13/13 participants for CMC, 13/15 for beta MEG SM1 power, and 13/13 for beta EMG power. In light of this, further analyses will be conducted only on modulations in relation to force increase and force plateauing low . Given that beta CMC and power present modulations in the same direction (increases in relation to force increase/decrease; decreases in relation to force plateauing high/low), it is not possible to tell whether CMC modulations are attributable to changes in the signal-to-noise ratio in the MEG SM1 and EMG signals. In fact, the modulations we observed were related in multiple ways. Figure 7 shows the existence of strong associations between beta CMC and MEG SM1 power modulations ( Fig. 7 A), and between each of these two measures assessed in relation to force increase versus force plateauing low ( Fig. 7 B). 3.5. Relevance of identified beta modulations for force steadiness We next determine the relevance for force steadiness of the most salient modulations in beta CMC and MEG SM1 , in relation to force increase and force plateauing low . For that, we sought a linear association 7 S.J. Mongold, H. Piitulainen, T. Legrand et al. NeuroImage 261 (2022) 119491 Fig. 6. Modulation of beta CMC (A), beta MEG SM1 power (B) and beta EMG power (C) in relation to force features. Traces are in gray for every participant and in black for the group-average. The time intervals associated with significant modulation (p < 0.05 corrected for multiple comparisons; permutation statistics) are highlighted in red; exact p-values are provided in the top-right corners, along with the number of participants showing statistically significant modulation (p < 0.05; surrogate-data-based statistics). Bootstrap confidence intervals for the peak values are indicated with cyan shaded areas. Individual coherence traces (part A) were translated vertically so their baseline value (mean across 1 s < |t| < 2 s) aligns with that of the group-average. Power modulation is expressed in percentage of change relative to baseline. between individual values of CMC or power increase/decrease and force SD in different frequency ranges (0.5–3 Hz, 3–15 Hz and 15–30 Hz). The unique global permutation test revealed a significant association between modulations —in beta CMC and MEG SM1 power in relation to force increase and force plateauing low —and force SD at 0.5–3 Hz ( p = 0.043), 3–15 Hz ( p = 0.042) and 15–30 Hz ( p = 0.0024). Figure 8 presents the associations for force plateauing low . There was a significant negative correlation between force SD in all tested frequency ranges and the modulation in both beta CMC and MEG SM1 power. The association was especially salient with force SD at 15–30 Hz. Associations for the force increase signal were all positive but were significant only for the force SD at 15–30 Hz ( r = 0.58–0.62, p = 0.011–0.018) and not in the two other frequency ranges ( r = 0.31–0.40, p = 0.12–0.24). 3.6. Relation between identified beta modulations and global beta power depression We next asked if inter-individual variability in the amplitude of the observed dynamic modulation in beta CMC and MEG SM1 power relates to inter-individual variability in global depression in the beta MEG SM1 power during contraction compared with task-free power. Specifically, we aimed to determine if weaker dynamic modulations reflect a state where the beta sensorimotor rhythm is ( i ) unaffected by the contraction task or ( ii ) continuously attenuated. A unique global test revealed a non-significant trend of association between modulations —in beta CMC and MEG SM1 power in relation to force increase and force plateauing low —and the ratio between global beta power during isometric contraction and while not performing the task ( p = 0.083; non parametric permutation test). Figure 9 presents the associations for the modulations in relation to force plateauing low . Although the associations were not significant, they were more in line with the second alternative: the weaker the dynamic modulation (in beta CMC and MEG SM1 power), the deeper the global depression in beta MEG SM1 power. 4. Discussion In the present study, we examined how beta CMC, MEG SM1 power, and EMG power modulate in relation to different force features (rate of force change and plateaus) during submaximal isometric contractions. We found consistent temporal modulations that preceded changes in force features by ∼40 ms ( Fig. 6 ). The amplitude of these modulations, reflecting the extent of fluctuations in cortical involvement, was associated with force variability, which indicates behavioral relevance. These results are key to understanding the role of the SM1 cortex in dynamically maintaining steady muscle force. 8