Sustained fatigue assessment during isometric exercises with time-domain near infrared spectroscopy and surface electromyography signals.

The effect of sustained fatigue during an upper limb isometric exercise is presented to investigate a group of healthy subjects with simultaneous time-domain (TD) NIRS and surface electromyography (sEMG) recordings on the deltoid lateralis muscle. The aim of the work was to understand which TD-NIRS parameters can be used as descriptors for sustained muscular fatigue, focusing on the slow phase of this process and using median frequency (MF) computed from sEMG as gold standard measure. It was found that oxygen saturation and deoxy-hemoglobin are slightly better descriptors of sustained fatigue, than oxy-hemoglobin, since they showed a higher correlation with MF, while total-hemoglobin correlation with MF was lower.


Introduction
The phenomenon of muscular fatigue can be defined as a progressive decline of muscle performances [1]. It depends on different factors based on both muscle type and exercise typology. In general, the endurance time, during an exercise, is defined as the ability of an individual to maintain a contraction at a given torque until exhaustion [2]. The exercise-induced muscular fatigue, i.e. peripheral fatigue, is considered as a reversible loss of muscle force during work over time, which occurs as a safety mechanism to avoid structural damages to muscle cells [3].
The assessment of muscular fatigue is performed with the analysis of ATP metabolism (e.g. lactate, ammonia and oxipurines), oxidative stress (e.g. lipid and protein peroxidation biomarkers and others) and inflammatory biomarkers (e.g. leukocytes). In this framework, the best fatigue descriptor parameters to use are still under debate, because of their dependence on a wide range of experimental settings and on the range of population groups [3]. Furthermore, with these techniques, a biological sample is required and it is not possible to perform a continuous monitoring during the whole exercise. A series of non-invasive techniques for the monitoring of muscle fatigue are listed as an alternative to the invasive analysis: mechanomyography (MMG), surface electromyography (sEMG) [4], near-infrared spectroscopy (NIRS) and sonomyography (SMG) for both isometric and non-isometric contractions [5].
In many cases, the sEMG technique is employed. sEMG is a versatile and non-invasive method that allows the extraction of myoelectric properties of the neuromuscular activation associated with muscle contraction. These features and their variation provide insight into both biochemical and physiological changes in the investigated muscle and can be used for a real-time fatigue monitoring task [6]. In fact, the myoelectric signal displays time-dependent changes prior to any force modification, and this is the reason why it has the potential to both study contractile fatigue in isometric conditions, and predict fatigue onset. Past studies investigated the progression of fatigue resulting from isometric contractions using various sEMG time and frequency signal processing methods [7]. sEMG time-domain assessments such as the mean absolute value (MAV) and the root-mean-square (RMS) are rarely the only measures used for fatigue assessment: most research papers perform a combined evaluation of both signal amplitude estimates and spectral properties. In fact, it is commonly accepted that the sEMG spectral features show better performances than other-domain properties in assessing fatigue [8]. It was proven that sustained muscle contractions cause a frequency compression in the power spectral density (PSD) of the sEMG signal, long before the muscle arrives to exhaustion [9]. Besides a downward shift in the frequency spectrum, other changes have been found in both signal power (high-frequency decrease and low-frequency small increase) and spectrum slope (high-frequency increase and low-frequency decrease); such changes are considered as myoelectric manifestations of localized muscular fatigue [10][11][12][13][14][15]. For these reasons, the mean frequency (MNF) and median frequency (MF) are the most commonly used spectral variables and are considered as the gold standard for muscle fatigue assessment, especially in static contractile conditions, where the sEMG signal can be assumed as stationary [16]. In particular, a drop in the MF (or MNF) course and a simultaneous increase in sEMG amplitude reveal the onset of local muscular fatigue [7]. Though these spectral variables display similar behaviors, the estimation of MF is less affected by random high-frequency noise, and more affected by the fatigue phenomenon [17]. To this day, the effects of muscle force and muscle geometry on MF are subject dependent and the results across different papers not always consistent. However, it can be clearly assumed from the previous literature that MF is a golden-standard variable to identify muscle fatigue, particularly for isometric muscle contraction.
Since the fatigue process can be also characterized from the oxidative metabolism point of view, NIRS technique can be considered an interesting tool in the assessment of muscular performances during exercises [18]. Near infrared light (600-1000 nm) is selectively absorbed by oxy-hemoglobin (O 2 Hb) and deoxy-hemoglobin (HHb) allowing to evaluate total hemoglobin content (tHb = O 2 Hb + HHb), and the tissue oxygen saturation (SO 2 =O 2 Hb/tHb) changes in the working muscle. NIRS parameters provide information about the energy requirement rate, which should be maintained during the whole exercise, giving information about local skeletal muscle microvascular function and capacity to utilize oxygen [19,20]. The most spread NIRS modality is the continuous wave (CW), since the instrumentation has limited dimension and it is commercially available with moderate prices [21]. CW-NIRS was already employed successfully in the fatigue assessment alone or together with EMG [22]. However, this approach has some limitations in terms of depth selectivity and quantification of absolute hemoglobin content, which we have overcome with the time-domain (TD) NIRS modality. We found, in literature, few examples of use of TD-NIRS for muscle fatigue monitoring. Yamada et al. [23] monitored the vastus lateralis muscle during 50% maximum voluntary contraction (MVC) sustained isometric contractions, but in the results only relative variations for O 2 Hb, HHb and total hemoglobin (tHb) are presented. Neither absolute values nor SO 2 time courses are presented. Furthermore, the examples of application of TD-NIRS on muscle are limited, probably because of the absence of commercial devices, outside Japan, based on this technique. TD-NIRS was also applied during functional electrical stimulation (FES) on the calf muscle [24], on the vastus lateralis during incremental cycling exercises [25] and on the rectus femoris of post-stroke patients [26]. Some other examples are reported in order to study the effect of arterial occlusion [27] or cutaneous vasodilation [28].
In this paper, we study the combined use of TD-NIRS and sEMG in muscle fatigue assessment. When a muscular exercise is performed, the hemodynamic response as measured by NIRS can be roughly divided in two parts: the initial one is a "fast phase", where there is a rapid change in the HHb and O 2 Hb parameters which represent the arising of the early fatigue. Then, a "slow phase" can be identified, which reflects the local energy turnover during isometric contractions [29]. This behavior will be identified by sEMG and TD-NIRS parameters. In particular, we will investigate the possible correlations between sEMG and TD-NIRS parameters during the "slow" phase, simultaneously measured, which, to our knowledge, were never presented during sustained isometric exercises.

