Novel Metrics for High-Density sEMG Analysis in the Time–Space Domain During Sustained Isometric Contractions

Goal: This study introduces a novel approach to examine the temporal-spatial information derived from High-Density surface Electromyography (HD-sEMG). By integrating and adapting postural control parameters into a framework for the analysis of myoelectrical activity, new metrics to evaluate muscle fatigue progression were proposed, investigating their ability to predict endurance time. Methods: Nine subjects performed a fatiguing isometric contraction of the lumbar erector spinae. Topographical amplitude maps were generated from two HD-sEMG grids. Once identified the coordinates of the muscle activity, novel metrics for quantifying the muscle spatial distribution over time were calculated. Results: Spatial metrics showed significant differences from beginning to end of the contraction, highlighting their ability of characterizing the neuromuscular adaptations in presence of fatigue. Additionally, linear regression models revealed strong correlations between these spatial metrics and endurance time. Conclusions: These innovative metrics can characterize the spatial distribution of muscle activity and predict the time of task failure.

This observational study consisted of three separate sessions per participant each performing different types of contraction (i.e., isometric, pseudo-dynamic and dynamic).Each session lasted approximately 2 hours and sessions were spaced at least three days apart.The analysis presented in the main manuscript focuses on data collected from a single session, specifically on the fatiguing isometric contraction.The study was conducted at the Centre of Precision Rehabilitation for Spinal Pain (CPR Spine, University of Birmingham, UK) and in accordance with the Declaration of Helsinki.The study was approved by the Ethics Committee of the University of Birmingham (approval number: ERN_2022-0514).Written informed consent was obtained from all participants before participating in the experiment.For the precise capture of physiological indicators, our study employed the High-Density Surface Electromyography (HD-sEMG), while for torque measurements an isokinetic dynamometer, provided with a specialized chair attachment, was employed.This acquisition chain was thoughtfully designed to ensure a holistic evaluation of muscle performance and fatigue.The use of High-Density sEMG contributes to providing a detailed map of muscle activation patterns by analyzing the electrical activity generated by the muscle fibers during contractions.The increased spatial resolution allows for a complete understanding of muscle behavior under stress, making it an invaluable tool for our analysis.Concurrently, the use of the Biodex system, facilitated the precise measurement of torque output, which is essential for quantifying the physical force exerted by the muscles, and consequently thereby offering a direct indicator of muscle fatigue and endurance over time.

B. Participants
Nine healthy men (mean ± SD: age 29.1 ± 3.04 years, weight 75.1 ± 18.4 kg, height 178.0 ± 9.50 cm) with no history of low back pain, hip pain, or sciatica, participated in the study.Other exclusion criteria were rheumatic or neuromuscular disorders, history of chronic respiratory or neurological problems, spinal deformity or surgery, and cardiovascular conditions.

C. Protocol
The protocol consisted of a sequence of three Maximum Voluntary Contractions (MVCs), each lasting 5 seconds, with a 1-minute rest period between them.During these three MVCs, we recorded the Maximum Voluntary Torque (MVT) from which we extracted the peak value that served to establish a baseline before initiating the fatiguing contraction.Following five minutes of rest period to ensure the participants were not fatigued, they were instructed to perform a fatiguing contraction at 60% of their MVT peak, aiming to maintain this level for as long as possible.Visual feedback, based on the exerted torque, was provided to maintain the required level of torque at the targeted level.When the exerted torque level of the participant dropped to 50% of MVT peak or less (i.e., more than 10%) for a minimum of three consecutive seconds, task failure was reached, and the experiment was stopped.

