Model-based analysis of fatigued human knee extensors

This study investigated the effect of isometrically induced fatigue on Hill-type muscle model parameters and related task-dependent effects. Parameter identification methods were used to extract fatigue-related parameter trends from isometric and ballistic dynamic maximum voluntary knee extensions. Nine subjects, who completed ten fatiguing sets, each consisting of nine 3 s isometric maximum voluntary contractions with 3 s rest plus two ballistic contractions with different loads, were analyzed. Only at the isometric task, the identified optimized model parameter values of muscle activation rate and maximum force generating capacity of the contractile element decreased from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20.8 \pm 8.4$$\end{document}20.8±8.4 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$11.2 \pm 4.1$$\end{document}11.2±4.1 Hz and from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18{,}137 \pm 150$$\end{document}18,137±150 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10{,}666 \pm 2139$$\end{document}10,666±2139 N, respectively. For all tasks, the maximum efficiency of the contractile element, mathematically related to the curvature of the force–velocity relation, increased from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.35 \pm 0.04$$\end{document}0.35±0.04 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.42 \pm 0.05$$\end{document}0.42±0.05. The model parameter maximum contraction velocity decreased from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.93 \pm 0.1$$\end{document}0.93±0.1 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.9 \pm 0.1$$\end{document}0.9±0.1 m/s and the stiffness of the serial elastic element from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1936 \pm 227$$\end{document}1936±227 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1432 \pm 245$$\end{document}1432±245 N/mm. Thus, models of fatigue should consider fatigue dependencies in active as well as in passive elements, and muscle activation dynamics should account for the task dependency of fatigue.


Introduction
Human movement is possible through the interaction of several structures including nervous system, muscles, tendons, and bones. As a consequence of high-intensity and longduration contractions, interdependent systems of the body will react to the physical strain caused.
Since the early work of Mosso (1906), it is known that fatigue affects the central and peripheral systems simultaneously (Bigland-Ritchie et al. 1978;Enoka and Duchateau 2016;Gandevia 2001;Taylor et al. 2016). In voluntary movements, the fatigue-related changes in both systems cause the reversible reduction of measures Communicated by Toshio Moritani.
* Harald Penasso harald.penasso@uni-graz.at 1 Institute of Sport Science, University of Graz, Mozartgasse 14, 8010 Graz, Austria of performance, e.g., power and force (Edman and Mattiazzi 1981;Hilty 2011;Williams and Ratel 2009). This is accompanied by an increased effort to maintain a sub-maximal task as well as by the sensation of exhaustion (Enoka and Stuart 1992;Noakes 2012). While the central system controls the level of voluntary activation, changes at or distal to the neuromuscular junction are ascribed to the peripheral system (Gandevia 2001). Furthermore, the magnitude and the effect of fatigue will depend on the modality of the movement, e.g., task dependency (Enoka and Stuart 1992;Taylor et al. 2016), contraction type (Babault et al. 2006;Griffin et al. 2000;Liu et al. 2003), contraction velocity (Harwood et al. 2011;Morel et al. 2015), and the level of achieved voluntary activation (Bigland-Ritchie et al. 1978;Gandevia et al. 1998;Sidhu et al. 2013).
The effects of fatigue on the interacting sub-structures of Hill-type muscle models, however, are not yet well studied. Models applied to fatiguing exercises should consider central as well as peripheral systems (Callahan et al. 2016;Hawkins and Hull 1993;Shorten et al. 2007). Central effects can be described by modeling the dynamics of motoneuron recruitment (Liu et al. 2003), and the peripheral force generation of the muscle-tendon complex ( MTC ) should consider changes at the contraction velocity-dependent contractile element (CE) (Jones et al. 2006) and the coupled serial elastic element ( SEE ) (Kubo et al. 2001).
The purpose of this study was to investigate the simultaneous effects of isometrically induced fatigue on Hill-type muscle model parameter values and muscle activation. Parameter identification tools were used to optimize model parameters within physiological limits to fit the measured force and position data. As some of the parameter values of activation dynamics, force-velocity relation (FV) and SEE were separately identified for isometric (ISO) and dynamic ballistic (BAL) maximum voluntary contractions (MVCs); also, task-dependent effects could be addressed.

Methods
The subjects performed isometric and dynamic ballistic knee-extension MVCs at an instrumented inclined leg press and completed an intermitted isometric fatigue protocol. At each isometric MVC, the subjects were encouraged with a pre-recorded loud '3-2-1-go!' sound (McNair et al. 1996). For all stages of fatigue, data from three contractions were simultaneously processed by a nonlinear parameter identification routine that fitted the kneeextension model output to the measured data (Penasso and Thaller 2017;Siebert et al. 2007).