Devices and probes
sEMG acquisitions were performed with a commercial device (Cometa, Italy) featuring 16 channels wireless electromyography. During the recordings, one sEMG probe was positioned in the belly of the deltoid lateralis muscle of the dominant arm of each healthy volunteer with Ag-Ag/Cl electrodes.
To perform TD-NIRS acquisitions, the class IIa medical device, previously developed by the authors at the Department of Physics of Politecnico di Milano, was employed. It has two diode lasers operating at 688 nm and 828 nm as light sources and two hybrid photomultipliers tubes as detectors. Laser pulses are delivered to the tissue by means of optical fibers according to the space multiplexing modality [30]. The electronic acquisition chain is based on time correlated single photon counting (TCSPC) technique. For further details, the reader can refer to the paper by Re et al. [31]. The optical bundle for the detection is 90°degrees bended, while the injection beam was bended by means of a glass prism. Both the fibers were host in a 3D custom probe printed by means of a filament printer (Sharebot NG, Sharebot s.r.l., Nibionno, Italy) and a black PLA filament (3Ditaly, Roma, Italy) compatible with diffuse optics measurements [32]. We then arranged one injection and one detection channel with an interfiber distance of 30 mm.
The sEMG electrodes (one channel) were placed perpendicular to the TD-NIRS probe, but on the same anatomical landmark, so that simultaneous acquisition on the same muscle volume were performed. In Fig. 1, (a) picture of the arrangement of the TD-NIRS optodes (one channel) and sEMG electrodes is presented. To fix the probes to the arm, we employed a black elastic bandage that guarantees good adhesion to the skin without compressing the muscle. It also acts as a shield from the ambient light, avoiding its entrance into the optical fibers. sEMG and TD-NIRS devices were also synchronized with a TTL signal sent by the sEMG (acting as master) to the TD-NIRS instrument (acting as slave).

