Participants

Sixty participants with low to moderate levels of physical activity were recruited for the study, including 30 young adults (15 females; mean ± SD, age: 24 ± 4 years, stature: 1.73 ± 0.10 m, mass: 67 ± 11 kg, body mass index: 22.4 ± 1.9 kg/m2) and 30 older adults (15 females, age: 73 ± 4 years, stature: 1.67 ± 0.10 m, mass: 73 ± 12 kg, body mass index: 26.4 ± 2.8 kg/m2). Older adults were shorter than young adults (p = 0.0216) but had greater mass (p = 0.0364) and body mass index (p < 0.0001). All participants were free from cardiovascular, neurological, or respiratory disorders, or musculoskeletal injury of the lower limbs. Participants were eligible for the young group if they were aged 18–35 years, and for the older group if they were aged over 65 years. Additionally, participants were ineligible if they were sedentary or engaged in systematic training or competition, including resistance training, endurance training and team sports. Written informed consent was obtained from all eligible participants. The study was approved by the Loughborough University Ethics Committee (2021-5361-4724), and all study procedures were performed in accordance with the Declaration of Helsinki except for registration in database.

All participants completed the International Physical Activity Questionnaire (IPAQ) [32] to estimate daily physical activity level, which was not different between groups (3059 ± 1696 vs. 3926 ± 2238 MET·min/week; p = 0.096). Older participants’ frailty status was determined by adapting the frailty phenotype guidelines of Fried et al. [33], using the following five criteria: (1) sarcopenia, determined by quadriceps peak cross-sectional area relative to young individuals (z-score ≤ −2) [19]; (2) exhaustion, determined via a questionnaire (answered yes ≥ 3–4 days per week to either: ‘I felt that everything I did was an effort’ or ‘I could not get going’); (3) slowness, determined by a 4 m walk test (≤ 0.8 m/s) [34]; (4) weakness, determined by grip strength test (in the lowest 20% by sex and BMI; Jamar Hydraulic Hand Dynamometer, set to position two); and (5) physical activity, determined via IPAQ (estimated expenditure for males < 383 kcals per week, or estimated expenditure for females < 270 kcals per week). Participants were considered robust if no criteria were met, pre-frail if 1–2 criteria were met, and frail if ≥ 3 criteria were met. Out of 30 participants, we found that none of the older adults were classified as frail, 18 were pre-frail and 10 were robust (n.b. in two individuals, calculation was not possible due to no imaging data; for details see below).

Study protocol

Participants visited the laboratory on four occasions: a familiarisation session, a dorsiflexion testing session, a knee extension testing session, and a magnetic resonance imaging (MRI) scan. In the familiarisation session, participants completed the same tasks as in the testing sessions to become accustomed to the maximal contractions and force-matching tasks during unilateral isometric dorsiflexion and knee extension with their dominant leg. For each session participants were instructed to arrive well rested, having refrained from strenuous activity the day prior, and to avoid caffeine on the day of the testing sessions.

MRI

T1-weighted axial MR images were taken from the twelfth thoracic vertebrae to the calcaneus of both legs with a 3 T scanner (Discovery MR750w; GE Healthcare, Chicago, IL, USA). Images were acquired in five overlapping blocks (time of repetition 600 ms, time of echo 8 ms, field of view 450 × 450 mm, image matrix 320 × 320, pixel size 1.4 × 1.4 mm, slice thickness 5 mm, interslice gap 5 mm). Participants lay supine, with the ankle at 90˚, the knee and hip fully extended, and the arms folded across the chest in a whole-body coil. Oil-filled capsules were placed laterally on the right leg at equidistant intervals to facilitate alignment of the overlapping blocks during analysis.

Two older and one young adult were unable to undergo MRI due to contraindications. Additionally, one young individual’s scan was excluded from quadriceps analysis due to an artefact in the images. Therefore, the results of MRI analysis are based on n = 57 for anterior compartment of the shank and n = 56 for quadriceps femoris. Additionally, for calculation of frailty phenotype (see ‘Participants’ section) the sample size was n = 28.

Force measurements