Experimental setup
The subjects were seated and fastened to an inclined leg press (leg-press inclination 10 • ; Tetra ® Ilmenau). The sledge of the leg press (32.5 kg + 2 × 20 kg additional load) could be allowed to slide or to be fixed at a desired position for isometric contractions and carried a force plate (Kistler type Z17068AA1; range F x,y 0-2.5 kN, F z 0-10 kN; sampling rate (SR) 500 Hz; Fig. 1). During ISO and BAL contractions, the sledge was positioned such that the initial knee-extension angle ( ) corresponded to 120 • (cf. Fig. 1). To prevent additional force transmission from unwanted plantar flexion during the knee extensions, a 5 cm block of wood, used as a spacer, was mounted onto the force plate and served as heel contact. The sledge position was measured with a high-precision rotary encoder (Stegmann 4096 ppr, SR 500 Hz). Motion capture data obtained from reflective markers, placed at the greater trochanter, the lateral epicondyle of the knee, and the lateral malleolus (6 Vicon M2 cameras; SR 250 Hz), were used to correct subject displacements relative to the 1D leg-press coordinate system. Following Sust et al. (1997b), anthropometric measures of the thigh (greater trochanter to lateral knee joint line), the shank (lateral knee joint line to lateral malleolus), the knee circumference (CIR, measured at the height of the center of the patella to estimate the knee moment arm r = CIR∕(2 ) ), and the approximated length of the patellar tendon ( l p , patella center to tuberositas tibiae, measured at the fully extended and relaxed leg) were taken with a tape measure (Fig. 1). Superficial muscle activities (sEMG) from the right m. vastus lateralis (VL), m. vastus medialis (VM), and m. rectus femoris (RF) were Fig. 1 Sketch of the muscle-tendon complex ( MTC ) and of the geometric representation of the leg including the knee moment arm (r) as well as the lengths of thigh ( l thigh ), shank ( l shank ), and patellar tendon ( l p ) (proportions not in scale). Using the position (X) of the sledge (m), the knee-extension angle ( ) is calculated and is calculated numerically. In the model, a movement is initiated by activating the contractile element (CE) via activation dynamics (A(t)). Thus, the force of the MTC , produced by the interaction of CE and serial elastic element ( SEE ), is transferred to the environment ( F mod ) acquired according to SENIAM guidelines for bipolar configurations (myon 320; SR 4000 Hz) (Hermens et al. 2000). The activity of the deep vastus intermedius muscle was not measured. All data were synchronized and recorded on one measurement system (Dewetron).

Procedure
Subjects were familiarized with the testing procedure one week before the experiment.
On the measurement day, the anthropometric measurements were taken, the sEMG electrodes as well as the reflective makers were placed and the subjects warmed up (5 min self-paced, self-rated low-intensity ergometer, 1 min severe ergometer exercise, 30 s crouching at 90 • knee-angle, and ten squats). Subsequently, they sat down at the instrumented inclined leg press and were fastened to the seat. The measurement started with two ballistic MVCs against low loads (B20, 32.5 + 20 kg) and continued with two BALs against higher loads (B40, 32.5 + 40 kg; Fig. 2). Between these short contractions, the subjects recovered for about 10 s.
The initial procedure was followed by the fatiguing protocol consisting of ten identical sets (Fig. 2). Thereby, the number of fatiguing sets was not revealed to the subjects who agreed to go all out at each contraction. Each set included nine 3 s ISO which were separated by 3 s rest periods. The ISO tasks were automatically initiated and terminated by a pacing-trigger signal that played the 3 s long '3-2-1-go!' sound. The corresponding trigger signal was also shown as green (contract) or red (relax) background color at a visual feedback display where the current force was shown as bar graph. Each 9th ISO was immediately followed by two BALs [(1) B40, (2) B20]. The resting period prior to the BALs was about 5 s. Together with the manual unlocking and fixation of the sledge, the duration of one fatigue set was about 1 min.

Subjects
Fifteen healthy male students of sport science completed the testing protocol correctly, which was checked using a usual documentation video (200 fps). The parameter identification results of six subjects had to be discarded because the final parameters values had converged toward solver boundary constraints. Consequently, the datasets of nine subjects remained ( 27 ± 3 years, 183 ± 8 cm, BMI 22.4 ± 1.6 kg/m 2 ) and were analyzed throughout the fatigue task.

