Determination of individual knee-extensor properties from leg extensions and parameter identification

ABSTRACT Neural commands control skeletal muscles that act on passive structures and regulate voluntary movements. Mathematical models can simulate such movements and therefore require the knowledge of neuromuscular properties. In contrast to scaling these properties to the individual, we present a non-linear parameter identification method to determine them non-invasively and in vivo. The classic model A describes an excitable contractile element (CE) embedded in a geometrical representation of the leg. Its extension model B is used to study the effects of the force–length relationship and the serial elastic element (SEE). We show the validation of model B and discuss the quality of neuromuscular properties identified from simulations and experiments. The main finding is that identifications that consider CE–SEE dynamics result in increased and more realistic curvatures of the force–velocity relations. This shows that CE and SEE work interdependently and we recommend to co-ordinate the parameter values of muscle–tendon units.


Introduction
In the last 20 years, mathematical models of human movements have been used with increasing frequency [1]. They contain several joints or just one joint, many muscles or a few muscles, sometimes even only one 'model muscle' that represents a certain muscle group (e.g. [2][3][4][5][6][7][8][9]). Models help to plan efficient exercise programmes [10], help to optimize sports equipment [11] and are used to calculate internal loads [7,12]. During movements, there are many non-linear interactions between active and passive structures of the body. In simplified terms, contractile elements (CEs) are recruited by activation dynamics and generate force, thus lengthening the linked elastic elements. The forces of those muscle-tendon complexes then move the segments of the body.
Every solution of a mathematical model depends on input parameters. For simulating human motion, it is therefore important to use proper values for the input parameters. Especially if explosive movements are investigated, the interaction between muscles and tendons becomes important. Energy can be stored in the elastic elements and is released later [1,13]. These interactions influence the resulting movement crucially [1,[13][14][15][16]. However, muscle-tendon parameters are still very difficult to measure in vivo [17][18][19]. To find average values for the input parameters, there exist many different ways. For the CEs, single fibre measurements are often scaled to the subjects' height, weight, sex, age or muscle diameter [19]. Muscle-tendon input parameters are also obtained from measurements on cadavers or are found using ultrasound [18,20]. All these methods have in common that they may depend Furthermore, most of them do not reflect the great interindividual differences of the parameter values [21,22]. For example, it was shown that Hill's parameter b varied by ±28% interindividually. Furthermore, the value of b normalized to muscle length predicted fibre distribution (r ¼ 0:92) [23]. Also, connections between individual values of parameter combinations and angiotensin-converting enzyme gene polymorphism were established [24]. To overcome the limitations of measuring parameter values one by one, parameter identification routines can be used to estimate whole parameter sets simultaneously and therefore are able to represent the dynamic functioning of the subsystems modelled. It is important to mention that the parameters we consider for parameter identification (e.g. parameters of CE and serial elastic element [SEE]) are simultaneously identified from different conditions of movement and thus are movement-independent constants. Nevertheless, models that include many nonlinear interactions between structures (e.g. many muscles and joints) can become very complex. Used in parameter identification routines, such models are more likely to cause principle problems than simpler models (e.g. [25,26]). Thus, we rely on a simple model to increase the possibility for getting a unique solution of the system, even if it is applied to noisy, measured data.
The aim of this study is to extend a simple muscle model A without elasticities and length dependencies [6,21,23,24,27] to a more detailed model B with force-length characteristics and serial tendon elasticity. Both models are used to identify individual neuromuscular parameter values from maximum voluntary contractions (MVCs) non-invasively and in vivo. By adding force-length characteristics and serial elasticities, we indent to extract more realistic parameter values and thus aim to enhance this type of analysis of athletic performance. Furthermore, the comparison of both model versions will highlight the impact of force-length characteristics and serial elasticities on the identified properties.
The paper is organized as follows: in Section 2 we describe the measurement device and the measurements and give detailed information on the equations of motion of model A and model B. The process of parameter identification as well as the verification and validation methods is described. The literature values that are used in the model B configurations are shown, and the statistical methods for comparing the parameter identification results from models A and B are expressed. Section 3 presents results of the parameter identifications, simulations and comparisons of the models. It provides statistical data on the goodness of fit and the variations observed. In Section 4 we discuss if we can identify the muscle-tendon properties from experimental data sufficiently. To show that model B is an improvement of model A, the additionally identified parameters of the model B configurations (e.g. the properties of elastic elements) have to lie in the range of measured values from literature. Furthermore, model B simulations have to predict phenomena observed in reality better than model A.

Methods
All numerical calculations are performed in MathWorks Matlab®2015b. Symbolical computations are done using Wolfram Mathematica 10.3.

Measurements
We included 14 male students of sports science who voluntarily participated in the study (24.7 ± 3.8 years, 1.83 ± 0.07 m, 75.5 ± 5.9 kg and BMI 22.5 ± 1.5). They performed four (2 × 2) dynamic MVCs at an inclined leg press against movable masses of 52.5 and 72.5 kg, respectively. These were followed by three isometric MVCs. All contractions started at a kneeextension angle of 120°. Dynamic contractions were repeated twice for each load and subsequently three isometric contractions were completed. The duration of the break in between the contractions was chosen by each subject individually. They were advised to wait until full recovery. Each subject was familiarized with the procedure in the week prior to the test. The study was approved by the local ethics committee of the University of Graz.

Data acquisition
The force F meas was measured at an inclined leg press (γ ¼ 10 ) (Tetra® Ilmenau) with a Kistler® force plate (500 Hz). The force plate is mounted onto a sledge, and its position X is recorded by a high-frequency position sensor (500 Hz). Data from a 3D motion capture system (Vicon M2 cameras) with markers on both legs (250 Hz) were synchronized. We placed the markers at the greater trochanter, the lateral femoral epicondyle and the lateral malleolus. The anthropometric lengths of the right leg were manually measured with an accuracy of ±5 mm using a tape measure. The length of the thigh l t was defined as the distance from the palpated greater trochanter to the lateral knee joint line. The length of the shank l s was measured from the lateral knee joint line to the lateral malleolus. The circumference of the knee CIR and the patellar tendon length l p were obtained at the fully relaxed and extended knee. CIR was measured at the level of the centre of the patella and l p was defined as the distance from the centre of the patella to the tibial tuberosity (cf. Figure 1(a)) [23,24]. At the beginning of the measurement X was shifted to match l t þ l s when the knees were fully extended. A position mark was set via the cosine rule to enforce an initial kneeextension angle of 120°at each contraction.

Post-processing
Force and position data were low-pass filtered at 15 and 30 Hz (zero-lag, butterworth, 4. O) for isometric and dynamic movements, respectively.
The comparison of the knee angles calculated from motion capture data and position sensor data revealed displacements of the subjects within the seat. During isometric contractions, the (a) ( b) Figure 1. Illustration of models A and B within the inclined leg press. (a) Shows the modelled leg and muscle-tendon complex within the leg press. (b) Points out the structural differences between model A and model B. Model A has no length dependencies and thus requires no CE length. Model B uses l CEopt ¼ 0:09 m [5,8] (reached at α opt ¼ 120 [4]), and l SEE;slack ¼ l p (reached at the fully extended knee). (a) Geometrical relations (proportions are not in scale) of the leg within the leg press. The angles in the sketch are: γ the inclination of the leg press, α the knee-extension angle, β the angle between the force vector of the modelled muscle-tendon complex and the knee moment arm r, β′ the angle between r and the patellar tendon.
The forces are: f int the force of the muscle-tendon complex and F ext the external force at the force plate. The lengths are: X the distance from the greater trochanter to the lateral malleolus, l t the length of the thigh, l s the length of the shank, l MTC the length of the muscle-tendon complex, l p the distance from the centre of the patella to the tibial tuberosity,l p the distance from the centre of the knee to the tibial tuberosity. (b) Schematic representation of the muscle-tendon complex f int of models A and B. f int of model A consists of a contractile element CE. Model B is an extension to model A adding a serial elastic element SEE and a force-length relation.
subjects deformed the device and soft tissue and reached a knee angle of 124°on average. At the end of the dynamic contractions, the inertia of the legs slightly pulled the fixed subjects out of the seat and motion capture data revealed that the feet lost contact with the force plate at 161°on average. As a consequence, X was corrected in each condition with respect to the average angle. Furthermore, force data of dynamic contractions were corrected to exactly match F meas ð0Þ ¼ Àmg sinðγÞ, using g ¼ À9:81 ms-2.
To extract single contractions from the continuous data file, force onsets were manually determined. The endpoints of isometric contractions were automatically set when the force dropped below 97% of the maximum of the contraction. As the parameter identification routine was sensible to the very last milliseconds of the dynamic contractions (maybe due to effects of contact elasticities), their endpoints were automatically set when the force dropped below the value at onset. The maximum velocity had already been reached at this point.
Submaximal isometric contractions were detected by comparing the rate of force development RFD from dynamic and isometric conditions [28], where RFD ISO had to be greater than RFD DYN . Within dynamic contractions the RFD was compared visually for equal loads and only the contraction with the greatest RFD per condition was kept. Thus, each final data set per subject included three contractions from different conditions. If at least one of the RFD DYN was greater than the steepest RFD ISO , the subject was excluded.

Forward simulation of the models of movement
To simulate the experiments, we used two mathematical models. Both models A and B describe the concentric or isometric leg-extension task with ordinary differential equations (ODEs) written as equations of motion. The numerical solution of these non-linear ODEs is calculated with ODE45: ('MaxStep', 0.005, 'RelTol', 1e-8, 'Refine', 1, 'Events', @events). Both models use a geometrical transfer function GðXÞ to convert the internal force f int to the external force F ext (Equation (2)) [29]. Thereby, the individual measures of l t , l s , the circumference of the knee CIR ¼ 2rπ and l p are used to calculate leg and knee moments ( Figure 1). F ext moves the distally located point-mass which represents the sledge moving in the gravitational field [6]: whereg ¼ À sinðγÞ Á g. Models A and B differ in the definition of f int only (cf. Figure 1(b)). GðXÞ is identical for both models and is defined as the ratio between F ext and f int . Both F ext and f int act on the moment of the knee. Assuming that the direction of r is moving (simplified patella movement) such that β ¼ β 0 (defined in Figure 1(a)), and using further geometrical simplifications, GðXÞ is given by [6,21,23,24,29]: The knee-extension angle α can be expressed as a function of β as well: where l MTC is the length of the muscle-tendon complex, andl p is the projection of l p onto the shank. As the lengths of l MTC andl p cannot be measured externally, we use the following anthropometric simplifications: l MTC ¼ l t andl p ¼ l p (cf. Figure 1(a)). GðXÞ is determined with the following numerical calculation: The function of muscle activation AðtÞ is identical in both models and is modelled with two coupled ODEs that represent the excitation uðtÞ and recruitment nðtÞ of inactive muscle fibres dnðtÞ dt ¼ uðtÞ Á ðn max À nðtÞÞ; where n max is the maximum number of available fibres and u max is the maximum excitation. This system is solved by using the boundary conditions nð0Þ ¼ A pre and uð0Þ ¼ 0, where A pre is the pre-activation level [30]. As AðtÞ is a ratio between 0 and 1, we set n max ¼ 1 (100%). To reduce the number of parameters, we set u max ¼ U. Finally, only U characterizes the slope of the function [6]: where t 0 is the onset of muscle contraction (cf. Figure 2).

Model A
The original model was introduced by Sust [31] and an improved version is described in [6]. This model is slightly modified with a smooth onset of AðtÞ (Equation (7)) and is denoted as model A. f int;A is equal to the force of the muscle f M A which is represented by the force-velocity relation f CE ðv M A Þ of Hill [32] coupled with AðtÞ: f CE ðv M A Þ depends on the contraction velocity of the muscle v M A . Its shape is characterized by a set of three subject-specific Hill parameters: a (N), b (ms -1 ) and c (W), that are convertible to the maximum isometric muscle force, or other equivalent sets of parameters [27]. To calculate AðtÞ, A pre;A has to be calculated initially. Therefore, the equation of motion (Equation (1)) together with Equation (8) is rearranged, and GðX 0 Þ and f CE ð0Þ ¼ f iso are used: where X 0 ¼ Xð0Þ. Using the assumption that the internally produced mechanical energy is losslessly transferred to the environment, we get [6]: This is used to get a formulation of the equation of motion that can be solved with ODE45 [6]: This forward simulation of model A stops at the condition where v M A ! v max . The solver returns X and _ X, and thus the simulated external force of model A, F mod;A , as well as the behaviour of the substructures can be calculated. The time evolution of GðXÞ is determined from the resulting X and the known individual anthropometric lengths. AðtÞ is calculated using the known t 0 and A pre (Equation (9)). With a, b, c and, _ X we get f int;A .

Model B
Model B is an extension to model A adding an SEE, and a force-length relation FLðl M B Þ ( Figure 1(b)).
where f int;B represents the force of the whole muscle-tendon complex f MTC . Except for v M , Figure 3(a), Equation (11)). GðXÞ is equivalent to Equation (2) and AðtÞ is identical to Equation (7).
The SEE is expressed within the contraction velocity of the muscle v M B : The velocity of the SEE, _ l SEE , is derived from the function of the tendon force that includes a nonlinear toe region. The non-linear part is created by a force dependency of the stiffness of the SEE, k SEE ðf MTC Þ (Figure 3(b)): The rearranged derivative of Equation (14) gives _ l SEE where we use that f SEE ¼ f MTC (Equation (A8)). After substituting k SEE and _ k SEE (Equation (A6)), f MTC and _ f MTC can be eliminated by using Equation (A4). Finally, we obtain where _ GðXÞ is the analytic derivative of GðXÞ. The detailed formulation of the non-linear toe region is described in Appendix 1.1.
The force-length relation FLðl M B Þ depends on the length of the muscle l M B and is identical to At the initial state of the simulation, A pre;B and the initial length of the muscle l M B ;0 must be calculated.
A pre;B is calculated by rearranging the equation of motion (Equation (12)): l M B ;0 is calculated using Equations (A3) and (A4), the initial displacement of the SEE, Δl SEE;0 , and the patellar tendon slack length l SEE;slack : The SEE in our model is a representation of all SEEs reaching from the origin of the MTC to the segment representing the tibiae. As there is no procedure to determine the correct, merged properties of all SEEs of the quadriceps, we set l SEE;slack ¼ l p . Further interactions between model B elements are described in Appendix 1.3 (cf. Figure 1(b)).
To get a solvable formulation of the equation of motion (Equation (12)), we use Equation (15) in Equation (13) and rearrange the merged equations for X. After the solver terminates at either α > 179 or m € X þ mg ¼ 0, it returns X, _ X and € X. Using these, the time evolution of each substructure is determined. GðXÞ, _ GðXÞ and l M B are calculated using X, _ X and the individual anthropometric lengths. From Equation (A4), f MTC is known and we get Δl SEE from Equation (A3). From these, FLðl M B Þ, l SEE and l CE are calculated. Post hoc, we can