All force measurements were made with two rigid custom-built isometric dynamometers, with straps tightly fastened over the pelvis and shoulders. Dorsiflexion force recordings were taken with the participants seated with 120° of hip flexion and the knee fully extended (180°) in a custom-built chair fitted with an ankle dynamometer (NEG1, OT Bioelettronica, Torino, Italy) with the ankle at 10° plantarflexion (0° = neutral). The knee joint was strapped down to minimise any knee flexion during dorsiflexion. The foot was secured to a rotating wooden footplate, with proximal (over the dorsum of the foot) and distal (on the distal third of the metatarsals) foot straps and with the heel supported by a small cup. The ankle joint centre was aligned with the axis of rotation of the footplate. The footplate was fixed in place by a strain gauge (CCT Transducer s.a.s., Torino, Italy) arranged in series. The analogue force signal was amplified (× 200; Forza-B, OT Bioelettronica, Torino, Italy) and sampled at 2048 Hz using an analogue-to-digital (A/D) converter (Micro 1401–4; CED Ltd., Cambridge, UK).

Knee extension recordings were performed with participants seated with the knee and hip flexed at 115° and 125°, respectively (180° = full extension). Force recordings during knee extension were measured with a calibrated S-beam strain gauge (linear range 0–1.5 kN; Force Logic, Swallowfield, UK). During knee extension the strain gauge was attached perpendicular and posterior to the tibia with a strap (35 mm reinforced canvas webbing) placed ~ 3–4 cm above the medial malleolus. The analogue force signal was amplified (× 370) and sampled at 2048 Hz using an A/D converter (Micro 1401-4; CED Ltd., Cambridge, UK).

The force signal during both knee extension and dorsiflexion was acquired in Spike2 software (version 10; CED Ltd., Cambridge, UK). Force was also concurrently sampled (2048 Hz) and recorded in OTBiolab + software (Quattrocento; OT Bioelettronica, Torino, Italy) using a multichannel amplifier (16-bit, Quattrocento; OT Bioelettronica, Torino, Italy) for synchronisation with high-density surface electromyography (HDsEMG) recordings. Visual force biofeedback was provided during all experimental procedures on a monitor ~ 1 m in front of participants.

Electromyography recordings

Each testing session began with identification of the motor point on the VL for knee extension and TA for dorsiflexion. The motor point was identified as the site on the muscle that produced the greatest twitch in response to low intensity percutaneous electrical stimulation (maximum output 400 V, 0.2 ms pulse width; DS7R; Digitimer, Welwyn Garden City, UK).

Preparation of the skin over the VL and TA involved shaving and abrasion. After that, 64-channel HDsEMG electrode grids (13 × 5, 1 mm electrode diameter, 8 mm interelectrode distance; OT Bioelettronica) were covered with a bi-adhesive foam layer (Spes Medica, Battipaglia, Italy) with the cavities filled with conductive paste (Ten20, Weaver and Co., Aurora, CO, USA). The HDsEMG electrodes were placed on the muscle bellies aligned with the presumed muscle fibre orientation. The VL electrode grid was placed distal of the motor point at 30° from a line drawn between the anterior superior iliac spine and the lateral patella [35]. The TA electrode grid was placed distal to the motor point over the belly of the TA muscle. Reference electrodes were placed over the fibular head for knee extension and over the medial malleolus for dorsiflexion, with a dampened strap around the non-dominant ankle used to ground the signal. HDsEMG signals were recorded in a monopolar configuration, bandpass filtered (10–500 Hz), sampled at 2048 Hz using a 16-bit multichannel amplifier (Quattrocento; OT Bioelettronica, Torino, Italy), and acquired using OT Biolab + software (OT Bioelettronica, Torino, Italy).

Experimental procedures

Identical dorsiflexion and knee extension tasks were completed. Participants began with a brief warm up consisting of 3 × 50%, 3 × 75% and 1 × 90% of perceived maximum effort contractions, each lasting ~ 3 s and separated by 15–30 s rest. After that, four contractions with maximal effort were performed separated by ≥ 30 s rest, in which participants were asked to perform the contractions “as hard as possible” for 3–5 s. During the maximal effort contractions, participants received strong verbal encouragement and visual force biofeedback was provided. Participants were instructed to try and beat their previous highest score (indicated by a horizontal cursor on the monitor) during successive efforts. The MVF was taken as the highest instantaneous force achieved during any contractions with maximal effort.