Data processing
The data processing was performed in Matlab ® 2016a. Force and position data were low-pass filtered at 15 Hz and 30 Hz for ISO and BAL contractions, respectively (zero-lag, butterworth, 4th ) (Maffiuletti et al. 2016). The onset of force was automatically detected by a threshold which was set to the mean of the manually annotated silent period between the contractions plus 80 N (which is less than 5% of the max. force).
Some sEMG data showed frequency artifacts which were removed at 200 and 400 Hz using an adaptive notch filter. To reduce the overall standard deviation of the sEMG amplitude value, it was defined as the mean of two adjacent 0.05 s time windows which enclosed the maximum force per contraction (Merletti et al. 1990). At each period, the mean sEMG amplitude was calculated using wavelet analysis which represents a band pass filter between 6.9 and 395.5 Hz (Von Tscharner 2000) and was defined as the mean of the sum of the rectified wavelets (Borg 2010). For ISO, the sEMG muscle onsets and terminations were detected automatically by using the Teager-Kaiser energy operator and a threshold method (Kaiser 1990;Solnik et al. 2010). For BAL, the sEMG onsets were determined in the same way as for ISO, but due to the ongoing muscle contraction after push-off, the BAL terminations had to be set manually when force dropped to zero. Due to re-positioning and other little activities of the subjects during the breaks, the baseline signal had to be annotated by hand. Onsets which could be detected automatically were set by hand via a graphical interface.

Mathematical model of the knee-extension task
The mathematical model and the parameter identification routine are discussed in detail in Penasso and Thaller (2017). This section provides a brief overview.
The model is formulated as an equation of motion that describes the time evolution of concentric and isometric kneeextension tasks as an ordinary differential equation which corresponds to the experimental setup ( Fig. 1, Eq. 1). Each The measurement protocol started with two ballistic contractions against low load (B20, 32.5 + 20 kg) and continued with two ballistic MVCs against higher load (B40, 32.5 + 40 kg) separated by ∼ 10 s breaks. This was followed by the fatiguing protocol consisting of n = 10 identical sets. Each set included nine 3 s ISO which were separated by 3 s rest periods. Each 9th ISO was followed by two BALs [(1) B40, (2) B20]. The resting period prior to the BALs was about 5 s. Together with the manual unlocking and fixation of the sledge, the duration of one fatigue set was about 1 min. The preceding familiarization, preparation, and warm-up procedures are not shown movement is initiated at the onset ( t 0 ) of the time-dependent muscle activation dynamics (A(t)) which consider the progressing recruitment of fibers after onset (Eq. 2). A(t) considers the pre-activation level ( A 0 ), the number of maximally available fibers ( n max ), and the rate constant (U) (Eq. 2). Once activated, the dynamics of the muscle-tendon complex ( MTC ) produce force ( f MTC ), where the actively force-generating element is the contractile element (CE). The CE is modeled as a contraction velocity ( ̇l CE )-dependent force-velocity relation (FV) giving the CE force ( f CE (l CE ) ) (Hill 1938) (Eq. 3). To consider length dependencies, its force is scaled by a parabolic force-length relation ( FL(l CE ) ), which is a function of the CE length ( l CE ) (Van Soest and Bobbert 1993) (Appendix A). It is assumed that the optimal length of the CE, representing a muscle fascicle, is 0.09 m which is reached at a knee-extension angle of 120 • (Hoy et al. 1990). To be able to determine the value of the linear stiffness of the SEE, it had to be set linear. The SEE depends on its length ( l SEE ) (Eq. 4) where the slack length is defined as the distance from the center of the patella to the tuberositas tibiae (cf. Sect. 2.1.1). The ratio of the external force of the model ( F mod ) and f MTC is given by G(X), which was published earlier (Sust et al. 1997b;Siebert et al. 2007;Wagner et al. 2006;Thaller et al. 2010) (Appendix B). Thus, F mod (t) changes the position (X) of the point mass (m) during dynamic contractions. The displacement of m represents the movement of the sledge on a 10 • -inclined leg press, where g = sin(10[ • ]) × 9.81 [m∕s 2 ] . Thereby, A(t) is the solution of two coupled first-order differential equations, which consider time (t)-dependent excitation and recruitment of inactive muscle fibers.
which represents a normalized smooth sigmoid-like function ranging between the pre-activation level ( A 0 ) and full activation ( n max = 1 ) with the slope controlled by the parameter of the rate constant of muscle activation (U). The FV is modeled according to Siebert et al. (2007) using Hill's parameters a, b, and c where G(X) is used to convert the external velocity ( Ẋ ) to the velocity of the MTC ( v MTC = G(X) ⋅Ẋ ), which is then used to calculate ̇l CE by considering the velocity of the SEE where the parameter k SEE,lin represents the linear stiffness.