Parameter identification
Contrary to the forward simulation, the observed movement is already known. Thus, the parameter identification routine uses X and F meas as inputs and minimizes the weighted sum of squared errors of the identification model output, F opt vs. F meas (cf. Figures 4 and 5). It processes data sets containing isometric and dynamic contractions from different loads at the same time.
The mass of the empty sledge including the force plate m sledge was determined using m sledge ¼ F meas Ág À1 ¼ 32:5 kg. One or two additional weights m add ¼ 20 kg were attached during the experiments. Thus, the overall load was m ¼ m sledge þ i Á m add , i ¼ 1; 2. The identification process uses and As multiple curves with different numbers of datapoints are merged into one data set, the influence of the number of datapoints must be removed from the identification process. This is represented by the weight factor W, which scales the sum of errors of each curve to the sum of errors of the first curve loaded. This is applied to the cost function of the solver: min bndcon2½Par min ;Par max ; nonlconðParÞ 0 To improve the computational speed for calculating β, we solve Equation (3) for β 2 ½0:5; 1:4 at an increment of 5 Á 10 À5 and fit a piecewise cubic Hermite interpolating polynomial to the result. The resulting function works as a look-up function to get βðαÞ within α 2 ½100 ; 180 . With these arrangements, the GlobalSearch class is initialized using the constraints, the objective function of the model and the fmincon interior point algorithm ('StartPointsToRun', 'bounds-ineqs', 'TolFun', 1e-18, 'TolX', 1e-18, 'NumStageOnePoints', 600, 'NumTrialPoints', 3000). Figure 4 illustrates the performed steps.
After identifying the set of parameters, the routine warns if the solver returns values close to one of the bound constraints (Equation (22)). Affected runs were repeated with a new set of the automatically randomized starting points. If the warning remained, the subject was excluded.

