Network of muscle fibers activation facilitates inter-muscular coordination, adapts to fatigue and reflects muscle function – Communications Biology

Subjects

Eighteen healthy young male adults (age 24.87 ± 2.14 years, height 174.65 ± 10.09 cm, and mean body mass 74.88 ± 8.34 kg) were recruited for this study. Participants were strictly recruited according to the following inclusion criteria: a) aged 20-30 years; b) BMI (in kg/m²) > 18.5 and < 30; c) normal physical activity > 5 and < 10 h/week, but without sport specialization (not active athletes), and d) blood pressure < 140/90 mmHg. The following exclusion criteria were applied: a) intake of prescribed drugs that could negatively affect muscle strength, such as corticosteroids; b) no current or previous injury that could prevent performance during the experimental protocol test, and c) any other physical condition (cardiac, respiratory etc.) that may have prevented the performance of a test protocol involving push-up exercise until exhaustion. The experimental protocol was approved by the local ethical committee (Clinical Research Ethics Committee of the Sports Administration of Catalonia), and was carried out according to the Helsinki Declaration. Before taking part in the study, participants read the study description and risks, and signed an informed consent. All ethical regulations were followed. To determine the sample size for this study, a power analysis was conducted using G-Power 3.1. Previous studies assessing fatigue effects on repeated exercise performed until exhaustion (same research design; 1 group, 2–3 repeated measures) have reported large effect sizes (0.8 and higher; Garcia-Retortillo et al., 2019). Thus, using an effect size of 0.8, α < 0.05, power (1 − β)  =  0.80, we estimated a minimum sample size = 15.

Study design and test protocol

After obtaining signed consent, participants visited the laboratory only one time. First, participants were be familiarized with the push-up movement — they practiced until they were able to execute the movement according to the protocol instructions (see study test protocol below). After at least a 5-minute rest, participants performed the study test protocol consisting of 3 sets of push-ups until exhaustion (Exercise 1, Exercise 2 and Exercise 3), interspersed by 5-min resting periods. The push-up sets will be performed according to the following instructions: 1) participants lied down (prone) with their knees on the floor and the hands on the preset marks (inter-acromial distance); 2) they did full elbow extension and place the arms perpendicular to the floor; 3) then, they lowered their bodies (straight line from shoulders to the knees) until the sternum was at 3 cm from the floor where a rope was placed; 4) when the rope was touched, they returned to standing position, and 5) repeated the push-up movement until exhaustion (Fig. 1a, b)^31,32,33. The push-ups pace was controlled by a metronome (Metronome Beats, Stonekick, UK), using a 3:3 tempo with three seconds down movement and 3 seconds up. Thus, one single push-up lasts 6 seconds. Each push-up set finished when participants were not able to do the next push-up or, alternatively, when they could not maintain the required 3:3 tempo. This means, that the observed changes in the intermuscular networks of muscle fiber interactions (measured by cross-frequency coupling) are induced by fatigue (both accumulated during an exercise bout, and residual across exercise bouts) provoked by the push-up movement (Figs. 3–5).

The repetition of three consecutive exercise bouts performed until exhaustion with short 5-min resting period between bouts allows to identify the effects of: (i) acute fatigue (accumulation of fatigue within exercise segments), and (ii) residual fatigue across Exercise 1, Exercise 2 and Exercise 3. Whereas acute fatigue occurs when the energy consumption exceeds the muscle aerobic capacity and a large fraction of the required energy has to come from anaerobic metabolism⁷² residual fatigue is characterized by neuromechanical and biochemical alterations (e.g., decrease in maximal force) provoked by previous exercise (i.e., Exercise 1 and Exercise 2)⁷³.

Two markers were utilized to assess the levels of fatigue and global physiological function in this manuscript. The first one is the amplitude of the EMG signal. Henneman’s Size Principle states that motor units within a motor unit pool are recruited in an orderly manner, based on their excitability, from the smallest to the largest unit^74,75. This means that slow type-I, low-force, fatigue-resistant muscle fibers are activated before fast type-II, high-force, less fatigue-resistant muscle fibers. As fatigue increases, more motor units need to be activated and recruited to maintain the same force output during push-ups. Consequently, surface electromyography (sEMG) amplitude increases as more motor units are recruited to maintain force output⁷⁶. As show in Fig. 1c, the amplitude of the EMG signal in both chest and arm muscles increases (i) within Exercise 1 (accumulation of fatigue) and across consecutive exercise bouts (residual fatigue; note the higher EMG amplitude at the beginning of Exercise 3 compared to Exercise 1). Same behavior is observed for the spectral power time series of individual muscle fiber types (frequency bands) where spectral power gradually increases with progression of the exercise bout — three-to-five-fold increase in spectral power amplitude comparing the Beginning (no fatigue) and the End (high fatigue) of the exercise bout (Fig. 1d).