Parameter identification
To calculate the model parameters for all stages of fatigue, the parameter identification routine minimized the sum of squared errors between F z of the force plate ( F meas ) and the external force of the model ( F mod (t) ). It processes datasets containing isometric and dynamic ballistic contractions from different loads (ISO, B20, B40) at the same time, equally weighted. Since ISO and BAL contractions are differently affected by fatigue (Griffin et al. 2000;Morel et al. 2015), a contraction-type-dependent version of the identification routine was first tested using simulated data and then used to identify parameter sets of the fresh and fatigued states. This eight degrees of freedom least squares problem (including one t 0 for each of the three contractions per dataset) was solved numerically using the Matlab ® global search class together with the interior-point algorithm of the fmincon solver.
Initially, the parameters of the fresh state are determined (Fig. 3). The dataset consists of three contractions. Out of the first two B20, two B40, and three ISO contractions, each steepest and/or strongest contractions was selected automatically. It was assumed that the number of maximally available fibers corresponds to n max = 100% and thus the parameter maximum isometric force ( f iso ) could be identified (Siebert et al. 2007). In the next step, model parameter values from all datasets are calculated. At each nth fatigued state, the 9th ISO contraction was processed together with the two subsequent dynamic MVCs (B40, B20). To assess if isometrically induced fatigue affects ISO and BAL tasks differently, f iso Fig. 3 Flowchart illustrating the calculation of the direct parameters and the procedure of parameter identification of the fresh state ( f iso,fresh ) was set constant to enable the identification of model parameter values which are specific for BAL and ISO tasks (cf. "Optimized parameters"). Subsequently, each identified parameter set was used to simulate the corresponding movement and the adjusted correlation coefficient ( r 2 a ) was calculated between the resulting simulated F mod and measured data F meas .

Extracted parameters
Fatigue-related parameter trends of the Hill-type model muscle coupled with activation dynamics (optimized parameters) and values calculated from the measured force, position and sEMG data (direct values) were analyzed.

Optimized parameters
Contraction-type-dependent values were identified for the function activation dynamics of the muscle (Table 1). As fatigue causes motor unit activation decrease during maximum efforts (Bigland-Ritchie et al. 1979) while cellular inorganic phosphate reduces contractility (Stutzig et al. 2017), the maximum CE force generation capacity cannot be assigned to either f iso or n max . f CE,max = f iso ⋅ n max , f iso was set to the identified value of f iso,fresh for computational reasons and the parameters n max,ISO , n max,BAL , U ISO , U BAL were identified at subsequent stages of the fatigue protocol ( Fig. 3 and Eq. 2). Besides f iso , two further parameters were required to characterize the FV (Thaller and Wagner 2004). One parameter represents the maximum contraction velocity ( v max ), the other the maximum efficiency ( max ) which is mathematically related to the curvature of the FV. Both were identified independently from contraction type. Furthermore, k SEE,lin was identified to characterize the SEE behavior.
Parameters of FL(l CE ) were not optimized because no substantial fatigue dependency is expected in normalized force-length relations (MacNaughton and MacIntosh 2006). The fatigue-related downward shift of the FL(l CE ) is reflected in the parameters n max,ISO and n max,BAL .

Direct values
Using filtered force ( F meas ) and position data ( X meas ), the external maximum force (F MAX ∶= max {F meas }) , the maximum external power (P MAX ∶= max {F meas ⋅Ẋ meas }) , and the maximum external velocity (V MAX ∶= max {Ẋ meas }) per Spring constant Elasticity of aponeurosis/tendon contraction were obtained ( Table 2). The rate of force development (RFD) was calculated at 50 ms after the detected force onset (Buckthorpe et al. 2014;Maffiuletti et al. 2016). The mean sEMG amplitudes (A) were calculated from the muscles of the right leg (e.g., A VL for m. vastus lateralis). The periods from the initiation and termination signal of the pacing trigger to the detected sEMG onset and termination events defined the reaction times. The duration of each contraction was defined as the duration between sEMG onset and termination.