Subjects and protocol
We recruited 12 healthy subjects (28.9 ± 3.2 age), 6 males and 6 females. All the enrolled subjects had body mass index (BMI) ranging from 20 to 25, within normal weight range. The study was approved by the Ethical Committee of Politecnico di Milano and was conducted in compliance with the Declaration of Helsinki. Before the study each subject was informed about the experiment modality and gave a written consent. All the subjects had no neurological or orthopaedical diseases affecting their performance and were not trained in fatigue tasks.
During the measurements, subjects were sitting on a chair holding their dominant arm elevated laterally at the shoulder level in order to keep the arm horizontal. The positioning was guaranteed by placing a mark on a vertical support, which the subject had to reach (without touching it) with the hand. The protocol was conceived to elicit fatigue during an upper-limb isometric trial. The position was held until exhaustion, or for a maximum time of 330 seconds (task period). The subject was verbally motivated to maintain the position. The task period was preceded by 30 s of baseline and followed by 210 s of recovery, during which the arm was lying along the body. For TD-NIRS, the instrument response function (IRF) with a "reference" modality [33] was acquired for each subject.
In order to speed up the onset of the fatigue, subjects held a water bottle full of water for the whole experiment. The amount of water in the bottle was determined in order to reach a selected static torque at the shoulder (M w ). The torque to be added to the physiological one was proportional to the static maximum shoulder torque (M s ), which was computed with a biomechanical model of the upper limb. M s depended on the arm mass, forearm mass and upper-limb barycenter locations, computed with anthropometric tables starting from body mass and height ( [34]), according to the following Eq. (1): Being M s computed as follows Eq. (2): where m b is the mass of the bottle, m w is the mass of the added water, b is the forearm length, m a is the mass of the arm, m f is the mass of the forearm, g is gravity acceleration, a l is the shoulder-arm segment barycenter length, a is the arm length and a f is the elbow-forearm barycenter segment length. The overall torque to be held was M w = M s (1+k), with k>0 that in this study was set at 0.5. We calculated the customized amount of water m W for each subject as Eq. (3): We explicitly note that m w depends on the anthropometric data of each subject.

sEMG analysis
Raw sEMG data were sampled at 2000 Hz. The pre-processing was performed in Matlab (Mathworks Inc., Natick, USA). The raw signal was notch-filtered to erase the 50 Hz frequency interference, then it was bandpass filtered with a 6 th order Butterworth filter with cut-off frequencies 5 and 500 Hz. Since the lifting arm phase lasted about 1.5-2.5 s, we decided to use sEMG signal starting from 35 s until the end of the lifting phase to eliminate the transient period related to arm elevation. Then, for each trial, the sEMG was divided into adjacent consecutive windows of 5 s -starting from the time instant 35 s. The number of windows depended on the subject's capability to complete the trial; the maximum number of windows was 65 in case the subject could hold the weight in isometric contraction until the end (we removed the last window to avoid the transient when the limb was brought back to the resting position). We stopped the windowing when the trial was completed, or when the subject dropped off. In each window, we computed the average RMS value [15,35] and the MF [17,36]. Then, for each subject, we created histories by plotting the RMS and MF as a function of time-windows progressing in time. In order to compare performances across subjects, we normalized both RMS and MF histories to their maximum value within each subject, obtaining adimensional, normalized RMS and MF values ranging from 0 to 1.
To perform frequency analysis, we preliminary tested sEMG signal's stationarity within each window with the augmented Dickey-Fuller (ADF) and the reverse arrangement (RA) tests [37,38].

TD-NIRS analysis
The acquired TD-NIRS reflectance curves at each wavelength were fitted with the solution of the diffusion equation for a semi-infinite homogenous medium after the convolution with the IRF [39]. From the absolute values of the absorption coefficient, we calculated O 2 HB and HHb with the Beer's law, knowing the extinction coefficients at the given wavelengths [40]. Starting from these values, we calculated the tHb and SO 2 . These values were calculated for all the temporal points of the acquisition. Furthermore, for the initial baseline, we considered the average value between 5 and 25, interval where no movements were performed. We checked the quality of the data acquired in terms of detected photons number (> 350 kcounts/second during the whole acquisition) and fitting parameter (χ 2 ≤ 2). We also check the stability of the instrument through the acquisitions on the different subjects, evaluating the fluctuations of the IRF in terms of curves barycenter.
For the analysis of correlations between TD-NIRS and sEMG signals, we removed from the TD-NIRS signal all the parts where the sEMG indicated an absence of movement. This procedure was personalized for each subject, since not all the participants were able to elevate the arm for the whole allowed time (330 seconds). Then we kept TD-NIRS signal starting from 35 s until the end of lifting phase. Furthermore, according to what was done for the sEMG signal, we divided also the time courses of the TD-NIRS parameters into adjacent consecutive windows of 5 s. For each window we calculated the average value and the standard deviation of the TD-NIRS parameters. With these procedures, TD-NIRS and sEMG data for each subject were not only synchronized but also of equal length.