Participants then completed several submaximal contractions in which they were asked to follow a force target as closely as possible. Specifically, participants were tasked to perform triangular-shaped ramp contractions by following a target trace with a 10 s linear ascending and then 10 s linear descending slope, peaking at 15, 30, 50, and 70% MVF. Two contractions were performed at each force level and separated by 30–60-s rest.

Data analysis

Offline analyses were performed using MATLAB (R2021b; Mathworks Inc., Natick, MA, USA).

Force signal

The force signal was filtered using a fourth-order, zero-lag low-pass Butterworth filter with a cut-off frequency of 20 Hz and gravity corrected.

HDsEMG analysis

High-density surface EMG signals were inspected, and poor-quality channels were removed based on their amplitude and area under the power spectrum using a semi-automated tool in MATLAB. The signals were then decomposed using the Convolution Kernel Compensation algorithm to obtain the individual MU spike trains [36]. This algorithm uses blind source separation principles, inverting the EMG mixing model to estimate the MU filter, which represents weights of the linear spatial-temporal combination of HDsEMG channels and generates an estimation of each individual MU spike train [37, 38]. Decomposed signals at the same force level were then concatenated and MU filters transferred to enhance MU yield [39, 40] followed by manual inspection and editing by trained operators as previously described [41]. Any duplicate MU spike trains that were detected as sharing 30% of the same MU discharges (discharge tolerance of 0.5 ms, i.e. 1 sample) were deleted based on lower estimated decomposition accuracy assessed by pulse-to-noise ratio. Motor unit spike trains were kept for further analysis if they had a pulse-to-noise ratio of 30 dB or greater and a consistent discharge pattern throughout the contraction [42].

Motor unit tracking

The discharge patterns of MUs were analysed independently across contraction force levels (see supplementary material) and with MU tracking to account for the potential replicates across contraction levels. A similar procedure to that described by Škarabot et al. [31] was used to track MU spike trains across contraction levels. Briefly, separation vectors (MU filters) were used to detect spike trains of MUs that were identified across multiple contraction levels. Specifically, the decomposition results from each contraction level were concatenated, and the MU filters were applied from signals recorded during lower to those recorded during higher contraction levels to identify common MU spike trains. Contractions at 15% MVF were excluded from this part of analysis, as it was unlikely any MUs could be tracked from 15 to 70% MVF due to the inherent sampling bias of HDsEMG decomposition (i.e. difficulty in identifying low-threshold MUs in the presence of higher-threshold MUs with large potentials) [37]. As part of the tracking process, duplicate MU spike trains were deleted using the same criteria as above (see HDsEMG analysis). Tracking low threshold MUs whilst high threshold MUs are concurrently active can be challenging [37], and therefore, we tracked MUs from neighbouring contraction levels 30–50% MVF and 50–70% MVF as well as 30–50–70% MVF to enhance MU yield (Fig. 1). The results were found to be similar with or without tracking of MUs across contraction intensities; therefore, we chose to principally present results of tracked MUs, whereas the results of non-tracked MUs are located in the Supplementary Material.

Fig. 1

Examples of motor unit discharge behaviour across contraction levels in young and older adults. For each group, i.e. young (upper panels, blue) and older adults (lower panels, orange) a raster plot of tracked motor unit (MU) discharges is shown in the centre during triangular contractions to 30, 50, and 70% of maximal voluntary force (MVF). In the upper tier of each group’s panel block, the ascending discharge rate of representative units expressed as a function of force is shown, along with the calculated brace height, and the acceleration and attenuation slopes. In the lower tier of each group’s panel block a representative pair of smoothed MU discharge rates is shown, with onset-offset hysteresis (∆F) annotated for each test unit, and with κ representing the maximal discharge rate modulation of the control unit between test unit recruitment and control unit derecruitment, to which ∆F were normalised

Motor unit discharge properties