Statistics
Missing values at the time series were replaced by using the EM algorithm (Dempster et al. 1977). To analyze the time series of each parameter, the optimized parameters and direct values at every contraction number were first tested for the normality of the distribution by using the Kolmogorov-Smirnov test. At least once, the identified values of all subjects were not normally distributed; thus, the nonparametric Friedman test was used to assess if a time evolution changed significantly. To investigate if parameter trends were differently between contraction types, the slopes of the linear trends were calculated. These were normally distributed and thus were analyzed with a one-way ANOVA followed by Tukey HSD pairwise comparisons. The resulting p values are indicated as ***, **, *, t and reflect p ≤ 0.001 , p ≤ 0.01 , p ≤ 0.05 , and p ≤ 0.1 , respectively. The software package G*Power was used to assess the statistical power of each test (Faul et al. 2007), where a value of 0.8 was considered as sufficient statistical power. Effect sizes of the Friedman test were calculated using Kendall's coefficient of concordance ( W K ) and ANOVA effect sizes were calculated using Hedge's g * . Effects where 0.2 ≤ W K or g * < 0.5 are classified as small ( S ), 0.5 ≤ W K or g * < 0.8 as medium ( M ), and W K or g * ≥ 0.8 as large ( L ) and are reported only if the Friedman test or the Tukey HSD was significant.

Results
The identification routine was successfully tested by re-identifying the input parameters of previously simulated datasets (relative errors < 2.3 × 10 −4 %). Thus, the identification of different values of n max for ISO and BAL contractions is in accordance with previous model tests (Penasso and Thaller 2017). In 12.1% of all datasets, the identification routine converged toward bound constraints and thus these data had to be replaced with estimations based on the EM algorithm. The goodness of fit did not change in the course of fatigue (mean r 2 a = 0.99 ± 0.009 , p = 0.23 ; mean normalized rootmean-squared error = 0.91 ± 0.02 , p = 0.15 ) implying that the model was able to fit fresh and fatigued data in equal measure.
To test if the subjects rested sufficiently between the fresh contractions, P MAX of the fresh first BAL contraction was compared to the subsequent second contraction which was 99.7 ± 12.6% of the initial one. During the fatigue protocol, the subjects reacted 41.6% faster to the stop signal (initial 0.88 ± 0.19 s, final 0.37 ± 0.31 s) (0.2% replaced, * * * S , power > 0.99). Thus, the duration of the ISOs changed from 3.6 ± 0.2 s to 3.1 ± 0.3 s (to 86%) (0.2% replaced, * * * M , power > 0.99) and remained greater than the 3 s intended.

Optimized parameters
Using the identified value f iso,fresh = 18,375 ± 3070 N, the optimized parameters of fresh and fatigued states were identified (Table 3). The model parameters U and n max , Table 3 Mean values and standard deviations ( ) of optimized parameters at the fresh and final fatigued states At the isometric task (ISO), the number of active motor fibers ( n max ) and the activation rate constant (U) declined, whereas isometrically induced fatigue did not affect these parameters at ballistic tasks (BAL). Reductions in the parameter values of max. efficiency ( max ), max. contraction velocity ( v max ), and linear stiffness of the serial elastic element ( k SEE,lin ) were found. The parameter description is given in Table 1 Statistical significance, power, and effect sizes are shown in the last row ***p ≤ 0.001 , **p ≤ 0.01 , * p ≤ 0.05 ; † † † power > 0.99 , † † power ≥ 0.9 , † power ≥ 0.8 ; S small effect, M medium effect, L large effect describing the activation dynamics of the model muscle, were identified independently for ISO and BAL contractions. While the number of maximally available fibers n max did not decrease significantly for the BAL tasks, n max,ISO decreased to 58.8 ± 11.9% (Fig. 4a). Considering f iso = f iso,fresh ⋅ n max , the behavior of n max is also valid for f iso . Similar to n max , parameter changes of U were only detected for ISO ( U ISO ), which was reduced to 53.9 ± 22.1% (Fig. 5a).
For all contraction types, fatigue caused an increase of the FV curvature, which is mathematically related to the increase in max to 119.4 ± 19.7% . Since max = p max ∕c, the increase in max is also related to p max . Using abc (Table 1) (Thaller and Wagner 2004), a decrease in p max to 82.2 ± 13.5% is seen (Fig. 4b).
The third parameter characterizing the FV is v max , which was reduced to 96.5 ± 5% (Fig. 6a).

(b) (a)
The linear stiffness of the SEE, identified for ISO and BAL contractions together, and was reduced to 74 ± 7.5% with fatigue (Fig. 5b).