Model
A _ X is calculated numerically by using the derivative of a cubic spline fit of the filtered data of the measured X. Using X, the individual anthropometric lengths, the numeric solution of βðαÞ and GðXÞ are determined. With Equation (11), € Xð0Þ ¼ 0 and v M A ð0Þ ¼ 0, we calculate A pre;A for each curve loaded: Using the relative time, t, A pre;A and two solver parameters U and t 0 , AðtÞ can be evaluated at every solver iteration. According to Equations (10) and (11), f CE ðv M A Þ is parameterized with f iso , v max and η max . Thus, F opt;A is given by 2.3.2. Model B F opt;B depends on X; _ X; € X and € X _ (cf. Equations (12)(13)(14)(15)). To calculate € X we use Equation (20): Using A pre;B , the time vector, as well as t 0 and U from the solver, AðtÞ is determined. f int;B is parameterized with f iso , v max , η max and k SEE , and if enabled n SEE;nl is identified too (cf. Equation (12)). We obtain F opt;B at every solver step: Finally, F opt;B is used together with F meas to minimize the cost function (Equation (21)).
To test the added substructures of model B, we run the identification procedure with different configurations ( Table 1). The configurations use seven to nine degrees of freedom dof and/or fixed values from the literature. To test the effect of the additional dof on the uniqueness of the optimization problem we first simulate data and then re-identify the parameters (cf. Section 2.5).