Electromyography acquisition and EMG signal processing

Participants were asked to wear appropriate clothing that would not obstruct EMG electrode placement sites. Before the mounting of the EMG electrodes, participants’ skin was shaved and cleaned using alcohol, and left to dry for 60 s to reduce the myoelectrical impedance, according to the SENIAM guidelines⁷⁷. EMG signals from the following muscles were recorded simultaneously during the entire exercise protocol: left and right pectoralis major (Chest-Left and Chest-Right); left and right triceps brachii lateral head (Arm-Left and Arm-Right). The location of the surface electrodes (pre-gelled Ag/AgCl bipolar surface electrodes, DORMO LF-50) on each muscle was also carried out according to the recommendations of SENIAM organization. More specifically, pectoralis major electrodes were placed at over the center of the muscle belly along in the principal direction of muscle fibers of the pectoralis major (sternocostal part),  halfway between the sternal notch and anterior axillary line, approximately 2 cm apart in-line with muscle fibers^78,79, and the triceps brachii electrodes were placed at 50 % on the line between the posterior crista of the acromion and the olecranon at 2 finger widths medial to the line (Fig. 1a, b). After the electrodes were secured, a quality check was performed to ensure EMG signal validity. The aforementioned Chest and Arm muscles were selected since they present the highest myoelectrical activity during the push-up movement^38,39.