Direct values
Using the maximum value out of the first three ISO contractions and the maximum value of the first two contractions of B40 and B20, the fresh states (100%) were calculated. The majority of the directly measured values calculated from force and position data as well as the sEMG amplitude declined with fatigue (Table 4). The slopes of the linear parameter trends calculated from force and position data indicate that the isometrically induced fatigue affects faster contractions less than slower ones (Table 5). Figure 7a shows that F MAX dropped to 53.2 ± 13 , 80.1 ± 14.8 and 84.8 ± 12.3% at ISO, B40 and B20 contractions, respectively. Thereby, the values of the slopes of the linear trends showed greater reductions in the ISO task than at BAL. The reduction of P MAX at B40 and B20 to 68.6 ± 21.2 and 82.4 ± 15.6% is shown in Fig. 4b. The decline at the slower and higher force requiring B40 task was greater than at the B20 task. Differences in the trends of V max at B40 and B20 contractions are illustrated in Fig. 6b, which declined to 84.3 ± 11.4 and 93.8 ± 7.5% . RFD was reduced to 25 ± 13.4 , 57.4 ± 38.6 and 65.6 ± 27.1% for ISO, B40, and B20, respectively, where ISO showed stronger fatigue effects than BAL (Fig. 7b).
Most of the sEMG amplitudes of the knee-extensor muscles decreased with fatigue. Figure 4a shows the reduction of A VM at ISO contractions to 63.5 ± 19.9% . At B40, A VM was reduced to 72.4 ± 25.3% . A VL decreased to 59.8 ± 15% for ISO and a tendency was found for B40 69.7 ± 23.3% . At ISO, A RF showed a reduction to 53.4 ± 16.2%, and at B40, a tendential reduction to 70.4 ± 37.9% was found. The slopes of the linear regressions of the sEMG amplitudes A VM , A VL , and A RF were not different for the ISO, B40, and B20 contractions (Table 5).

Comparison between optimized parameters and direct values
Our model-based analysis allowed the comparison between optimized parameters and values calculated directly from measured data. Using isometric MVCs at assumed optimal CE lengths, only G(X) determines the discrepancy between the maximum force-generating capacity of the CE ( f iso and/or n max ) and F MAX . Thus, their normalized values should fit (Figs. 4a, 7a).
The values of P MAX and the FV parameter p max can be directly compared (Fig. 4(b)). However, P MAX must be greater than p max at ballistic contractions, because energy is initially stored in the SEE and released later (Komi 1984).
The difference in the values of the measured maximum velocity per contraction ( V MAX ) and the FV parameter maximum velocity ( v max ) arises from the G(X) as well as from the interplay of CE and SEE. Thereby, V MAX represents the transferred sum of the velocity of CE and SEE, whereas v max represents the maximum possible velocity of the CE and even may not be reached during the movement (Fig. 6).
Considering that fiber activation thresholds increase with fatigue (Grabowski et al. 1972) and that calcium sensitivity decreases (Westerblad and Allen 1991), the normalized trends of the values of the rectified and smoothed sEMG and the value of f iso and/or n max are expected to separate with fatigue (Fig. 4a). Thus, sEMG amplitudes of fatigued contractions cannot be easily converted to the active state of activation dynamics, which becomes even more difficult for dynamic contractions (Disselhorst-Klug et al. 2009). These effects also apply to RFD which also highly depends on the interval used for calculation. Thus, neither RFD nor slope values calculated from sEMG amplitude can be compared to the slope parameter of the activation dynamics U (Fig. 5a).