Confidence interval and correlations between identified parameters
After the identification-solver terminates, the goodness of the fit (F opt vs. F meas ) is assessed by the adjusted coefficient of determination r 2 a . Furthermore, the confidence interval, the parameter correlation matrix and the error propagation are computed. The statistics of the final set of parameters are based on [33,34] and [35, pp. 257-262]. The details of these calculations are described in Appendix 2.

Subsequent forward simulation
To ensure that not only the measured movement determines the identified parameters but also a simulation using these parameters results in the same motion, we generate a simulated data set containing each condition of the movement measured, e.g. isometric and dynamic contractions with different loads (cf. Figure 5). The measured and simulated data sets are compared with the r 2 a of F meas vs. F sim and this measure is denoted as r 2 a sim. As the simulation has a different time vector than the measurement, we fit a piecewise cubic Hermite interpolating polynomial to the simulated data and evaluate this function with the time vector of the measurement. This evaluation is done for the overlapping period only.

Manipulating B to reproduce A
To check if the added model B elements work as expected, we set their parameters to values that cause their ineffectivity. The influence of the FLðl M B Þ is disabled by setting WIDTH ¼ inf and thus FLðl M B Þ ¼ 1 (Equation (A7)). Theoretically, the SEE is ineffective if it is indefinitely stiff. Setting k SEE;lin ¼ inf causes too long computation times and a suitable threshold of k SEE;lin has to be found (Equation (14)). To quantify this value, the normalized mean squared error (NMSE) of F mod;A vs. F mod;B is calculated. Model A and model B outputs are judged to be equal if NMSE > 0:9999.