Data were recorded using a mDurance system (mDurance Solutions S.L., Granada, Spain), which includes a portable sEMG system (Shimmer3 EMG unit, Realtime Technologies LGd, Dublin, Ireland⁸⁰, and processed by means of Matlab (Mathworks, Natik, MA, USA). Raw data was recorded at a sample frequency of 1024 Hz and filtered online using a 5–250 Hz band-pass filter. A Notch filter was used with a width of 1 Hz at the frequency of 50 Hz (i.e., 49.5 – 50.5 Hz) to remove line interference from the grid. All EMG recordings were visually inspected and only signals without noise artifacts were utilized in the analysis.

Note that in this study we employ the same method to assess and quantify intermuscular interactions as in a previous work²⁴ (i) spectral decomposition (2.4.), (ii) cross-correlations between time series of spectral power (2.5.), (iii) Fourier phase randomization surrogate test (2.6.), as well as (iv) same data visualization method (cross-correlation matrices and intermuscular interaction networks; 2.7. and 2.8.).

Spectral decomposition

To uncover how distinct muscle fibers dynamically coordinate and synchronize their activation and integrate as a network across different muscles during the push-up movement, we first segment the EMG signals from the left and right pectoralis major (Chest-Left and Chest-Right) and the left and right triceps brachii (Arm-Left and Arm-Right) into 2-sec non-overlapping time windows with 1-sec overlap across the three consecutive Exercise sets. Within each 2-sec time window, we extract the spectral power S(f) from each EMG signal using the ‘pwelch’ function in Matlab, based on the discrete Fourier transform (DTF) and the Welch’s overlapped segment averaging estimator. For each time window we obtain a spectral power value in bins of 0.5 Hz for the range [5–250 Hz], that is, N = 500 is the number of spectral power data points for each window of 2 seconds. To probe specific contributions from different frequency bands F_i to the spectral power within each 2-sec time window of the EMG signal, we consider 10 frequency bands with equal width of 19.5 Hz. These frequency bands correspond to the range of activity of different types of muscle fibers in each Chest and Arm muscle, i.e., F1 = [5-24.5 Hz], F2 = [25-44.5 Hz], F3 = [55-74.5 Hz], F4 = [75-94.5 Hz], F5 = [95-114.5 Hz], F6 = [115-134.5 Hz], F7 = [135-154.5 Hz], F8 = [155-174.5 Hz], F9 = [175-194.5 Hz] and F10 = [195-214.5 Hz].

We calculate the sum of the power \(\widetilde{S}\left(f\right)\) across all frequency bins within each frequency band: \(\widetilde{S}\left(f\right):=\mathop{\sum }\nolimits_{i=1}^{n}S({f}_{i})\), where \({f}_{i}\) are all n = 39 frequencies considered in each frequency band F_i. Thus, we obtain 10 time series of EMG band power \(\widetilde{S}\left(f\right)\) with 1-sec resolution for each muscle during the three exercise sets, representing the dynamics of all representative EMG frequency bands. Frequencies below 40-60 Hz (corresponding to our frequency bands F1, F2 and F3) could be attributed to the activity of small alpha-motor neurons that innervate type-I slow muscle fibers. The frequency range 60-120 Hz (bands F4, F5, F6 and F7) could be attributed to medium alpha-motor neurons that innervate type IIa intermediate (oxidative) fibers, and high frequencies in the range 170-220 Hz (bands F8, F9 and F10), correspond to large alpha-motor neurons that connect to type IIb fast (glycolytic) fibers^{81,82,83,84,85}. The [45-54.5 Hz] range is filtered out by the Notch filter centered at 50 Hz to remove interference from electric grid, which modifies the EMG signal altering the spectral power of frequencies around 50 Hz. The obtained ten times series of EMG spectral power for each band F_i are then normalized to zero mean (µ = 0) and unit standard deviation (σ = 1). The obtained time series of EMG power \(\widetilde{S}\left(f\right)\) in each frequency band F_i reflects the micro-architecture (in 1-sec resolution) of synchronous modulation in the amplitude of muscle activation, and allows to track variations in coupling and network interactions of EMG frequency bands F_i across different muscles with accumulation of fatigue (Fig. 1).

Cross-correlations between time series of EMG spectral power in different frequency bands

To investigate cross-frequency interactions among muscle fibers (corresponding to EMG frequency bands), we consider all pairs of muscles involved in the push-up movement: (i) interactions of all frequency bands from same-type muscles (ChestL-ChestR and ArmL-ArmR) and (ii) different-type muscles (ChestL-ArmL, ChestL-ArmR, ChestR-ArmR and ChestR-ArmL). For each protocol set (Exercise 1, Exercise 2 and Exercise 3) and for each muscle pair, we calculate the bivariate equal-time Pearson’s cross-correlation for all pairs of time series representing EMG spectral power \(\widetilde{S}\left(f\right)\) in the frequency bands F_i where i = 1,…,10. This leads to 10×10 = 100 cross-correlation values C_i,j for each pair of muscles, as shown in the intermuscular cross-correlation matrices (Fig. 2a) — i.e., C_i,j quantifies the degree of coupling EMG frequency band F_i from one muscle with the frequency band F_j from another muscle. The cross-correlation values range from C_i,j = −1 (fully anti-correlated) to C_i,j = 1 (fully positively correlated), with C_i,j = 0 indicating the absence of linear relation between the power \(\widetilde{S}\left(f\right)\) time series of two EMG frequency bands. When cross-correlating the spectral power time series of two EMG frequency bands for an entire exercise bout, we found that the highest coefficient is the one corresponding to the zero lag (no delay). Our analyses of entire exercise bouts showed that the larger the delay in the cross-correlation, the lower the cross-correlation coefficient and thus, weaker interactions and network links.

Fourier phase randomization surrogate test and significance threshold for links strength in networks of intermuscular interactions

To demonstrate the physiological significance of the networks of intermuscular interactions obtained from our empirical analysis, we perform a Fourier phase randomization surrogate test on the EMG signals recorded from the four muscles (ChestL, ChestR, ArmL and ArmR). The test preserves the global spectral power in the different frequency bands F_i within the EMG recording but destroys Fourier phase information related to nonlinear EMG characteristics. Thus, the test eliminates heterogeneities in the fine temporal structure associated with short time scale modulations in the spectral dynamics of EMG frequency bands F_i within a muscle, which reflect activation patterns of different muscle fiber types. Note that in the Fourier-phase-randomized surrogate EMG signals the relative ratios among the average spectral power of muscle activation in the frequency bands F_i are preserved, while synchronous modulations in EMG frequency bands that underlie effective cross-frequency coupling and account for the nonlinear characteristics of EMG signals are eliminated. As a result, the degree of coupling between EMG frequency bands F_i from different muscles is remarkably reduced after the Fourier Phase Randomization procedure, since physiologically relevant information regarding coordinated activation of distinct frequency bands is lost.

Further, to test the statistical significance and physiological relevance of the observed hierarchical network organization and its reorganization with fatigue, we perform an additional step in our surrogate test to establish the significance threshold for network links strength. Specifically, for each network link, surrogates are generated considering signals from each combination of two randomly chosen subjects. Since we have 18 subjects in the database, 153 pairs of random subject combinations are generated. Each muscle pair involves 100 links within the corresponding subnetwork between 10 frequency bands F_i (muscle fibers) of each muscle. Thus, combining the 6 subnetworks representing all muscle pairs we obtain a distribution of 54,600 surrogate links (cross-correlation values) for each protocol segment (Exercise 1, Exercise 2 and Exercise 3)— i.e., 6 muscle pairs x 100 links x 91 subject combinations = 91,800 surrogate links. For each protocol segment distribution the mean µ_surr and standard deviation σ_surr are obtained. Thus, the significance threshold at 95% confidence level for the network links strength is defined as µ_surr + 2σ_surr for each protocol segment. We find that the significance threshold for network links strength is Th_exercise = 0.2 (corresponding to the highest Th value during the exercise segments. The physiological significance threshold is represented by horizontal dashed black line in Fig. 2b.

Intermuscular cross-correlation matrices

Group-averaged cross-correlation matrices represent pairwise coupling strength between the 10 frequency bands F_i of one muscle with the same bands derived from another muscle (i.e., 6 distinct muscle pairs: ChestL-ChestR, ChestL-ArmL, ChestL-ArmR, ChestR-ArmL, ChestR-ArmR and ArmL-ArmR) during a given protocol segment (Exercise 1, Exercise 2 and Exercise 3) (Fig. 2a). Matrix elements indicate the group-averaged (for all 18 subjects in the database) coupling strength between a frequency band that represents the activation of a specific muscle fiber type in a given muscle and a frequency band from another muscle. The dynamics of each muscle is represented by 10 EMG frequency bands with equal width of 19.5 Hz in the interval [5-215 Hz] (“Spectral decomposition”, Methods), and thus, each matrix consists of 100 elements (cross-correlation coefficients) representing interactions for each pair of muscles during a given Exercise set (Fig. 2a).

Intermuscular interaction networks

A multiplex network of subnetworks is obtained to visualize interactions from all pairs of muscles and their hierarchical organization within the network (Fig. 3a. 4a). We map the group-averaged matrices in Fig. 2a into different networks for the Beginning (first third) and End (last third) segments during Exercise 1 (Fig. 3a; accumulation of fatigue), and across Exercise 1, Exercise 2 and Exercise 3 (Figs. 4 and 5; residual fatigue). The graphical approach we employ is essential to identify universal patterns in the intermuscular network structure, the hierarchical organization of subnetworks and modules, and to track the transition in network characteristics with accumulation of fatigue. Each muscle is represented by a semicircle where color nodes represent distinct frequency bands F_i corresponding to different muscle fiber types in the muscle (“Spectral decomposition ”, Methods). Network links correspond to the values of cross-correlation matrix elements C_i,j in Fig. 2a and reflect the coupling strength between the frequency bands from two different muscles, where the frequency bands nodes and links for each pair of muscles form a separate subnetwork (Figs. 3, 4, 5). Links strength is marked by line color and width – we use distinct link color-code to demonstrate network reorganization with fatigue. To represent how network organization changes with accumulation of fatigue during Exercise 1 (Fig. 3) and with residual fatigue across repeated exercise bouts (Exercise 1, 2, 3; Figs. 4 and 5) we use the following links strength classification: weak links (0.20 < C_i,j < 0.35; very thin grey lines), intermediate links (0.35 < C_i,j < 0.5; thin green lines), strong links (0.5 < C_i,j < 0.65; dark blue thick lines) and very strong links (C_i,j > 0.65; magenta very thick lines). To illustrate the contrast among different subnetworks, we use a different links strength classification. For same-type muscle subnetworks (Figs. 7 and 8): weak links (0.4 < C_i,j < 0.5; very thin grey lines), intermediate links (0.5 < C_i,j < 0.6; thin green lines), strong links (0.6 < C_i,j < 0.7; dark blue thick lines) and very strong links (C_i,j > 0.7; magenta very thick lines). For different-type muscle subnetworks (Fig. 9): weak links (0.25 < C_i,j < 0.3; very thin grey lines), intermediate links (0.3 < C_i,j < 0.35; thin green lines), strong links (0.35 < C_i,j < 0.4; dark blue thick lines) and very strong links (C_i,j > 0.4; magenta very thick lines). Links corresponding to cross-correlation values C_ij < 0.35 (Figs. 3, 4, 5), C_ij < 0.4 (Figs. 7, 8), C_ij < 0.25 (Fig. 9) are not shown in the network maps (see Methods 2.6.).

To visualize the hierarchical organization of intermuscular network interactions among muscle fibers (frequency band F_i) from different muscles, the entire multiplex networks for the Beginning and End segments of Exercise 1 (Fig. 3a), as well as for Exercise 1, Exercise 2 and Exercise 3 (Fig. 4a) are presented separately as subnetwork maps for all pairs of muscles and distinct modules within the subnetworks. Each subnetwork map corresponds to a muscle pair for a given Exercise set and follows the same color-code as in the original network. We obtain two types of subnetworks: (i) same-type muscle subnetworks (ChestL-ChestR and ArmL-ArmR; Figs. 3b and 4b), and (ii) subnetworks of different muscle types (ChestL-ArmL, ChestL-ArmR, ChestR-ArmL and ChestR-ArmR; Fig. 5b). For all subnetworks, we show selected network modules for frequency bands F1, F3, F5, F8 and F10 that represent the activation of type-I, type IIa and type IIb muscle fibers.

To further dissect intermuscular interactions, we obtain stratification profiles of links strength for each subnetwork, and we track how the stratification profiles are modulated with accumulation of fatigue during Exercise 1 (Fig. 3) and with residual fatigue across repeated exercise bouts (Figs. 4 and 5). For each protocol segment/bout and all subnetwork representing different pairs of muscles, we calculate the average links strength for distinct modules within the subnetwork (corresponding group-averaged cross-correlation matrix in Fig. 2a). Considering that each subnetwork represents the interactions between EMG frequency bands for a given pair of muscles, we analyze the following network modules of cross-frequency coupling in each subnetwork: Low-low frequencies ([F1-F2]—[F1-F2]) bands interactions: average value of the coupling (matrix elements in Fig. 1d) between frequency bands [F1-F2] from one muscle with [F1-F2] from another muscle; Low-intermediate frequencies ([F1-F2]—[F3…F7]) bands interactions: average value of the coupling between frequency bands [F1-F2] from one muscle with [F3, F4, F5, F6, F7] from another muscle; Low-high frequencies ([F1-F2]—[F8…F10]) bands interactions: average value of the coupling between frequency bands [F1-F2] from one muscle with [F8, F9, F10] from another muscle; Intermediate-intermediate frequencies ([F3…F7]—[F3…F7]) bands interactions: average value of the coupling between frequency bands [F3, F4, F5, F6, F7] from one muscle with [F3, F4, F5, F6, F7] from another muscle; Intermediate-high frequencies ([F3…F7]—[F8…10]) bands interactions: average value of the coupling between frequency bands [F3, F4, F5, F6, F7] from one muscle with [F8, F9, F10] from another muscle; High-high frequencies ([F8…F10]—[F8…F10]) bands interactions: average value of the coupling between frequency bands [F8, F9, F10] from one muscle with [F8, F9, F10] from another muscle (Figs. 3c, 4c and 5b).

Statistics and reproducibility

Statistical analyses were performed using SPSS (version 23; SPSS, Inc.). We find a non-gaussian distribution for the cross-frequency coupling values obtained for the individual subjects in our database (see boxplots for the links strength in different network modules presented in Fig. 6). Thus, to assess effects of accumulation of fatigue during Exercise 1 (Beginning vs End segment), and residual fatigue across Exercise 1, Exercise 2 and Exercise 3 on average links strength for the individual network modules we perform a non-parametric Wilcoxon matched-pairs test. To account for multiple tests correction in Figs. 4, 5, 8 and 9 (Exc 1 vs Exc 2, and Exc 1 vs Exc 3), we have applied a Wilcoxon matched-pairs test with Bonferroni correction (alpha level 0.025).

Our study aims to identify hierarchical organization in the intermuscular network of muscle fiber interactions where subnetworks are composed of distinct modules. Thus, we quantify average links strength for the individual network modules, as represented by the bar charts in Figs. 3d, 4d, 5d, 7, 8, 9, where we compare two conditions (Beginning vs End segment, and Exercise 1 vs Exercise 3) to test for the effects of accumulated and residual fatigue. Correspondingly, we focus only on pair wise comparison where the average links strength of a given module at the Beginning is compared with the links strength of the same module under fatigue for the End segment of Exercise 1.

To verify the reproducibility of the experimental findings, we compared the data: (i) within each participant and (ii) across participants. The findings are consistent within participants and across participants. To promote reproducibility among other researchers, we utilized Matlab_2022a’s analytic tools to obtain spectral power distribution and cross-correlation values (see Methods for a detailed explanation).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.