Discussion
Especially in modeling dynamic movements, the geometric configuration of the legs (Bobbert 2012), the dynamics of contractile and elastic elements (Penasso and Thaller 2017), as well as the state of muscular activation (Siebert et al. 2007) determine the resulting movement (Bernstein 1967). For example, the external maximum power ( P MAX ) must be greater than the parameter maximum power ( p max ) of the CE, because the externally measured values represent the combined dynamics of the CE and SEE. Thereby, a parameter is a property of a structure, e.g., a muscle, which does not depend on the measurement condition or any other variable. Thus, movement-dependent values, e.g., RFD which is different for ISO and BAL contractions (Table 4), should not be used as model parameters.
Parameters reflecting the geometric relations can be measured directly and may be used within geometric functions to calculate internal forces. As the geometric function incorporates the knee as a hinge joint with constant lever arm, a more detailed description of the knee, e.g., Van Eijden et al. (1986), would change the identified values. The resulting parameter trends would be shifted. However, their trends would follow the same patterns because the range of motion was constant. Using parameter optimization methods, these internal forces can be reconstructed by a model of force generation, e.g., a Hill-type muscle model. These are characterized by parameter values that, within the range of validity of the model, do not depend on the movement modality. The curvature of the muscle FV is such a parameter and can be determined using the modeling approach. It is related to the individual fiber-type distribution (Sust et al. 1997b;MacIntosh et al. 1993), which cannot be extracted from simple measurements without an underlying model. For isometric tasks (ISO), the maximum value of the first three contractions and for the ballistic contractions with 40 kg (B40) and 20 kg (B20) extra load, the maximum of the first two contractions per subject is compared to the value calculated from the last contraction. The parameter description is given in Table 2 Significance and effect size of time series analysis and differences between linear slopes for ISO, B40, and B20 tasks of max. force per contraction ( F MAX ), max.  Table 5 Fresh   To determine the internal parameter values within physiological limits, a constrained parameter identification routine was used to optimize the parameters to fit the measured knee-extension data of fresh and continuously fatigued subjects. This allowed us to establish fatigue-related model parameter trends of A(t), FV, and SEE, which also revealed task-specific effects of isometrically induced fatigue regarding ISO and BAL tasks.
The analysis of the linear trends showed that most of the externally measured values declined more for ISO compared to the subsequent BAL tasks. In general, isometrically induced fatigue effected slow high-force tasks more than fast low-force tasks (Tables 3,4,5). Similar results were already reported by Morel et al. (2015), who showed that the magnitude of muscle activation was more affected at ISO tasks than at fast dynamic tasks, which seemed to be mainly affected by peripheral fatigue factors. As the parameters of A(t) also showed such a behavior, our findings could be explained by similar mechanisms. Thereby, the rate constant of muscle activation (U) and the maximum capacity of force generation of the CE (related to the number of maximally available fibers n max ) declined at ISO only (Fig. 5a). Thus, it is likely that the BAL tasks were mainly affected by fatigue-related changes in the peripheral system (e.g., excitation-contraction coupling, cross-bridge dynamics), whereas the ISO tasks had been influenced by a combination of central and peripheral factors.
At very low load dynamic movements performed with maximum effort, it is possible that higher-threshold motor units are not recruited (Milner-Brown et al. 1973). Furthermore, if a dynamic movement lasts less than the time needed for full activation, only a part of the whole amount of motor units is activated. Thus, it is likely that several higher-threshold motor units are not recruited.
The fatigue effects on the CE resulted in an increase in the FV maximum efficiency ( max ), while the maximum contraction velocity ( v max ) was reduced. As the weaker lowerthreshold motor units (slow-type I) muscle fibers also have a higher value of max (MacIntosh et al. 1993;Sust et al. 1997b), this shift can be explained by the fatigue-related deactivation of high-threshold motor units (fast-type II) (Fitts 2008;Johnson et al. 1973;Westad et al. 2003). Another explanation is the single fiber force reduction due to inorganic phosphate accumulation (Stutzig et al. 2017).
Together with the active elements, the passive SEE became less stiff ( 74 ± 7.5% , Fig. 5b). A similar increase in the compliance of the SEE was reported by an ultrasound study of Kubo et al. (2001), who showed that the m. vastus lateralis tendon and aponeurosis stiffness decreased to 72% after 50 sets consisting of 3 s MVC and 3 s relaxation. Furthermore, their results indicated that the increase in compliance was mainly related to the duration of the contraction and not to the functioning of the muscle.