Parameter re-identification from simulated data sets
To show that the identification routine is correct and which parameters can be identified at which quality, we re-identify the input parameters of simulated data sets (cf. 'In' rows of Table 2). The consistency of the identification routine is assessed for each parameter by using the relative error δ ¼ Param Out =Param In À 1 and by the coefficient of variation c v ¼ σ Param id f g =Param id (Equation (B4)).

Comparison of parameters identified with model A and model B using measured data sets
From all subjects we calculate the average μ and the standard deviation σ of the parameters identified. The reported c v s are the averages of all subjects. Statistics are computed in MathWorks® Matlab 2015b. We compare the values of the final parameters from different models and their configurations by using a one-way ANOVA which is followed by a post-hoc pairwise comparison (Tukey HSD, if the main effect is 0:05). The effect size is determined with Hedges' g, g Ã [36] and is reported only if the Tukey HSD is significant. Effects in the range 0:2 g Ã <0:5 are classified as small, (y), 0:5 g Ã <0:8 as medium (yy) and g Ã ! 0:8 as large (yyy) [37]. Normality of data sets is assessed by using the Shapiro-Wilk test. If normality is violated, a Wilcoxon rank sum test is conducted prior to the ANOVA. The ANOVA results are reported only if the Wilcoxon rank sum test is also significant (p 0:05).   [3,13,45]. We set TH % to f MTC;TH ¼ 0:4 Á f iso . The toe region is usually modelled as a power function and is set to the power of 2.5 (n SEE;nl ¼ 2:5) [46].