Coupling TD-NIRS and sEMG analysis
To provide cross-domain correlations on signals, some segmentation criteria were needed, in order to clearly identify the slow and fast phase in TD-NIRS parameters. In this study, we executed a two-step procedure to provide accurate automated segmentation of the NIRS signal. First, as a compromise between smoothing and avoiding altering the original signal, we implemented a 5 th order Butterworth low pass filter with cut-off frequency f c = 0.3*f NYQUIST (since this application is novel, f c was empirically determined). Then, we found the local minima and maxima of the TD-NIRS signals; this was achieved implementing a first derivative algorithm on the filtered signals. The first local minimum/maximum was consequently appointed as reference cutpoint to separate the fast phase from the slow phase of the TD-NIRS signal. This allowed us to separate the initial transient of the signal, and to consider the rest as slow phase. Of course, the cutpoints could vary from subject to subject, but also among TD-NIRS parameters within the same subject.
We then performed a correlation analysis using the Pearson linear correlation coefficient between all the possible TD-NIRS and sEMG measure pairs. This analysis was performed using different "driver signals". The driver signal for an analysis was defined as the one currently controlling the segmentation of the epochs. The purpose of using driver signals was to investigate if different segmentations could lead to higher/lower correlation levels with the RMS average value and MF. The role of driver signal was assigned, alternatively, to three out of four TD-NIRS variables: O 2 Hb, HHb and SO 2 . In this work, we used the following notation: SO 2 →HHb(MF) means that we are correlating the HHb with MF with the segmentation done with SO 2 as driver signal. Since in the results we will motivate the use of MF only, as well as not using tHb as driver, we omit MF and simply write "SO 2 →HHb".

Statistical analysis
TD-NIRS and sEMG correlation data were merged in a 12-by-12 matrix where each column contained the MF-NIRS correlation values calculated for each subject (12 subjects, each one with 12 correlations -3 driver signals x 4 MF-NIRS i combinations). In this way, it was possible to compare the effect of different segmentations on correlations and to search for statistically significant differences in our dataset. Statistical analysis was conducted with Matlab (MathWorks Inc., Natick, USA). To verify which TD-NIRS parameter would correlate better with the MF, we removed the outliers from each dataset (rmoutliers Matlab built-in function) and then proceeded with the computation of the absolute value for each correlation to allow comparisons. We verified that data were now normally distributed and ran 3 separated one-way ANOVA test (one for each driver signal) with correlations between MF-NIRS data as factor. The ANOVA test compares the means of different independent groups to find statistically significant differences between the variables. In all tests, the significance level was set at α = 0.05.

EMG parameters
From a visual inspection of the EMG parameters, we found that MF (reported in Fig. 2), was mostly monotonous and with a fairly intra-subject constant decay rate. However, for some subjects, we noted a decrease of the decay rate in the last part of the trial, when close to the end of the considered time window for detecting fatigue. Inter-subject average slope could vary. The average slope values are reported in detail in Table 1, together with the task time and the number of temporal windows. On the other hand, in our dataset we could not find a common trend in the RMS time-course, as shown in Fig. 2. In most cases, the RMS showed no clear variations, or was not monotonic; instead, some subjects showed increased RMS values towards the end of the recording when they were close to exhaustion, or vice-versa.
The ADF and the RA tests revealed that the hypothesis of stationarity could not be rejected on the chosen windows.

TD-NIRS parameters
The TD-NIRS signal quality check confirmed that the minimum count rate was reached during the whole acquisition for all the subjects. The χ 2 value of the inversion procedure was always ≤ 2 except for some few seconds for subject 11 (anyway it was < 2.8) because of some direct light, which we verified didn't affect the hemodynamic time courses. The standard deviations of the IRF barycenter positions of each subject, was always <10 ps.
In Table 2, TD-NIRS parameters obtained for the initial baseline (standard deviation negligible, not reported in table and in the graphs neither) are reported. In Fig. 3(a) the complete time courses for the absolute values of O 2 Hb (red) and HHb (blue) are shown for a representative subject. During the baseline period (0-30 s) the signals are almost flat for each subject. When the subject starts to move the arm in order to displace it to the horizontal position (starting around at time instant 30 s), we observe that the hemodynamic parameters undergo a fast variation with an increase for O 2 Hb, and a decrease for HHb. In the inset, the detail of these hemodynamic parameters' behavior is highlighted. After few seconds, the trends swap. Since the aim of our work is to consider the oxygen metabolism during a sustained contraction, this initial transient was then removed, as described in Section 2.4, together with the baseline period. In Fig. 3(b), we show the time course as obtained in this way. We used a new time scale, so that now the origin of the time is set 35 s after the start of the acquisition. Here we can clearly notice a rapid change of HHb (increase) and O 2 Hb (decrease). After 30 seconds (65 from the beginning of the acquisition), the two parameters start to settle around a plateau value; then, they show again changing trends but slower. The final part of these signals is more variable among the subjects: signals do not show always a plateau, but sometimes they start to vary again (increase for O 2 Hb and decrease for HHb) but with slower slopes. At the end of the exercise, observing again Fig. 3(a), we can find the typical hyperemic phase of the recovery period. The pattern of SO 2 of a typical subject recall the one for O 2 Hb, while the behavior of tHb is not uniform among the subjects. In some cases, it is almost flat, in some others it has a significant increase.