D. Surface EMG and Torque acquisition
HD-sEMG signals were detected from the right lumbar Erector Spinae (ES) with two semi-disposable adhesive grids of 64 electrodes each (OTBioelettronica, Torino, Italy, model GR08MM1305).The grids are composed of 13 rows and 5 columns of electrodes (1-mm diameter, 8-mm Inter-Electrode Distance (IED) in both directions).Notably, the upper left corner does not contain an electrode; this absence is deliberate, serving as a reference point to establish the orientation of the grid.HD-sEMG signals were amplified by a factor of 150 (Quattrocento amplifier, OTBioelettronica, #3 dB bandwidth, 0.8-500 Hz), sampled at 2048 Hz and converted to digital form by a 16-bit A/D converter with a resolution of 152 μV (input range ±5 V).Before applying the grids, an investigator identified the L5 vertebra and marked the skin 2cm lateral to this point, extending until the ~T10 vertebra following the direction of the muscle fibers.Then, the skin was prepared by shaving the area (when necessary), cleaning it with an abrasive paste (SPES Medica, Italy) to reduce skin impedance and, finally, washing it with water and drying it.The two grids were applied one immediately above the other to cover and record signals from most of the lumbar ES muscle.A reference electrode was placed on the left posterior superior iliac spine (PSIS), as shown in Figure 1.The lumbar ES was chosen because it plays an important role in lifting activities.The right side of the muscle was selected based on the dominance side of the participants, who were all right-handed.Furthermore, a study by Varrecchia et al. [25] demonstrated no significant differences in bipolar and HD-sEMG parameters between the

Supplementary Materials
Novel Metrics for High-Density sEMG Analysis in the Time-Space Domain during Sustained Isometric Contractions Giovanni Corvini, Michail Arvanitidis, Deborah Falla, Silvia Conforto right and left sides of the lumbar ES during lifting tasks.This approach has also been effectively used in studies involving healthy individuals [22][31].
For the torque acquisition, an isokinetic dynamometer with the chair attachment (Biodex Corporation, Biodex System 3 Pro model, USA) was used.The torque signal was acquired during the back extension movements through a BNC cable attached to one of the auxiliary inputs of the sEMG amplifier.For this reason, being directly acquired by the sEMG acquisition software, the sampling frequency of the torque was equal to that of the sEMG signals, that is 2048 Hz.Before the beginning of each experiment, this device was calibrated.The torque signal was used to provide the participant with visual feedback of the exerted torque.The visual feedback, displayed on a screen placed in front of the participants (approximately 1.5 meters), showed the exerted torque, represented by a dynamic blue bar (increasing or decreasing according to the participant's muscle output), and the reference torque, represented by a constant red line at 60% of MVC.The participant should exert torque until the blue bar reached the red line and keep it constant as much as possible.

E. Data Processing
To calculate the Region of Activation (RoA) as the weighted centroid of the electrical activity greater than 70% of the local maximum, we follow this approach.First, we denoted (, , ) as the electrical activity at a point (, ) at time , where  and  are the spatial coordinates and  indicates the time instant.We refer to  and  as Fiber Transverse (FT) and Fiber Parallel (FP), respectively.At each time instant , we found the local maximum activity level   (), and we defined the active muscle Region (, , ) at time  as the set of points where the electrical activity is greater than 70% of the local maximum: (, , ) = (, , ) > 0.7 *   () (1) All the points falling below this threshold (0.7 *   ()) were set to zero.To calculate the weighted centroid along the two directions (FT and FP), we used the following formulas: where X and Y are the total number of columns (x) and rows (y) of the HD-sEMG grid, respectively.
The x-and y-coordinates of the original time series (RoA  0 , RoA  0 ) were smoothed with a low pass filter at 10Hz (3 rd order Butter) before deriving the spatial metrics.For the calculation of the metrics using the entire length of the signal, we used the discretized version of the RoA coordinates, referring to the time instant n (also called time point).From the entire length of the signal, we extracted three 10-second segments (S) at the beginning (T 1), the middle (T 2) and the end (T 3) of the contraction for the statistical analysis.In the metric definitions we will use K to refer to the total number of points within one of these segments S, to differentiate it from the total number of points N corresponding to the original length of the signal.

F. Metric Definitions
Using the formula from Prieto et al. [20], we derived and adapted several parameters typically calculated from the Centre of Pressure (CoP) to study postural control.In the following, we derived all the metrics from the RoA coordinates (RoA  0 , RoA  0 ), but the same definitions would apply to the coordinates of the Center of Gravity (CoG), defined as the weighted centroid of the total electrical activity.
First, we define the Activation Intensity Distance (AID), a vector representing the distance of each point of the muscle's Region of Activation along both directions (RoAFP, RoAFT) from a point representing the initial muscle activity.This initial centroid can be determined either at rest or at the beginning of a muscle contraction, depending on the objective of the study.In our case we aimed to assess changes relative to the beginning of muscle activity, therefore we arbitrarily chose to use as an initial centroid the CoG averaged in the first 20% of the duration of the contraction (that is different for each subject according to the total endurance time).This AID distance vector is computed for the entire length of the contraction at each time point n, and it is mathematically expressed as: where RoAFT and RoAFP are the demeaned coordinates of the RoA centroid, and N is the total number of points of the original signal length.These two coordinates are calculated according to the following equations: where RoA  0 is the original x-coordinate of the RoA (filtered) and CoG  ̅̅̅̅̅̅̅̅ is the mean of the x-coordinate extracted by the Centre of Gravity in the first 20% of the contraction.The same holds for the RoA  0 and CoG  ̅̅̅̅̅̅̅̅ being the y-coordinate.
Then, we define the Mean Activation Intensity Distance (Mean AID), which is the average distance of the RoA displacement from the initial centroid in the three selected time segments S. The Mean AID is computed according to: where K is the total number of points contained in a specific time segment S (in this paper S=3, each of 10 seconds).The AID metric could also be computed along each single direction (Transverse or Parallel) as: where AIDFT is the mean absolute value of the FT time series calculated in a selected time interval consisting of  points, representing the average FT Activation Intensity Distance (AID) from the initial centroid (i.e., the mean CoGFT in the first 20% of the contraction).The same applies to the Fiber-Parallel component.
The Total Activation Intensity Path (AIP) represents the entire path travelled by the RoA along both directions and it is defined as follows: This global measure can be calculated singularly for the FT coordinate (as well as for the FP coordinate by substituting the FT with FP).The total AIPFT (or AIPFP) indicate the total length of the distance travelled by the RoA along the FT (or FP) direction and it can be approximated by the sum of the distances between consecutive points extracted from the FT (or FP) temporal series, according to the formula: (10) where the two original RoA coordinates (RoA  0 ̂, RoA  0 ̂) were adjusted by removing the trend of the time series within the specific time window considered.For the definition of this metric, it was essential to use the original time series to avoid any potential biases that might arise from mean removal during calculations.However, the trend of the time series, as determined within each respective time window where the parameters were computed, was subtracted.This adjustment ensured that any differences between the RoA coordinates due to the magnitude did not influence the observed shifts in both directions, thereby preserving the integrity of the parameter calculations.
Using the Total AIP, it is possible to define and calculate the Mean Activity Intensity Velocity (AIV), which represents the average speed of the RoA over a specific time interval composed of  points.It is defined as: where K is the total number of points in the considered time interval for the parameter calculation.The AIP metric can also be computed on a single direction as follow: The Mean AIVFT represents the average velocity of the FT components of the RoA in a specific time interval composed of K points extracted from the total length of the time series.The same applies for the FP component.Then we defined two metrics able to identify the total space covered by the travelling RoA coordinates.Initially, we introduced the concept of Zone, an area enclosing most of the points of the AIP.We define the Activation Ellipse Zone (AEZ), which represents the area of the 95% bivariate confidence ellipse that shall enclose approximately 95% of the points of the RoA AIP.For this definition, the assumption is that that the number of data points used is large enough such that the ratio The formula for AEZ is given by: In this equation, F0.05[2, q-2] represents the F-statistic at the 95% confidence level for bivariate distribution with n points, and for large sample sizes (n > 120) this value is set to 3.00 [20].The   and   are the standard deviations of the FT and FP time series, respectively, while   is the covariance.By using the equation 5 and 6, these variables can be computed as: where   can be computed by substituting in equation 14 the RoAFT with RoAFP.Similarly, we also define the Activation Circle Zone (ACZ).This zone denotes the projected area of a circle whose radius is the 95% percentile of the AID; in simpler terms, this circle includes approximately 95% of the points of the RoA AIP.For this metric the assumption is that all Activation Intensity Distances follow a normal distribution.The formula to compute the ACZ is: In this equation the z0.05 refers to the z-score associated with a significance level of  = 0.05 in a standard normal distribution and it is equal to 1.645.The standard deviation of the AID can be computed as: Finally, we define the Activation Movement Zone (AMZ), which estimates the total region covered by the RoA AIP per unit of time.This metric is computed by summing the area of the triangles formed by two consecutive points on the AIP and the initial centroid.It can be thought of as being proportional to the combination of Mean AID and Mean AIV.The formula to compute the AMZ is: where  is the total number of points  in the selected time interval in which the metric is computed and, as for the Total AIP, we use the original RoA coordinates without the trend in the specific time interval .Finally, by using the standard deviations along the two directions FT and FP, we define the Activation Intensity Spread (AIS) as: This metric provides information on the amount of the RoA fluctuations around its initial position over time.A higher value indicates more variability in the muscle activation patterns.All these metrics have been normalized to the height of the subject.These metrics could be grouped according to the information they provide.Specifically, we have AID that represents the distance reached from the RoA within the HD-sEMG grid.It indicates the different points in the grid where Motor Units have been recruited relative to the center of activation at the beginning of the contraction.
Then we have the AIP and AIV, which are directly proportional, that express the cumulative distance travelled by the RoA and its velocity.These metrics do not provide information on the specific location reached but quantify the extent of fiber recruitment across different muscle regions.The AIV might provide indirect information on the rate of the selective fibers recruitment in different muscle regions.For example, if the entire muscle was completely active throughout the contraction, these values would be low, indicating more uniform activation of the muscle below the grid and suggesting a more static distribution of recruited motor units.
Finally, we have metrics related to the Zones in which they are active, and specifically ACZ, AEZ, AMZ and AIS.The first two metrics provide a global view of the RoA movements within the HD-sEMG grid.ACZ captures the overall shift in all directions, while AEZ highlights movement in specific directions, providing insight into the directional nature of muscle activation.
AMZ instead provides information on the oscillation of the RoA with respect to an initial point, and AIS quantifies the variability and spread of activation intensity.The latter, for example, could reveal asymmetries in muscle activation patterns that might be related to specific functional or pathological conditions.Together, these metrics might offer a robust toolset for analyzing muscle function, coordination, control, and adaptation in both research and clinical settings.
In addition to the calculation of these metrics, we also computed the parameters commonly analyzed in previous studies, which are the averaged value of the RoA coordinates in the specific time intervals, and the RoA shift calculated as the modulus of the RoA displacement in FT and FP directions at the beginning (1% endurance time) and at the end (99% endurance time) of the contraction [17].Additionally, we also determined the amount of the exerted torque (normalized with respect to the peak value extracted from the MVT), the Coefficient of Variation of the torque (CoV of torque), the envelope, (mean value across all channels of the grid) and the AIF (mean value across all channels of the grid).
Given the variability of the contraction duration among individuals, we will refer to different phases of the contraction in terms of their relative time percentages.For clarity, 0% and 100% correspond to the two "quasi-zero activity" window, 1% of the contraction time refers to its beginning, while 99% denotes its conclusion.In this study, based on the minimum total duration (around 70 seconds), we arbitrarily selected three 10-second intervals to capture the essential phases of the contraction.From these windows (i.e., initial 10 seconds (1%), middle 10 seconds (50%) and the last 10 seconds (99%), we extracted all the described metrics and performed the statistical analysis.

Figure 1
Figure 1 HD-sEMG grids applied on the right part of the lumbar Erector Spinae and the reference electrode.