Identifying SEE properties (model B k id )
To identify properties of the SEE, the number of parameters to identify is increased. We use a linear and a non-linear configuration of the SEE and identify k SEE;lin (cf. Table 1). The linear configuration is parameterized by using TH % ¼ 1 (100%) and n SEE;nl ¼ 1. In model B k id n id the number of parameters to identify increases again, and k SEE;lin is identified together with n SEE;nl .

The influence of the force-length relation
Finally, we investigate the influence of the FLðl M B Þ by comparing the identified parameter values from models B k id with those of the corresponding configurations B ? k id where the FLðl M B Þ is disabled (cf. Section 2.4).

Accuracy of identified parameters
The mean accuracy μ acc f g of a parameter is calculated via the mean coefficients of variation μ c v f g in the cohort: μ acc f g ¼ μ Á μ c v f g, where μ is the mean parameter value. To assess the detectability of the individual differences in each parameter we use μ acc f g with respect to the interindividual range (Par max À Par min ): μ acc % f g¼ μ acc f g Par max À Par min ¼ μ Á μ c v f g Par max À Par min : To avoid effects of outliers, we excluded data outside 99.3% of the distribution.

Parameter re-identification from simulated data sets
The parameter re-identification process (cf. Figure 4, dotted arrows) was performed on simulated data of model B. The input parameters of the simulations were set to the values of the 'In' rows of Table 2. To re-identify these, we used the first three configurations of model B listed in Table 1.
The re-identifications of the input parameters led to an r 2 a sim ¼ 1. When simulations using n SEE;nl ¼ 2:5 were re-identified (upper section of Table 2), we had to cut the simulated data where F mod;B dropped below its initial value to get the correct values returned. To re-identify simulations using n SEE;nl ¼ 1 with model B k id n id , we had to use the same procedure (bottom row of Table 2).
Model B k id n in¼1 returned the input values with the lowest mean relative error μ δ f g ¼ 6:05 Á 10 À6 . However, this configuration showed the greatest mean coefficient of variation μ c v f g ¼ 0:46. The lowest μ c v f g ¼ 0:15 was found for the non-linear SEE configuration, model B k id n in¼2:5 which showed a much higher μ δ f g ¼ 213:46 Á 10 À6 (cf. Table 2).  (Table 3). Relative errors were not different either (Table 4).

Identifying SEE properties (model B k id )
The tests with model B k id revealed that the parameter identification with model B k id n 2:5 resulted in high variations of k SEE;lin , whereas the linear SEE configuration model B k id n 1 showed a much lower c v for k SEE;lin (Tables 3 and 4). Among the c v s of all parameters, model B k id n 1 exhibited the highest μ c v f g ¼ 0:26 that was still acceptable. The comparison of the correlation matrices of those two model configurations showed that model B k id n 2:5 produced higher parameter correlations between the activation onsets (t 0;20 , t 0;40 and t 0;ISO ) and the other parameters. Except from a 4.7 times higher value for corr(v max ; k SEE;lin ), the muscle-tendon complex parameters exhibited lower correlations.
Using model B k id n id on measured data sets resulted in problems with the bound constraints. The solver tried to converge towards n SEE;nl < 1 and caused CE lengths beyond the physiological range and thus led to errors.
As the high μ c v f g of k SEE;lin in model B k id n 2:5 cannot be accepted, only the results of model B k id n 1 will be discussed below.

3.3.2.1.
Identified SEE properties with model B k id n 1 . The comparison of parameters identified with model A and model B k id n 1 revealed differences in p max and in curvature-related parameter values η max and a/f iso (cf. Figure 6, Table 3). The values of the Hill-constants a ÃÃÃ yyy , b ÃÃ yyy and c ÃÃÃ yyy * p 0:05, ** p 0:01, yyy g * ! 0:8. Models B are compared to A, and models B ? id are compared to B id : Figure 6. Boxplots of identified a=f iso from measured data sets and different models. Note that only models B kidn1 and B ? kidn1 are significantly different from the others. differed as well. The c v increased for f iso and U (Table 4). Table 5 illustrates the changes in the correlation matrix from model B k id n 1 relative to model A. The correlation matrix changed for all parameter combinations except for corr(v max ; f iso ), corr(v max ; U), corr(v max ; t 0;20 ), corr(v max ; t 0;40 ) and corr(f iso ; t 0;40 ).
The identification of k SEE;lin led to a mean value of k SEE;lin ¼ 1777 AE 324 N/mm (c v ¼ 0:3 AE 0:36). This is the average of all subjects and represents both legs.