sEMG and TD-NIRS Assessments
As a first result, we portray in Fig. 4 the typical graph of the absolute values of the TD-NIRS variables together with the MF achieved with the sEMG. In this figure, all the signals are reported according to the 5 s windows. In the TD-NIRS data, the standard deviations are also shown. We underline two main results that can be noticed among all the subjects. First, positive correlation between MF and HHb were found; negative correlation between MF and SO 2 and O 2 Hb were also spotted in most of the subjects; while no clear correlation between MF and tHb was detected. This pre-processing procedure also revealed that, in several cases, tHb had no transient phase: either the signal was mostly monotonic with many small fluctuations, or the trend changes (concavity change) occurred at a late time sample. For this reason, tHb was not used as driver signal.
Moreover, tHb showed the lowest correlation with the other TD-NIRS variables, and with a high variability (with O 2 Hb: 0.715 ± 0.397; with HHb: 0.062 ± 0.528; with SO 2 : 0.411 ± 0.493). These preliminary results, along with the lack of transient phase in many subjects, justify our choice not to use tHb as a driver signal. Following these choices, we could provide 12 NIRS-sEMG correlations: 4 correlations using O 2 Hb as driver, 4 using HHb and 4 using SO 2 . Furthermore, since RMS showed poor correlation and varying trend with both NIRS measures and the MF, it was excluded from further analysis procedures being discarded a-priori as a possible reliable descriptor of fatigue.
From the preliminary correlation test on not-segmented signals (detailed results are reported in Fig. 5), we can see that three out of four TD-NIRS measures showed high average correlation values between themselves, as it can be seen on histogram (a), (b) and (c). The O 2 Hb->SO 2 correlation was the one with the lowest standard deviation value. The correlation between all the TD-NIRS and sEMG parameters of the dataset (considering the not-yet segmented signals) is reported on Table 3.  As for the sEMG domain, the mean correlation value found between the RMS and the TD-NIRS variables was 0.079 ± 0.549, while the RMS-MF correlation value was -0.112 ± 0.623. Both the poor cross-domain correlation and the high variability of this data caused for the exclusion of the RMS from further analyses in this study.
On the contrary, the MF showed good levels of cross-domain correlation, albeit with rather high standard deviation values.
Afterwards, NIRS data were segmented -using O 2 Hb, HHb and SO 2 as drivers. In Fig. 6, we portray graphically the results of the segmentation for one subject in our dataset. Once these three sets of cutpoints were found for all subjects, the cross-domain correlation test was repeated only on the slow phase, and again using Pearson's linear correlation coefficient.

TD-NIRS and sEMG correlations
Given the high levels of standard deviation due to outliers, we removed outliers from each distribution and then re-run the correlation test. As can be seen on the bar graphs in Fig. 7, this further step produced a notable reduction of the standard deviation values and increased the average MF-NIRS correlation values. The new values are reported in Table 4.   Three one-way ANOVA were performed, one for each segmentation. Results are reported in Fig. 8 as three different boxplot sets. All tests revealed statistically significant differences. For the O 2 Hb-driven dataset (p value = 0.008); for the HHb-driven dataset, (p value = 0.016); for the SO 2 -driven dataset (p value = 0.006). Fig. 8. Results of the ANOVA test on the normalized MF-NIRS correlation of the dataset without outliers (from left to right: driver signal is O 2 Hb with p-value = 0.08; driver signal is HHb with p value = 0.016; driver signal is SO 2 with p value = 0.006).
From the Tukey-Kramer HSD post hoc test applied to ANOVA stats, statistically significant differences were found when considering the following comparisons. The comparison between the MF-tHb Pearson coefficient was always significantly different from both the MF-HHb and MF-SO 2 correlation values (for the O 2 Hb-driven dataset, p value = 0.037 and p value = 0.008; for the HHb-driven dataset, p value = 0.025 and p value = 0.034; for the SO 2 -driven dataset, p value = 0.018 and p value = 0.009).
As can be seen in all boxplots reported in Fig. 8, we could conclude that the worst-performing measure across the different segmentations is the MF-tHb correlation. Instead, whatever driver signal was chosen, SO 2 and HHb performed averagely better than tHb, even if they could not be distinguished from O 2 Hb. Thus, we could conclude that 3 out of 4 TD-NIRS parameters correlate highly with the MF during the slow phase.