For all MU discharge-related outcome variables, MU spike trains from triangular ramp contractions were smoothed with support vector regression to produce a continuous estimate of MU discharge patterns using parameters as previously described [43]. Peak discharge rate was calculated as the peak of the smoothed discharge rate. The relative force produced at the time a given MU was recruited (i.e. the first spike in the spike train was identified) was considered the recruitment threshold.

Analysis of motor unit discharge patterns

A hallmark of PIC-mediated actions on motoneurons is the prolonged effect of synaptic input. In cases of linearly increasing and decreasing excitatory synaptic input, the prolongation of MU discharge is evident as the discharge onset-offset hysteresis, whereby a motoneuron is derecruited at lower levels of synaptic input than was required for its recruitment [44]. In human recordings, onset-offset hysteresis can be quantified as the difference in discharge rate of a lower-threshold (control) unit at the instants of higher-threshold (test) unit recruitment and derecruitment (ΔF) [26, 45], which we calculated in the present study. In the case of non-tracked MUs, we used a ‘unit-wise’ approach which provides a single, mean ΔF value for a given test unit which may be paired with multiple control units [30, 46]. However, in the case of tracked MUs, the calculations of ΔF were made for individual pairs (Fig. 1) in order not to bias the estimation of hysteresis (e.g. to avoid a potential case of different number of tracked test units being deemed suitable at different contraction levels due to subtly meeting the inclusion criteria). A pair of MUs was deemed suitable for calculation of ΔF if (1) the test unit was recruited > 1 s after the control unit, to ensure full activation of control unit PICs [47], (2) the control unit discharge rate was > 0.5 pps whilst the test unit was active, to reduce the influence of MU saturation on ΔF calculation [48], and (3) the rate-rate correlation of pairs of MUs was r2 > 0.7 [49], to increase the likelihood that the MU pair shared common input. To account for potential age-related differences in descending MU discharge rate modulation of control units and/or across contraction levels, we additionally normalised ΔF values to the maximal discharge rate modulation for each test-control unit pair (i.e. ΔF/κ, where κ represents the difference in control unit discharge rate between test unit recruitment and control unit derecruitment, Fig. 1) [31].

To further appreciate the contribution of PICs to MU discharge rate patterns, we examined the ascending MU discharge modulation via quasi-geometric analysis. Specifically, we quantified the non-linearity of the ascending MU discharge rate (termed ‘brace height’) to disentangle the relative contribution of neuromodulation and excitation-inhibition patterns to MU discharge [50]. That is, utilising in silico studies with biologically realistic motoneuron pools, it has been shown that the ascending discharge rate non-linearity is associated with neuromodulation, whereas onset–offset hysteresis is dependent on both neuromodulation and the pattern of inhibitory synaptic input relative to excitation [31, 50, 51]. For this analysis, the smoothed discharge rate of any given MU was expressed as a function of force; from there, the relative deviation from linearity (brace height) was quantified as the maximum orthogonal vector between the smoothed discharge rate and the theoretical linear discharge (a linear line between onset and peak firing; Fig. 1). To minimise the confounding effects of the absolute ascending discharge rate modulation (e.g. when comparing individuals of different ages that are likely to exhibit different peak discharge rate, or when examining different contraction levels), the brace height values were normalised to a right triangle with a hypotenuse between onset and peak discharge rate.

The quantification of brace height provides additional metrics which may elucidate a deeper insight of the control of MU discharge behaviour. Specifically, we measured the slope of the initial acceleration phase of MU discharge and the following attenuation phase (Fig. 1), which were segmented at the point of brace height occurrence [50]. The acceleration phase corresponds to the secondary range of MU discharge and features an increase in MU discharge rate, thought to be due to PIC amplification [52]. The subsequent attenuation phase corresponds to the tertiary range of MU discharge, which is likely when PICs are fully activated, meaning that further increases in synaptic input result in an attenuated increase in discharge rate [53]. In addition to brace height and ΔF metrics, these variables were quantified to potentially allow us to decouple the influence of neuromodulation and pattern of inhibition to motoneurons. Specifically, in simulations with a fixed excitation level, the attenuation slope has been shown to be sensitive to the pattern of inhibition whereas the acceleration slope may be more sensitive to the amount of neuromodulation and the total inhibitory input [31, 50].