Limitations
Our knee-extensor muscle model represents the combined function of all muscles involved and thus does not account for synergistic and/or antagonistic force sharing (Stutzig and Siebert 2015;Herzog and Leonard 1991). Subject-and sport-related adaptations in force-length relations were not considered , thus, such effects are absorbed by the optimized parameters, but should be negligible compared to the induced fatigue effects. As the model represents a more general mechanical level, cellular mechanisms (Allen et al. 2008;Fitts 1994;MacIntosh et al. 2012;Place et al. 2010;Rzanny et al. 2016), as well as effects of non-constant recruitment thresholds of muscle fiber types (Bigland-Ritchie and Woods 1984), were not modeled. Central and peripheral fatigue effects were considered in the data analysis, but an underlying model that could explain the task dependency of fatigue was not included.
Extracting model parameter values from maximum voluntary contractions is always subject to uncertainties related to the achieved magnitude of effort. For example, the first ISO was often performed more cautiously than the subsequent ones, where one of the second to fourth contractions showed the highest force and/or rate of force development (Fig. 7a). Previous works had shown that the model is able to simulate the knee-extension movement if the model requirements are met (Sust et al. 1997b;Siebert et al. 2008). By comparing a similar modeling approach and the electroencephalography mean spectra theta amplitude Sust et al. (1997a) showed that poor force fits (33% of their data) seemed to coincide with sub-maximal contractions.
As the solution space of our parameter identification routine is relatively flat, only one contraction, performed submaximally or with non-constant effort, is able to create false attractors for the solver. This would lead to unstable and unphysiological results that are excluded using bound constraints. Especially when fatigue sets in, the physical performance highly depends on the individual perception of fatigue (Noakes et al. 2005). Subjects might get used to repeated constant encouragement with the repetitions and pain increases. Thus, some parameter identification problems may be related to short-term motivational problems impairing the subjects' ability to reliably perform MVCs at fatigued states (MacIntosh et al. 1993). As the parameter identification was again successful after it had failed at a level of fatigue, it is reasonable to explain the 12.1% replaced datasets with unsuccessful MVC attempts. The amount of 40% excluded subjects seems not drastically higher than the reported 33% value of Sust et al. (1997a) if one considers that the strict protocol did not allow second MVC attempts and MVCs had to be performed under more difficult conditions. Thus, only data of subjects who, despite fatigue sensations, were able to very reliably perform MVCs could be included. Consequently, the reported values characterize only highly motivated active young healthy males. Reducing the amount of data loss by specifically training the subjects (Maffiuletti et al. 2016) and by using motivating interventions such as rewards would clearly improve the approach. However, it is known that the fresh parameter values vary individually and between different populations (Thaller and Wagner 2004;Wagner et al. 2006). Thus, in dependence of the initial state, the magnitude of the fatigue response should vary individually and between populations. The general behavior, however, should be the same.
History effects of contractions are known since the early 1940s (Abbott and Aubert 1952). Force depression, denoting reduced isometric force, occurs if an isometric contraction is performed after a shortening contraction against resistance (Siebert et al. 2008), but disappears if shortening continues (Till et al. 2014). Force enhancement describes the opposite, where following a lengthening contraction the isometric force is increased (Herzog 1998). These phenomena can be explained by calcium-controlled mechanisms where titin binds to the actin filament (Herzog 2014;Rode et al. 2009). Furthermore, calcium sensitivity of the myosin light chain is increased by preceding contractions, thus improving muscular performance (postactivation potentiation) (Sale 2002). As our protocol included isometric and ballistic contractions prior to further contractions, the effects of postactivation potentiation, fatigue, and to some extend also force depression overlap and might compensate and/or amplify each other. Consequently, our identified parameter values are influenced by these basic mechanisms. Only more detailed cross-bridge models would be able to consider such and further dependencies, e.g., length-dependent activation, but require more parameters. This, however, facilitates that the model output is only barely sensitive to some parameters (Rockenfeller et al. 2015) and as a consequence parameter identification is not applicable to these hardly measurable parameters.
To assess the effect of the rapid central and peripheral recovery within the first 2 min (Froyd et al. 2013;Gruet et al. 2014), we tested if the short recovery periods prior to ISO tasks ( ∼ 3 s) and to the B40 and B20 tasks ( ∼ 5 s) had substantial influence on the contractile properties of the muscles. A Wilcoxon signed-rank test showed that neither RFD nor F MAX of the ninth ISO per set was different to the first ISO of the subsequent set. Thus, the subjects' muscles seem not to have recovered significantly during the ballistic contractions. However, it is important to note that the subjects always looked forward to the BAL contractions because these where less painful than the ISO contractions in the fatigue conditions. This motivational boost in favor of BAL certainly contributed to the observed differences between ISO vs. BAL, but less to the differences between B40 vs. B20 (cf. Tables 4, 5). Furthermore, also the strict sequence of contractions (ISO-B40-B20) might have influenced the results. Despite these limitations, our results reflect how a Hill-type muscle model with activation dynamics should typically adapt to fatigue. We have demonstrated that muscle activation is crucially related to task dependency and that CE as well as SEE parameters change with fatigue.
Even though our parameter identification approach is limited to contractions performed with maximum effort, widespread simulation models of fresh movements, which were also validated in the submaximal case (e.g., OpenSim), principally rely on parameter values similar to ours. In contrast to models using parameter values of fresh muscles, for modeling fatigued movements the recovery process must also be taken into account. In long-duration submaximal movements, inactive fibers recover or fresh ones are activated additionally, thus changing the fatigue state of the whole muscle. However, in short submaximal movements, these effects may be neglected and therefore the identified muscle parameters may be used for simulations.

Conclusions
Using a modeling and parameter identification approach, it was possible to determine fatigue-related trends of Hill-type muscle model parameters and activation dynamics from simple measurements. The use of the mathematical muscle model allowed us to study the effects of fatigue at a macroscopic level. Therefore, our results can be used within other Hill-type muscle models for short-duration movements.
According to our results, models intended for simulating movements of fatigued subjects should consider fatiguerelated parameter dependencies in the sub-models of contractile element, and serial elastic element. Additionally, the muscle activation dynamics should account for contractiontype dependencies if simulations of mixed task modalities are intended. To reveal the underlying neurological control of the task dependency, further studies would be beneficial for creating detailed fatigue models. tutional and/or national research committee and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards.

Informed consent Informed consent was obtained from all individual participants included in the study.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.