TD-NIRS signal
The TD-NIRS quality check determined that no subjects or parts of the time courses had to be eliminated because of direct light or inadequate photons count rate. Furthermore, the standard deviations referred to the IRF barycenter show a good device stability.
The subjects' baseline values of the hemodynamic parameters (see Table 2), are not uniformly distributed around an average value, except for the SO 2 , which is a good candidate to represent the metabolic state of resting muscle. Our testing scenario involved healthy subjects and analysis on the deltoid lateralis muscle. Even though shoulder is implied in most of our movements, very few papers report examples of muscular fatigue studies with EMG and NIRS on this muscle. Ferguson et al. [41] and Jensen et al. [42] investigated deltoid and trapezius during a repetitive task; while Lusina et al. [43] investigated deltoid deoxygenation during arm cranking exercises. Unfortunately, no absolute values for its hemodynamic parameters are presented, but just their relative variations with respect to a baseline period. We find just an example, where absolute values for SO 2 of the deltoid muscle is reported as 85.4% ±4.4% [42]. Hamaoka [28]. Both these studies employed TD-NIRS and used homogenous models for the data analysis. A direct comparison between the hemodynamic parameters of different muscles compartments, might not be the best approach, since differences in oxygen status are reported in literature [44]. The values we find for SO 2 are almost homogeneous and in line to the expected one. The values for O 2 Hb and HHb respect the ratio of 2/3 that they should have in human tissues [45] but are not uniform. The heterogeneity could be due to the homogeneous model we have employed for analyze TD-NIRS data, which doesn't consider that the volume under the optical probe is a multilayered medium. In particular, we should consider the adipose tissue thickness (ATT), which is different among subjects and can create confounding effects in the estimation of the basal optical properties. Furthermore, there are evidences of differences in SO 2 between male and female in literature, always because of the different ATT [23]. Other approaches to analyze TD-NIRS muscle data were already proposed [46,47], but for the aim of this paper, the determination of the absolute values during the baseline period is not a critical point, i.e. we do not consider the baseline for the statistical analysis. Since the typical distance from the skin surface to the deltoid muscle is between 6 and 8 mm [42], we can consider the ATT thin enough to consider the sample as homogenous.

sEMG signal
In this study, we found that MF presented a monotonic decrease in all subjects. This result matches with previous findings in the literature and allowed us to consider MF as a reliable biomarker for fatigue. In our experiment, MF had a monotonic decrease that could be observed even if with different magnitude between subjects. In fact, during a sustained isometric contraction, the sEMG power spectrum shifts to lower frequency bands as local muscle fatigue develops. This compression is linked to variations in both spectral variables and conduction velocity (CV) [15], and may also be correlated with the type of motor unit (MU) recruited [48,49]. In fact, the typical MF decrease associated with fatigue onset has been attributed to the decrease in CV [50,51] meaning these two variables are highly correlated, and MU synchronization [52]. It has been shown that during a sustained submaximal isometric contraction, the decline rate of MF is highly correlated with both the endurance time and increased levels of metabolites implicated in the development of local muscle fatigue [53][54][55][56][57]. The decay rate of the MF and endurance time were found to be linearly correlated, and thus the first one was used to predict the latter, regardless of the force of isometric contraction being sustained. Furthermore, changes in MF are associated with the high frequency fatigue of a muscle and is highly correlated with a decline in force from the fresh state [9,57]. Interestingly, we found a MF slope of about 0.10-0.30%/s which is in line with the one found in previous works [58].
On the contrary, throughout literature on muscle fatigue, the RMS is often quite variable in terms of time-series fatigue-related patterns, and our data confirmed this unclear trend. The variability of RMS often leads to inconsistency in muscle fatigue evaluation. This is because the EMG signal amplitude can be easily influenced by experimental conditions (workload, contraction type, endurance time), and thus the use of this variable as fatigue index should be interpreted with caution [59,60]. Some studies based on isometric exercises showed that, in some cases, an increasing trend was found on the EMG signal amplitude [61,62], as it happened also in this study for some subjects. These higher values were associated with MU synchronization and increase of the firing rate, which counteracts the fatigue process of the MUs and contraction force decrements [63].