Motor units were excluded from the quasi-geometric analysis if they had a negative slope during the acceleration phase, a normalised brace height greater than 200%, or if peak MU discharge rate was reached after peak force.

MRI

The MR images of the dominant leg were analysed for the anterior compartment of the shank (tibialis anterior, extensor digitorum longus and extensor hallucis longus) which was segmented as one whole and quadriceps (rectus femoris, vastus lateralis, vastus medialis, and vastus intermedius) in which each constituent quadricep was segmented individually and then summed. This analysis was performed using an open access DICOM image-analysis software (HOROS, version 2.2.0, www.thehorosproject.org). From the axial images, the image with the subjectively largest anatomical cross-sectional area (ACSA) and two slices either side (proximal and distal) were first manually segmented for each muscle/compartment, with additional slices segmented according to the ACSA values to ensure that at least two slices were segmented either side of the largest ACSA. A cubic spline function was fit to these five ACSAs generating 1000 interpolated values (R2021b; Mathworks Inc., Natick, MA, USA). The peak ACSA was defined as the maximum interpolated value. In the quadriceps, the sum of each constituent muscle peak ACSA was defined as maximum quadriceps ACSA.

Statistical analysis

All statistical analysis was performed in R (RStudio, v 1.4.1106, R Foundation for Statistical Computing, Vienna, Austria).

To assess age-related differences in MU discharge characteristics, linear mixed effects models (lme4 package; ver. 1.1.35.3, [54]) were constructed with age group (young vs. older), muscle (TA, VL), contraction level (30, 50, and 70% MVF) and their interaction as fixed factors with participant identifier and MU identifier as random intercepts, which were nested in a hierarchical structure. We also included MU recruitment threshold as a covariate to control for its influence on MU discharge behaviour [55]. Additionally, although sex differences were not the focus of this study, to account for potential differences between males and females [46], sex was used as a covariate in the statistical model. For example, the Wilkinson notation for the model with peak discharge rate as the outcome variable was as follows: Peak discharge rate ~ Age Group × Muscle × Contraction Level + Recruitment Threshold + Sex + [1 | Participant ID] + [1 | Participant ID | MU ID]. For ΔF, in order to account for the influence of both control and test units the model was as follows: ΔF ~ Age Group × Muscle × Contraction Level + Test Recruitment Threshold + Control Recruitment Threshold + Sex + [1 | Participant ID] + [1 | Participant ID | Test MU ID] + [1 | Participant ID | Test MU ID]. Quantile–quantile plots and histograms were used to verify the assumptions of normality and linearity of residuals. Observations with studentised residuals exceeding ± 3 were identified as outliers and excluded from the models. Square root (sqrt) and log transformation of dependent variables were used when residuals were found to be non-normally distributed. For variables containing negative values, a constant was added, which was the smallest whole number required to shift all values into the positive range before transformation. The selection of transformation function was determined based on visual inspection of residuals (as described above) and on model fit using Akaike Information Criterion and Bayesian Information Criterion. The following variables were transformed: peak discharge rate (log), ΔF (sqrt + 8), brace height (log), acceleration slope (sqrt), and attenuation slope (log + 1).

Similarly, linear mixed effects models were used to assess if there were differences in MVF due to factors of age group and joint action (knee extension and dorsiflexion), and to assess differences in muscle peak ACSA due to age group and muscle group (quadriceps and anterior compartment), with participant identifier as a random intercept and sex as covariate in both instances. Both MVF and peak ACSA underwent log transformation. Linear models were used to assess differences in IPAQ score, MU counts and participant characteristics (height, body mass and body mass index).

The significance of the model fixed effects and interactions was assessed using Type II ANOVA based on Wald chi-square tests from the car package (ver. 3.1.3). Significant main effects and interactions were further investigated using simple effects, interaction contrasts and pairwise comparisons of estimated marginal means (emmeans package; ver. 1.11.1, [56]) with Bonferroni-adjusted p-values to control for multiple comparisons. Cohen’s d, calculated as the difference in estimated marginal means divided by the residual standard deviation of the linear mixed effects model, was used to estimate effect sizes of post hoc comparisons. Significance was set at an alpha level of 0.05. Data are presented as estimated marginal means [95% confidence interval].