The influence of the force-length relation
Switching the FLðl M B Þ off (models B ? id ) did not significantly change the values of the identified parameters whereas the c v of f iso returned to the lower value of model A (Tables 3 and 4

Accuracy of identified parameters
The assessed mean accuracy μ acc f g of parameters identified with model A and model B k id n 1 showed that the μ acc f g was smaller than the interindividual range (cf. Table 6).

Discussion
We used non-linear parameter identification strategies to identify parameters of the knee-extensor muscle-tendon complex MTC and investigated how a serial elastic element SEE would change the values of the parameters identified. Therefore, the identification results of the original model A were compared to the introduced Hill-type muscle model B. In model A the 'MTC' is defined by the force-velocity relation whereas model B adds a force-length relation and a serial elastic  Table 6. Parameter range (Par max -Par min without outliers), mean parameter accuracies μ acc f g and μ acc % f g defined as μ acc f g normalized to the range (cf. Equation (29)).
To test the identification procedure we used data sets generated from simulations and measurements. Each data set included isometric as well as concentric dynamic MVCs which were conducted bilaterally at an inclined leg press. This motion was modelled with ODEs representing the equations of motion of model A and different configurations of model B.
The principal problems of complex models (cf. [25,26]) are potentiated if they are used for parameter identification from noisy measured data. If reasonable computation times are needed, limitations in applicability arise. Our model is limited to monoarticular concentric movements with a constant value of the maximum voluntary drive (cf. Section 2.2, n max ). The parameter identification routine requires that the maximum voluntary drive is known for all simultaneously processed input data. Thus, using the assumption that n max ¼ 100%, only MVCs can be used. The model does not describe phenomena such as e.g. force depression, force enhancement [47], tendon hysteresis [48], fatigue/fatigability [49], electromechanical delay [50], moving axis of rotation of the knee [51,52] and mass distributions of segments and muscles [8,53]. Also, contact material deformations (e.g. foot contact elasticity) and force-sharing amongst synergistic and/or antagonistic are not included [54]. Thus, the results of this method represent net monoarticular properties of all muscles involved. Its application to other joints is possible if the function of geometry and the measurement procedure are adopted accordingly [29,55].
Our results showed that model B was able to simulate the dynamics of a contractile element CE that is linked to an SEE. At the beginning of explosive movements, the model stores energy in the SEE and releases it towards the end of the movement [56]. Thus, it is possible that the external power exceeds the maximum power of the CE (cf. Figure 7).
To test the consistency of model A and model B, the parameter values of model B were adjusted accordingly. Model B reproduced model A correctly (cf. Section 3.1).
If we used simulated data sets, the identification process resulted in sufficiently correct reidentified input parameter values. The configuration that was used to re-identify the stiffness of the linear region of the SEE k SEE;lin and the power of the non-linear region n SEE;nl was the most unstable one and had to be discarded (cf. last row in Table 2).
Using measured data sets revealed that k SEE;lin exhibited great variations if the non-linear part of the SEE was enabled in the model. This observed low sensitivity of k SEE;lin of model B k id n 2:5 could be explained by the fact that the MTC reached only about 25% of f iso during the dynamic contractions. To test this assumption, further measurements with additional load and/or an increased inclination of the leg press are needed. As the linear SEE configuration model B k id n 1 showed much less variation with measured data and had the smallest relative error at the reidentification process, it is the best candidate to identify k SEE;lin from measured data (cf . Tables 3,  4, and 5).
To test the influence of the SEE we used a fixed k SEE;lin ¼ 2 Á 3087 N/mm from the literature first [38][39][40][41][42][43][44]. The use of this value had no significant effect on the mean and variation of the identified parameter values. The literature value of k SEE;lin caused changes in the model output that were too small to be detected within the variation between the subjects (cf. rows 2 and 3 in Tables 3 and 4).
The identifications from measured data using model B k id n 1 resulted in an average value of k SEE;lin ¼ 1777 AE 324 N/mm representing the SEEs of both legs (cf. Table 3). As this model configuration describes a linear SEE behaviour, these results do not depend on the overall SEE length. Assuming that the SEEs of two legs are in parallel we get k SEE;lin;single ¼ 888:6 N/mm for a single leg. These are 28.8% of the mean patellar tendon stiffness we found in the literature. However, our result is an overall representation of the serial elastic stiffness of the leg-extensor muscle group. Therefore, the value of k SEE;lin cannot be easily compared to a measured stiffness of a single structure. Furthermore, it is unclear whether the tendon and/or the aponeurosis are in series with the CE [57]. Kubo et al. reported a stiffness of the vastus lateralis aponeurosis of 143 N/ mm (average of both groups) [58]. If we apply this value to all four muscles of the quadriceps we get a value of 572 N/mm for four parallel aponeurosis in series. This gives a range from low aponeurosis stiffness to high tendon stiffness. As our identified SEE stiffness represents the overall serial elastic behaviour, the identified value can indicate whether the aponeurosis or the tendon is more likely to be represented by the SEE model. The low SEE stiffness we identified supports the view that the characteristics of the SEE are more similar to the characteristics of an aponeurosis than to those of a tendon.
Regarding all subjects, comparisons of the identification results of model A and model B k id n 1 (cf. Table 3), which we used to identify k SEE;lin together with the parameters of the CE, showed changes in the maximum power p max and in the curvature related parameter values of the force-velocity relation (maximum efficiency η max and a=f iso ). The maximum contraction velocity v max and the maximum isometric force f iso were not different to model A. However, all Hill parameters a, b and c changed significantly. The maximum efficiency increased from model A (0:27 AE 0:04) to model B k id n 1 (0:38 AE 0:07). It is defined as p max /c, where the Hill constant c (W) represents the energy consumption of the muscle. This relation is considered to be equivalent to the definition of the thermodynamic efficiency η therm ¼ W m =ΔG where W m is the mechanical work and ΔG is the Gibbs free energy. Using this definition, the values of η therm are in the range from 0.45 to 0.68 [59]. Thus, η max from model B k id n 1 moved towards values reported in the literature. As η max is directly related to a=f iso [27], its value has also shifted towards the right direction.
To investigate the influence of the force-length relationship, we compared model B configurations with active and disabled force-length relations. Among the variance of the subjects, we found no significant changes in the values of the parameters identified (cf. Table 3). Only the coefficient of variation of f iso decreased to the level of model A (cf. Table 4).
To be able to distinguish between subjects, the uncertainty of the identified parameters with model A and model B k id n 1 must be smaller than the interindividual range. In Table 6 we see that all average parameter uncertainties are within these bounds. Note that the c v is calculated via the unconstrained behaviour of the model (cf. Appendix 2) and therefore underestimates the accuracy of the constrained parameter identification solver. Additionally, the decrease of accuracy of e.g. f iso in model B k id n 1 does not take into account that the data of the isometric measurement gives more weight to f iso within the solver routine. Thus, the accuracy of the identified results is better than given in Table 6. This can also be seen by comparing σ of model A and B k id n 1 in Table 3 where σ did not change considerably between the models.
Siebert et al. [6,Section 3] used a nearly identical version of model A to identify human neuromuscular properties in vivo, and in the same paper the importance of muscle-tendon dynamics was already pointed out for single muscles [6, Section 2] (and in more detail in [13]). Thus, the combination of both approaches was the next logical step that allowed us to identify the serial elastic stiffness together with the parameters of the CE and the activation dynamics. Considering the discussed above, model B k id n 1 is a feasible extension to model A.

Conclusion
With model B k id n 1 , η max shifted towards the values reported in the literature. Thus, model B k id n 1 identified a more realistic set of parameters than model A and should be favoured if the curvature of the subject-specific force-velocity relation is of interest.
The behaviour of the SEE was mainly influenced by the curvature of the force-velocity relation and vice versa. In more detail we see that the identified curvature of the force-velocity relation is underestimated if muscle-tendon dynamics are neglected. Thus, to achieve realistic MTC dynamics, we recommend to co-ordinate the curvature of the force-velocity relationship with the stiffness of the SEE.
(g) Finally, the confidence interval CI Paramid of each parameter is computed: (2) Parameter correlation matrix: (a) The Hessian matrix H is calculated using hessian.m [33]. Therefore, the analytical function of the model and the final set of identified parameters are used.
(b) With the inverse of H we get the covariance matrix Cov [35, p. 256]: Cov ¼ H À1 : (c) Finally, we compute the parameter correlation matrix Cor: (3) Error propagation: We convert the three parameters that characterize the force-velocity relation to other parameters (cf. [27]. for conversion details). The error propagation is calculated using the law of error propagation for standard deviations that includes also the parameter correlations [60].