TD-NIRS versus sEMG
In Fig. 3(b), we can observe the typical response, when a muscular sustained isometric exercise is performed. From a physiological point of view, it can be divided into two main parts. The initial one (fast phase), where there is a rapid change of HHb (increase) and O 2 Hb (decrease). This is due to the sudden energy demand at the onset of the exercise together with the increase in intramuscular pressure (IMP) [29]. In particular, the oxygen delivery does not match the oxygen demand even if the oxygen extraction increases. When static contractions are performed, the vascular bed is more compressed resulting in a temporary occlusion of the blood flow. Following this initial fast phase, as we can see in the figure after about 30 s, there is a second one (slow phase), which reflects the local energy turnover during isometric contractions. During the fast phase, the principal mechanism for energy production is the aerobic one (oxygen is consumed). During this phase, the type I fibers, the oxidative ones, are mostly recruited. In the slower phase, there is less oxygen available (decrease oxygenation) and higher occlusion due to the IMP. This lack of oxygen leads to the preference for anaerobic metabolism [64] and mostly the recruitment of the type II muscular fibers. During the fast phase, we observe a drop in SO 2 , which reflects the fact that large oxidative fibers are prevalent [65]. The deltoid muscle is about 33% type I fibers, i.e. it has a reduced aerobic capacity with respect to other muscles such as trapezius [66]. We acknowledge that the description of the mechanisms underlying muscular metabolism is more complex than the one presented above. In particular, in literature, muscular fatigue is described by the phosphocreatine system in the first 5-10 seconds, followed by glycogen and lactic acid conversion in the next couple of minutes, followed by oxidative (aerobic) metabolism for longer fatigue durations [67]. Since the aim of this paper is focused on the technical aspect of comparing two different signals, rather than going deep into the physiology of the muscular fatigue, we preferred to divide the timeline in two main blocks (fast and slow) as suggested from a preliminary analysis. Future works will investigate in TD-NIRS time-courses, other inflection points or significant changes, which could better describe all the phases of the muscular metabolism during the fatigue process.
For what concerns the interpretation of the sEMG signal, we can affirm that MUs are recruited in an orderly fashion, from type I, oxidative, to type IIb, mostly anaerobic, as the force output increases [65]. During the fast phase we should observe correlations between the sEMG and TD-NIRS parameters because higher energy requirements should lead to a faster fatigue process [29], in particular between the parameters' slopes. This behavior was already widely investigated in literature. Guo et al. [68], for example, tested the effect of the short term fatigue (80 s), on the extensor digitorum, and showed the effect of muscular fatigue both in NIRS (O 2 Hb) and sEMG (MF) parameters. In the study by Taelman et al. [69], during a static elbow flexion, a strong correlation between frequency content of the sEMG signal and TOI was found. Very recently, also study on elderly people were presented [70]. We also observed that correlations were noticed between NIRS and sEMG parameters during the fast phase, and used as predictors of medium-term fatigue.
One of the novelties of this work was to investigate deeply sustained fatigue, focusing on the slow phase. In literature, different objective criteria and ad-hoc algorithms were proposed to separate these two phases -focusing mainly on the determination of the initial transitory phase and, eventually, on how to correlate its slope with the slope of the MF [29,65]. Most of the examples of sEMG-NIRS combined studies are proposed for the "fast phase", while we focused on the slow one. We assessed the slow phase for some minutes of isometric contraction, when the TD-NIRS parameters reach a plateau condition. After the plateau we found different behaviors that we characterized with the analysis proposed in paragraph 3.3 and 3.4 in order to couple TD-NIRS with sEMG. During the slow phase, the MF decreases in a monotonic way.
As already depicted in the introduction, we could not find works where all the parameters that characterize the muscular oxidative metabolism are interpret together. For this reason, also the interpretation of the muscular fatigue and the identification of a unique NIRS representative parameter is nowadays controversial. Thanks to the fact that with TD-NIRS, it is possible, with one single measurement to extrapolate all the four hemodynamic parameters (O 2 Hb, HHb, SO 2 and tHb), we could show correlations between their absolute values and the sEMG signal in the same paper, for the first time. The main finding of this analysis is that there are 3 out of 4 NIRS parameters that show high correlation with the MF signal during the slow phase, and could potentially be used as biomarkers for muscular fatigue. In particular, our tests showed that SO 2 and HHb are a better biomarker than tHb, even if they cannot be distinguished from O 2 Hb. This result should not surprise since we also showed that TD-NIRS parameters are strongly correlated one with the other (especially O 2 Hb, HHb and SO 2 ). We also found that this result is not achieved on all the enrolled subjects and that there might be some outliers that differ from this trend (especially when considering HHb). A more detailed investigation on variability was not reported in this paper due to the relatively small number of subjects.
In our study, we also could exclude some sEMG and TD-NIRS parameters to be used as biomarkers for sustained fatigue. RMS showed high variability among subjects, not revealing any average relevant correlation with MF or TD-NIRS parameters. In the literature, there are contrasting evidences suggesting that probably RMS has inferior utility compared to other indices, indicating inconsistent performance to evaluate muscle fatigue [60]. This is due to some factors related to experimental conditions such as muscle contraction, workload, endurance time, and other factors, the use of this index as fatigue indicator should be interpreted with caution [60,71,72]. Coherently, in our study, we found high variability across subjects. In NIRS, while 3 out of 4 parameters correlated well with MF, tHb showed poorer performances. This can be explained considering that tHb reflects the blood volume under the sample. During the contraction, because of the IMP, it is possible to observe an occlusion of the vessels, which is not always complete also if the contraction is high (around the 80% of the MVC) [65]. Also Yamada et al., in a TD-NIRS study on the vastus lateralis during sustained isometric contraction of 1 minute, show different patterns for tHb for the fast phase, since the 50% of the MVC would not cause a complete suppression of blood flow [23]. The torque we added with the employment of the bottle full of a certain quantity of water, was not enough to reach a 100% MVC, so that we can affirm that a not complete occlusion of the blood flow was performed also during the slow phase. This fact is reflected by the different behaviors among the subjects of tHb, at least during the first 60-100 seconds of the exercise. After this interval we can recognize a general increase of the tHb, except for subject 4. The increase in ∆tHb with muscle fatigue and local muscle fatigue level was already reported in literature [73]. Furthermore, changes in muscle blood flow due to muscle fiber architecture during muscle contraction may affect significant relationships between the decrease in the slope of MF and maximal changes in O 2 Hb and HHb [23]. In order to overcome this problem, it is necessary in a future study, to integrate the acquisitions with a blood flow sensor, in order to understand what is really happening in each subject. For example the diffuse correlation spectroscopy (DCS) technique could be a good candidate, since it is a non-invasive optical technique which allows real time measurements of the blood flow [74]. This technique could be also used in combination with Doppler acquisitions over major arteries in order to estimate the input function for oxygenated blood flow rates [75].
In order to discriminate between fast and slow phase, different strategies are suggested in literature. For example, it is possible to consider the inflection point, i.e. the point when the initial fast change rate of HHb and O 2 Hb decelerates considerably. It is possible to consider also the inflection duration, i.e. the duration from the start of the exercise until the inflection point [29]. Felici et al. identified automatically the beginning and end of the fast phase with a 7-degree polynomial interpolation [65]. In this paper, we used TD-NIRS parameters for signal segmentation and in particular to identify the end of the transitory phase. This approach is reasonable since it is NIRS-driven, coherently with the objective of the work. However, for future developments, it could be equally reasonable to use MF as a signal to detect fatigue and distinguish phases. However, interestingly enough, in our trials MF showed a monotonic, linear decreasing trend, apparently showing a continuous augmentation of fatigue -except a slight decrease in MF slope for a few subjects when approaching the end of the test (close to complete exhaustion). This result was not favorable to the use of MF as a segmentation signal. In fact, MF did not allow us to clearly and univocally distinguish the transient phase found in TD-NIRS. It is thus possible that, despite the slower dynamics, some TD-NIRS signals may be appropriate to distinguish fatigue phases. Further experiments should be designed to investigate this issue in more detail.
We think that all these considerations can be very useful, in particular to drive the choice of the best NIRS parameters to employ for the assessing of the sustained fatigue during isometric exercises. This can help in more fields of research: from the sport sciences for evaluating the athletes training or in clinical setting, where it is necessary, for example, to assess the rehabilitation progresses of different kinds of populations (athletes after injuries or elder people after fractures) [76]. Considering also the fact that, thanks to the recent technological developments, cheaper, more compact and wearable TD-NIRS devices are now available [77], this analysis can be of interest in newer applications, where a portable system is required.

Future works and limitations
While providing novel evidences on sustained fatigue with the combined use of EMG and NIRS, this study leaves some open points that can be investigated and overcame in future works. First, it will be interesting to better assess the physiological outcome about muscular metabolism during the fatigue process from the TD-NIRS signal, comparing it with other gold standard techniques used for the assessment of fatigue. Some possible correlations can be investigated on the exhaled breath with a metabolic cart to assess the anaerobic and the ventilator threshold [22], on the blood metabolite with the analysis of a blood sample [78], or with the photoplethysmography and electrocardiographic technique [79]. Moreover, another open point is to determine a unique parameter for the detection of the onset of fatigue to perform uniform analysis on the dataset. Our future research will investigate the meaningfulness of estimating a single fatigue onset by combining the multi-modal metrics proposed in this study into a common denominator, possibly also including other fatigue-related signals. These synthesis process in the analysis might prove useful when contextualized in pragmatic clinical scenarios and for uniform and straightforward interpretation of the data.
Lastly, we acknowledge that this work has some limitations. The number of subjects enrolled for this study is quite low and can be increased for more reliable conclusions. Moreover, increasing the number of subjects would allow to assess in more detail whether relevant differences can be found in sustained fatigue assessment when considering subject-specific parameters (e.g. gender, age, BMI) that may affect the analysis.

Acknowledgment
This work was supported by Regione Lombardia and Fondazione Cariplo in the framework of the project EMPATIA@Lecco -EMpowerment del PAzienTe In cAsa -Rif. 2016-1428.
Disclosures D.C. and A.T.: pioNIRS S.r.l., Italy (I). Other authors declare no conflicts of interest related to this article.