Comparative Sensitivity Analysis of Muscle Activation Dynamics

We mathematically compared two models of mammalian striated muscle activation dynamics proposed by Hatze and Zajac. Both models are representative for a broad variety of biomechanical models formulated as ordinary differential equations (ODEs). These models incorporate parameters that directly represent known physiological properties. Other parameters have been introduced to reproduce empirical observations. We used sensitivity analysis to investigate the influence of model parameters on the ODE solutions. In addition, we expanded an existing approach to treating initial conditions as parameters and to calculating second-order sensitivities. Furthermore, we used a global sensitivity analysis approach to include finite ranges of parameter values. Hence, a theoretician striving for model reduction could use the method for identifying particularly low sensitivities to detect superfluous parameters. An experimenter could use it for identifying particularly high sensitivities to improve parameter estimation. Hatze's nonlinear model incorporates some parameters to which activation dynamics is clearly more sensitive than to any parameter in Zajac's linear model. Other than Zajac's model, Hatze's model can, however, reproduce measured shifts in optimal muscle length with varied muscle activity. Accordingly we extracted a specific parameter set for Hatze's model that combines best with a particular muscle force-length relation.


Introduction
Scientific knowledge is gained by an interplay between quantitative real world measurements of physical, chemical, or biological phenomena and the development of mathematical models for understanding the dynamical processes behind. In general, such phenomena are determined as spatiotemporal patterns of physical measures (state variables). Modelling consists of distinguishing the surrounding world from the system that yields the phenomena and formulating a mathematical description of the system, a model, that can predict values of the state variables. The calculations depend on model parameters and often on giving measured input variables. By changing parameter values and analysing the resulting changes in the values of the state variables, the model may then be used as a predictive tool. This way, the model's validity can be verified. If the mathematical model description is moreover derived from first principles, the model has the potential to explain the phenomena in a causal sense.
Calculating the sensitivities of a model's predicted output, that is, the system's state variables, with respect to model parameters is a means of eliminating redundancy and indeterminacy from models and thus helps to identify valid models. Sensitivity analyses can be helpful both in model-based experimental approaches and in purely theoretical work. A modelling theoretician could be looking for parameters to which all state variables are nonsensitive. Such parameters might be superfluous. An experimenter may inspect the model that represents his working hypothesis and analyse which of the model's state variables are specifically sensitive to a selected parameter. Hence, the experimenter would have 2 Computational and Mathematical Methods in Medicine to measure exactly this state variable to identify the value of the selected parameter.
In a biomechanical study Scovil and Ronsky [1] applied sensitivity analysis to examine the dynamics of a mechanical multibody system: a runner's skeleton coupled to muscle activation-contraction dynamics. They calculated specific sensitivity coefficients in three slightly different ways. A sensitivity coefficient is the difference quotient calculated from dividing the change in a state variable by the change in a model parameter value, evaluated in a selected system state [2]. The corresponding partial derivative may be simply called "sensitivity. " Therefore, a sensitivity function is the time evolution of a sensitivity [2]. Accordingly, Lehman and Stark [2] had proposed a more general and unified approach than Scovil and Ronsky [1], which allows systematically calculating the sensitivities of any dynamical system described in terms of ordinary differential equations. As an example for sensitivity functions, Lehman and Stark [2] had applied their proposed method to a muscle-driven model of saccadic eye movement. By calculating a percentage change in a state variable value per percentage change in a parameter value, all sensitivities can be made comprehensively comparable, even across models.
A sensitivity as defined so far is of first order. Methodically, we aim at introducing a step beyond, namely, at calculating second order sensitivities. These measures are suited to quantify how much the sensitivity of a state variable with respect to one model parameter depends on changing another parameter. By analysing second order sensitivities, the strength of their interdependent influence on model dynamics can be determined. In addition to this so-called local sensitivity analysis, we will take the whole parameter variability into account by calculating global sensitivities according to Chan et al. [3] and Saltelli and Chan [4]. This approach allows translating the impact of one parameter on a state variable into a parameter's importance, by completely comprising its interdependent influence in combination with all other parameters' sensitivities.
In this study, we will apply the sensitivity analysis to models that predict how the activity of a muscle (its chemical state) changes when the muscle is stimulated by neural signals (electrical excitation). Such models are used for simulations of muscles' contractions coupled to their activation dynamics. Models for coupled muscular dynamics are often part of neuromusculoskeletal models of biological movement systems. In particular, we want to try and rate two specific model variants of activation dynamics formulated by Zajac [5] and by Hatze [6]. As a first result, we present an example of a simplified version of the Zajac [5] model, in which sensitivity functions can in fact be calculated in closed form. Subsequently we calculate the sensitivities numerically with respect to all model parameters in both models, aiming at an increased understanding of the influence of changes in model parameters on the solutions of the underlying ordinary differential equations (ODEs). Additionally, we discuss which of both models may be physiologically more accurate. The arguments come from a mixture of three different aspects: sensitivity analysis, others' experimental findings, and an additional attempt to best fit different combinations of activation dynamics and force-length relations of the contractile element (CE) in a muscle to known data on shifts in optimal CE length with muscle activity [7].

Two Models for Muscle Activation Dynamics
Macroscopically, a muscle fibre or an assembly thereof, a muscle belly, is often mapped mathematically by a onedimensional massless thread called "contractile component" or "contractile element" (CE) [8][9][10][11][12]. Its absolute length is ℓ CE which may be normalised to the optimal fibre length ℓ CE,opt by ℓ CErel = ℓ CE /ℓ CE,opt . In macroscopic muscle models, the CE muscle force is usually modelled as a function of a force-(CE-)length relation, a force-(CE-)velocity relation, and (CE-)activity . Commonly the muscle activity represents the number of attached cross-bridges within the muscle, normalised to the maximum number available ( 0 ≤ ≤ 1). It can also be considered as the concentration of bound Ca 2+ -ions in the muscle sarcoplasma relative to its physiological maximum. The parameter 0 represents the minimum activity that is assumed to occur without any stimulation [6].
We analyse two different formulations of muscle activation dynamics, that is, the time (its symbol: ) evolution of muscle activity ( ). One formulation of muscle activation dynamics was suggested by Zajac [5], which we modified slightly to take 0 into account: with the initial condition (0) = ,0 . In this context, is supposed to represent the (electrical) stimulation of the muscle, being a parameter for controlling muscle dynamics. It represents the output of the nervous system's dynamics applied to the muscle which in turn interacts with the skeleton, the body mass distribution, the external environment, and therefore the nervous system in a feedback loop. Electromyographic (EMG) signals can be seen as a compound of such neural stimulations collected in a finite volume (being the input to a number of muscle fibres) over a frequency range and coming from a number of (moto-)neurons. The parameter denotes the activation time constant, and = / deact is the ratio of activation to deactivation time constants (deactivation boost).
An alternative formulation of muscle activation dynamics was introduced by Hatze [6]: We divided the original equation from Hatze [6] by the parameter = 1.37 ⋅ 10 −4 mol/L which represents the maximum concentration of free Ca 2+ -ions in the muscle sarcoplasma. Thus, the values of the corresponding normalised Computational and Mathematical Methods in Medicine 3 concentration are 0 ≤ ≤ 1. The activity is finally calculated by the function and the parameter is shifted to the accordingly renormalised function with = ⋅ 0 and ℓ = 2.9. Two cases have been suggested by Hatze [13]: 0 = 6.62 ⋅ 10 4 L/mol (i.e., = 9.10) for ] = 2 and 0 = 5.27 ⋅ 10 4 L/mol (i.e., = 7.24) for ] = 3, which have been applied in the literature [7,8,14,15]. By substituting (2) and (3) intȯ= ( , ℓ CErel )/ ⋅̇and resubstituting the inverse of (3) afterwards, Hatze's formulation of an activation dynamics can be transformed into a nonlinear differential equation directly in terms of the activity: with the initial condition (0) = ,0 . The solutions ( ) and ( ) of both formulations of activation dynamics (1) and (5) can now be directly compared by integrating them with the same initial condition ,0 = ,0 using the same stimulation .

Local First and Second Order Sensitivity of ODE Systems regarding Their Parameters
Let Ω ⊆ R × R × R and : Ω → R . We then consider a system of ordinary, first order initial value problems (IVP): where ( ) = ( 1 ( ), 2 ( ), . . . , ( )) denotes the vector of state variables, = ( 1 , 2 , . . . , ) the vector of right hand sides of the ODE, and Λ = { 1 , 2 , . . . , } the set of parameters which the ODE depends on. The vector of initial conditions is abbreviated by The first order sensitivity of the solution ( , Λ) with respect to the parameter set Λ is defined as the matrix Simplifying, we denote = ( , Λ), = ( , , Λ), and = ( , Λ) but keep the dependencies in mind. Because the solution ( ) might only be gained numerically rather than in a closed-form expression, we have to apply the wellknown theory of sensitivity analysis as stated in Vukobratovic [16], Dickinson and Gelinas [17], Lehman and Stark [2], and ZivariPiran [18]. Differentiating (8) with respect to and applying the chain rule yield with / being the gradient of state variables. Hence we obtain the following ODE for the first order solution sensitivity: or in short termṡ where = ( ) is the × sensitivity matrix and = ( ) is the × Jacobian matrix with = ( / ) ; furthermore, = ( ) denotes the × -matrix containing the partial derivatives = ( / ) and 0 × denotes the ×matrix consisting of zeros only.
By analogy, the second order sensitivity of ( ) with respect to Λ is defined as the following × × -tensor: with assuming = for all = 1, . . . , , therefore assuming that the prerequisites of Schwarz theorem (symmetry of the second derivatives) are fulfilled throughout. Differentiating with respect to and applying the chain rule lead to the ODĖ with (0) = 0. For purposes beyond the aim of this paper, a condensed notation introducing the concept of tensor (or Kronecker) products as in ZivariPiran [18] may be helpful. For a practical implementation in MATLAB see Bader and Kolda [19].
Furthermore, if an initial condition ,0 (see (7)) is considered as another parameter, we can derive a separate 4 Computational and Mathematical Methods in Medicine sensitivity differential equation by rewriting (6) in its integral form Differentiating this equation with respect to 0 yields and differentiating again with respect to results in a homogeneous ODE for each component ,0 ( ); namely, The parameters of our analysed models are supposed to represent physiological processes and bear physical dimensions therefore. For example, and 1/ are frequencies measured in (Hz), whereas is measured in (mol/L). Accordingly, = ( / ) would be measured in (Hz) and in (s) (note that our model only consists of ODE and therefore we do not need a second index). Normalisation provides a comprehensive comparison between all sensitivities, even across models. For any parameter, the value fixed for a specific simulation is a natural choice. For any state variable, we chose its current value ( ) at each point in time of the corresponding ODE solution. Hence, we normalise each sensitivity = / by multiplying it with the ratio / ( ) to get the relative sensitivitỹ = ⋅ .
A relative sensitivitỹthus quantifies the percentage change in the th state variable value per percentage change in the th parameter value. This applies accordingly to the second order sensitivitỹ= It can be shown that this method is valid and mathematically equivalent to another common method in which the whole model is nondimensionalised a priori [20]. A nonnormalised model formulation has the additional advantage of usually allowing a more immediate appreciation of and transparent access for experimenters. In the remainder of this paper, we are always going to present and discuss relative sensitivity values normalised that way. In our model the specific case = 1 applies, so (10) and (14) simplify to the case = 1 (no summation).

Variance-Based Global Sensitivity Analysis
The differential sensitivity analysis above is called a local method because it does not take the physiological range of parameter values into account. Additionally factoring in such ranges characterises the so-called global methods. The main idea behind most global methods is to include a statistical component to scan the whole parameter space C and combine the percentage change in a state variable value per percentage change in a parameter value with the variability of all of the parameters. The parameter space C can be seen as a -dimensional cuboid C = [ − 1 ; + 1 ] × ⋅ ⋅ ⋅ × [ − ; + ], where − and + are the minimal and maximal parameter values and is the number of parameters. We can now fix a certain pointΛ = (̂1, . . . ,̂) ∈ C and calculate the local gradient of the solution with respect toΛ. The volume of the star-shaped area, investigated by changing only one parameter at once and lying within a ball around Λ, vanishes in comparison to C for an increasing number of parameters [21]. For an overview of the numerous methods like ANOVA, FAST, Regression, or Sobol's Indexing, the reader is referred to Saltelli and Chan [4] and Frey et al. [22].
In this section we want to sketch just the main idea of the variance-based sensitivity analysis approach as presented in Chan et al. [3], which is based on Sobol's Indexing. We chose this method because of its transparency and low computational cost. This method aims at calculating two measurands of sensitivity of a state variable with respect to parameter : the variance-based sensitivity function denoted by VBS ( ) and the total sensitivity index function denoted by TSI ( ). The VBS functions give a normalised first order sensitivity quite similar tõfrom the previous section but include the parameter range. The TSI functions, however, additionally include higher order sensitivities and give a measurand for interdependencies of parameter influences.
For this family perform the following calculations. Note that the computations in Chan et al. [3] are done using Monte-Carlo integrals as an approximation. The VBS and TSI can be finally calculated as The normalisation entails additional properties of VBS and TSI (see [3, Figure 1]): In other words, VBS ( ) gives the normalised global first order sensitivity function of the solution with respect to in relation to the model output range. Accordingly, TSI ( ) quantifies a relative impact of the variability in parameter on the model output, factoring in the interdependent influence in combination with all other parameters' sensitivities. Chan et al. [3] suggested to denote the TSI ( ) value as the "importance" of .

An Analytical Example for Local Sensitivity Analysis including a Link between Zajac's and Hatze's Formulations
By further simplifying Zajac's formulation of an activation dynamics (1) through assuming a deactivation boost = 1 (activation and deactivation time constants are equal) and a basic activity 0 = 0, we obtain a linear ODE for this specific case sp , which is equivalent to Hatze's equation (2) modelling the time evolution of the free Ca 2+ -ion concentration: By analysing this specific case, we aim at making the above described sensitivity analysis method more transparent for the reader. Solving (22) yields depending on just two parameters (stimulation: control parameter) and (time constant of activation: internal parameter) in addition to the initial value 0 = ,0 . The solution ( ) equals the value after about .
We apply the more generally applicable, implicit methods (10) and (17) to determine the derivatives of the solution with respect to the parameters (the sensitivities), although we already know solution (23) in a closed form. Hence, for the and insert these partial derivatives into (10) and (17). Solving the respective three ODEs for the three parameters ( , , and ,0 ) and normalising them according to (18) give the relative sensitivities of sp with respect to , , and ,0 as functions of time (see Figure 1): , A straightforward result is that the time constant has its maximum effect on the solution (Figure 1, seẽ( )) at time = . In case of a step in stimulation, the sensitivitỹ( ) vanishes in the initial situation and exponentially approaches zero again after a few further multiples of the typical period . Note that̃( ) is negative, which means that an increase in decelerates activation. Thus, for a fixed initial value ,0 , the solution value ( ) decreases at a given point in time if is increased. After a step in stimulation , the time in which the solution ( ) bears some memory of its initial value ,0 is equal to the period of being nonsensitive to any further step in (comparẽ, 0 ( ) tõ( ) and (25) to (27)). After about /2, the sensitivitỹ, 0 ( ) has already fallen to about 0.1 and ( ) to about 0.9 accordingly.

The Numerical Approach and Results
Typically, biological dynamics are represented by nonlinear ODEs. Therefore, the linear ODE used for describing activation dynamics in the Zajac [5] case (1) is more of an exception. For example, a closed-form solution can be given. Equation (23) is an example as shown in the previous section for the reduced case of nonboosted deactivation (22).
In general, however, nonlinear ODEs used in biomechanical modelling, as the Hatze [6] case (5) for describing activation dynamics, can only be solved numerically. It is understood that any explicit formulation of a model in terms of ODEs allows providing the partial derivatives of their right hand sides with respect to the model parameters in a closed form. Fortunately, this is exactly what is required as part of the sensitivity analysis approach presented in Section 3, in particular in (10).
As an application for applying this approach, we will now present a comparison of both formulations of activation dynamics. The example indicates that the approach may be of general value because it is common practice in biomechanical modelling to (i) formulate the ODEs in closed form and (ii) integrate the ODEs numerically. Adding further sensitivity ODEs for model parameters then becomes an inexpensive enhancement of the procedure used to solve the problem anyway.
For the two different activation dynamics [5,6], the parameter sets Λ and Λ , respectively, consist of including the initial conditions. The numerical solutions for these ODEs were computed within the MATLAB environment (The MathWorks, Natick, USA; version R2013b), using the preimplemented numerical solver 45 which is a Runge-Kutta algorithm of order 5 (for details see [23]).
The memory (influence on solution) of the initial value is lost after about 2 , almost independently of all other parameters. This loss in memory is obviously slower than in that case ,0 = 0 (initial value) and = 1 (for = 1 and 0 = 0 exactly the final value; see Section 5 and Figure 1). In that extreme case, the influence (relative sensitivity) of the lowest possible initial value ( ,0 = 0) on the most rapidly increasing solution (maximum possible final value: = 1) is lost earlier.

Relative Sensitivitỹ.
The influence of the time constant on the solution is reduced with decreasing difference between initial and final activity values (compare maximum values in Figures 1 and 2) and, no matter the value, with jointly raised levels of initial activity ,0 and , the latter determining the final activity value if = 1. When deactivation is slower than activation ( < 1: right column in Figure 2),̃is higher than in the case = 1, both in its maximum amplitude and for longer times after the step in stimulation, especially at low activity levels (upper rows in Figure 2).

Relative
Sensitivitỹ. Across all parameters, the solution in general is most sensitive to . However, the influence of the deactivation boost parameter is usually comparable. In some situations, this also applies to the activation time constant (see below). For = 1 (Figure 2, left), the solution becomes a little less sensitive to with decreasing activity level (̃< 1), which reflects that the final solution value is not determined by alone but by 0 > 0 and ̸ = 1 as much. If deactivation is much slower than activation ( = 1/3 < 1: Figure 2, right), we find the opposite to the = 1 case: the more the activity level rises, the lesser determines the solution. Additionally, stimulation somehow competes with both deactivation boost and time constant (see further below). Using the term "compete" illustrates the idea that any single parameter should have an individual interest in influencing the dynamics as much as possible in order not to be considered superfluous.

Relative Sensitivities̃,̃,̃.
At submaximal stimulation levels < 1, the final solution value is determined to almost the same degree by stimulation and deactivation boost , yet with opposite tendencies (̃> 0,̃< 0). As explained, both parameters compete for their impact on the final solution value. Only at maximum stimulation = 1 (lowest row in Figure 2), this parameter competition is resolved in favour of . In this specific case, does not influence the solution at all. For = 1 the competition about influencing the solution is intermittently but only slightly biased by : sensitivitỹpeaks at comparably low magnitude around = . This influence comes likewise intermittently at the cost of influence: the absolute value of̃rises a little slower thañ. In the case < 1, this competition becomes more differentiated and spread out in time. Again at submaximal stimulation and activity levels, the absolute value of̃is lower than that of̃but higher than that of̃, making all three parameters , , and compete to comparable degrees for an impact on the solution until about = 4 . Also,̃does not vanish before about = 10 .

Results for Hatze's Activation Dynamics: Sensitivity Functions.
We also simulated activation dynamics for the parameter set Λ (29), leaving now four of the values constant ( 0 = 0.005, = 10 1/s, ℓ = 2.9, ℓ CErel = 1) and again varying three others (initial condition ,0 , stimulation , and nonlinearity ]), keeping in mind that the eighth parameter ( ) is assumed to depend on ]. Time courses of the relative sensitivities̃( ) with respect to all parameters ∈ Λ are plotted (see Figure 3). In the left column of Figure 3, ] = 2, = 9.10 is used, in the right column ] = 3, = 7.24. Here, the same pairs of the parameter values ( 0 = 0.005 ≤ ,0 ≤ 0.5 and 0.01 ≤ ≤ 1, increasing from top to bottom; see legend of Figure 3) are used as in Section 6.1 (Figure 2).
Hatze's activation dynamics (5) are nonlinear unlike Zajac's activation dynamics (1). This nonlinearity manifests particularly in a changeful influence of the parameter ]. Additionally, the parameter is just roughly comparable to the inverse of the exponential time constant in Zajac's linear activation dynamics.
6.2.1. Relative Sensitivitỹ. In Zajac's linear differential equation (1), establishes a distinct time scale independent of all other parameters. The parameter in Hatze's activation dynamics (5) is just formally equivalent to the reciprocal of : the sensitivitỹdoes not peak stringently at = 1/ = 0.1 s but rather diffusely between about 0.05 s and 0.1 s in both of the cases ] = 2 and ] = 3. At first this is not surprising because the scaling factor in Hatze's dynamics is ] ⋅ rather than just . However, ]⋅ does not fix an invariant time scale for Hatze's nonlinear differential equation. This fact becomes particularly prominent at extremely low activity levels for ] = 2 (Figure 3, left, top row) and up to moderately submaximal activity levels for ] = 3 (Figure 3, right, top two rows). Here, is negative, which means that increasing the parameter results in less steeply increasing activity. This observation is counterintuitive to identifying with a reciprocal of a time constant like . Rather than being expected from the product ]⋅ , the exponent ] does not linearly scale the time behaviour becausẽpeaks do not occur systematically earlier in the ] = 3 case as compared to ] = 2.

Relative Sensitivitỹ,
0 . Losing the memory of the initial condition confirms the analysis of time behaviour based oñ. At high activity levels ( Figure 3, bottom row), Hatze's activation dynamics loses memory at identical time horizons (no matter the ] value) seemingly slower for higher ] at intermediate levels ( Figure 3, two middle rows) and clearly faster at very low levels ( Figure 3, top row). The parameter still does roughly determine the time horizon in which the memory of the initial condition ,0 is lost and the influence of all other parameters is continuously switched on from zero influence at = 0.

Relative Sensitivitỹ0 .
As in Zajac's dynamics the solution is generally only sensitive to 0 at very low stimulation levels ≈ 0 (Figure 3, top row). At such levels, the ] = 3 case shows the peculiarity that the solution becomes strikingly insensitive to any other parameter than 0 itself (and ,0 ). The time evolution of the solution is more or less determined by just this minimum ( 0 ) and initial ( ,0 ) activities, and determining the approximate switching time horizon between both. The ℓ CE dependency, constituting a crucial property of Hatze's activation dynamics, is practically suppressed for ] = 3 at very low activities and stimulations. In contrast,̃ℓ CErel remains for ] = 2 on a low but still significant level of about a fourth of the three dominating quantities̃0, ,0 , and̃].

Relative Sensitivitỹ]. The sensitivity with respect to
] is extraordinarily high at low activities and stimulations around 0.1, both for ] = 2 and for ] = 3 (Figure 3, second row from top), additionally at extremely low levels for ] = 2 ( Figure 3, left, top row). At moderately submaximal levels ( Figure 3, third row from top), the solution is influenced with an already inverted tendency (̃] changes sign to positive) after around a 1/ time horizon for ] = 2. However, at these levels the solution is practically insensitive to ] for any ]. At high levels ( Figure 3, bottom row) we find that there is no change in the character of time evolution of the solution, despite the specific value of ]. The degree of nonlinearity ] is unimportant because the time evolution and the ranking of all other sensitivities are hardly influenced by ]. In both cases, the rise in activity is quickened by increasing ] (̃] > 0), as opposed to low activity and stimulation levels where rises in activity are slowed down (̃] < 0; see also above).

Relative Sensitivities̃,̃,̃ℓ
CErel , and̃ℓ . Of all the remaining parameters, stimulation , scaled maximum free Ca 2+ -ion concentration , relative CE length ℓ CErel , and the pole ℓ of the length dependency in Hatze's activation dynamics, the latter has the lowest influence on the solution. The influence characters of all four parameters are yet completely identical. Their sensitivities are always positive and coupled by fixed scaling ratios due to all of them occurring within just one product on the right side of (5).̃and̃are identical, while the sensitivity with respect to ℓ CErel is the highest, with ratios̃ℓ CErel /̃ℓ ≈ 3 and̃ℓ CErel /̃≈ 1.2. Except at very low activity (where 0 plays a dominating role) and except for the generally changeful ] influence, these are the four parameters that dominate the solution after an initial phase in which the initial activity ,0 determines its evolution. The parameter does not have a strong direct influence on the solution. As stated above, it defines the approximate time horizon in which the ,0 influence gets lost and all other parameters' influence is switched on from zero at = 0. Table 1 pools the lower and upper boundaries for every parameter in Λ and Λ used in our calculations. We refer to Hatze [24], Zajac [5], or Günther et al. [11] for traceability of our choices. The left hand side of Figure 4 shows the VBS functions of every parameter in Λ of Zajac's model. The plotted functions can be compared to our previously computed relative first order sensitivity functions from Figure 2: at first sight,̃, 0 and VBS ,0 look equal, but the VBS function indicates a slightly increased duration of influence of ,0 . Regarding , the VBS function peaks at the same time as̃, but with a smaller amplitude. Likewise, the courses of VBS and VBS are comparable tõand̃from the second and third row of Figure 2. The calculated VBS functions in the Zajac case show what would be expected intuitively: a VBS represents a parameter's mean influence averaged over its range of values. Additionally, we plotted the sum of all first order sensitivities. This sum indicates which amount of the total variance is covered by first order sensitivities. The closer the sum to 1 the smaller the impact of the second and higher order sensitivities.

Variance-Based Sensitivity (VBS) and Total Sensitivity Indices (TSI) for Zajac's and Hatze's Activation Dynamics.
The right hand side of Figure 4 shows the TSI functions of every parameter in Λ of Zajac's model. Generally, there are only minor deviations of the TSI functions from their counterparts VBS . That is, the influence of none of the parameters is significantly enhanced by an interdependent effect in combination with other parameters. According to both analyses, there are just four globally important parameters that govern the system's state throughout the whole examined solution space: the initial condition ,0 within a typical time horizon after a step in , the new stimulation level determining activity after about , the deactivation boost with smaller impact than , and determining the time horizon itself.
The left hand side of Figure 5 shows the VBS functions of every parameter in Λ of Hatze's model. Very similar to the Zajac case, the calculated VBS seemingly represent to a high degree a parameter's mean influence averaged over its range of values (compare Figure 3). As in the Zajac case, there are four globally important parameters, according to both VBS and TSI analyses. Compared to Zajac's model, the interdependent effect in combination with other parameters (TSI: right hand side of Figure 5) is more pronounced for two parameters: both the stimulation and the CE length ℓ CE importance are distinctly higher than their first order effects as expressed by VBS functions. Furthermore, the time horizon within the initial condition ,0 has an aftereffect in response to a step in globally a little higher in VBS as compared to local sensitivity analysis (Figure 3). In addition, the time horizon of ,0 is clearly enhanced by interdependencies with other parameters (TSI: right hand side of Figure 5).
Altogether, VBS versus TSI analysis substantiate local first and second order sensitivity analyses: for one thing, Hatze's model is more inert against steps in stimulation than Zajac's model. For another thing, the dynamics described by Hatze's model incorporates stronger nonlinear coupling effects from combinations of parameters than Zajac's model. These latter effects are better seen in detail when looking at local sensitivities, that is, analysing just small and selected volumes of the parameter space C. In turn, VBS and TSI provide a broad but coarse overview about first and higher order sensitivities of all parameters.

A Bottom Line for Comparing Zajac's and Hatze's Activation Dynamics: Second Order
Sensitivities. At first sight, Zajac's activation dynamics [5] is more transparent because it is descriptive in a sense that it captures the physiological behaviour of activity rise and fall in an apparently simple way. It thereto utilises a linear differential equation with wellknown properties, allowing for a closed-form solution. It needs only four parameters to describe the Ca 2+ -ion influx to the muscle as a response to electrical stimulation: the stimulation itself as a control parameter, the time constant for an exponential response to a step increase in stimulation, a third parameter (deactivation boost) biasing both the rise and fall times, and the saturation value | ∞ of activity which in turn depends on and the basic activity 0 being the fourth parameter. The smaller the < 1 is (deactivation slower than activation), the faster the very activity level | ∞ =1 = 0 + ⋅ (1 − 0 ) is reached, at which saturation would occur for = 1. Saturation for < 1 occurs at a level | ∞ = controller to deal with a muscle as the biological actuator.
Regarding the nonlinearity exponent ], solution sensitivity further depends nonmonotonously on activity level, partly even with the strongest influence, partly without any influence. We also found that the solution is more sensitive to its parameters , ℓ CErel , and ℓ than is Zajac's activation dynamics to any of its parameters.
This higher complexity of Hatze's dynamics becomes even more evident by analysing the second order sensitivities (see (13) or (19) for their relative values). They express how a first order sensitivity changes upon variation of any other model parameter. In other words, they are a measure of model entanglement and complexity. Here, we found that the highest values amongst all relative second order sensitivities Computational and Mathematical Methods in Medicine 13 in Zajac's activation dynamics are about −0.8 (̃) and 1.6 (̃). In Hatze's activation dynamics, the highest relative second order sensitivities are those with respect to ] or ℓ CErel (in particular for , , and ], ℓ CErel themselves) with maximum values between about −8.0 (̃ℓ CErel ] ,̃] ) and 13.4 (̃ℓ CErel ℓ CErel ,̃ℓ CErel ,̃ℓ CErel ,̃] ] at submaximal activity). That is, they are an order of magnitude higher than in Zajac's activation dynamics.
Yet, we have to acknowledge that Hatze's activation dynamics contains crucial physiological features that go beyond Zajac's description.

A Plus for Hatze's Approach: Length
Dependency. It has been established that the length dependency of activation dynamics is both physiological [7] and functionally vital [15] because it largely contributes to low-frequency muscle stiffness. It has also been verified that Hatze's model approach provides a good approximation for experimental data [7]. In that study, ] = 3 was used without comparing to the ] = 2 case. There seem to be arguments in favour of ] = 2 from a mathematical point of view. In particular, the less changeful scaling of the activation dynamics' characteristics down to very low activity and stimulation levels, at which some CE length sensitivity remains, seem to be an advantage when compared to the ] = 3 case. Up to this point, we have argued solely mathematically. It is, however, physiological reality that is eventually aimed at. We therefore repeated the model fit done by Kistemaker et al. [7] while now allowing a variation in ] and in force-length relations.

An Optimal Parameter Set for Hatze's Activation Dynamics
Plus CE Force-Length Relation. Sensitivity analysis allows rating Hatze's approach as an entangled construct. Additionally, Kistemaker et al. [7] decided to choose ] = 3 without giving a reason for discarding ] = 2. It further seemed that they did not perform an algorithmic optimisation across various submaximal stimulation levels to find a muscle parameter set, which best fits known shifts Δℓ CE,opt,submax = ℓ CE,opt − ℓ CE,opt,submax in optimal, submaximal CE length ℓ CE,opt,submax at which isometric force isom = isom ( , ℓ CE ) peaks. Accordingly, it seemed worth performing such an optimisation because isom generally depends on length ℓ CE and activity , and the latter may be additionally biased by an ℓ CE -dependent capability for building up cross-bridges at a given level of free Ca 2+ -ions in the sarcoplasma, as formulated in Hatze's approach: isom ( , ℓ CE ) = max ⋅ ( , ℓ CE ) ⋅ ℓ (ℓ CE ). Thus, a shift in optimal CE length Δℓ CE,opt,submax with changing can occur depending on the specific choices of both the length dependency of activation ( , ℓ CE ) (see (3) and (4)) and the CE's force-length relation ℓ (ℓ CE ).
Consequently, we searched for optimal parameter sets of Hatze's activation dynamics in combination with two different force-length relations ℓ (ℓ CE ): either a parabola [7] or bell-shaped curves [11,25]. For a given optimal CE length ℓ CE,opt = 14.8 mm [26] representing a rat gastrocnemius muscle and three fixed exponent values ] = 2,3,4 in Hatze's activation dynamics (all other parameters as given in Section 2), we thus determined Hatze's constant 0 and the width parameters of the two different force-length relations ℓ (ℓ CE ) (WIDTH in Kistemaker et al. [7] and van Soest and Bobbert [9] and Δ asc = Δ des = Δ in Mörl et al. [25], resp.) by an optimisation approach. The objective function to be minimised was the sum of squared differences between the Δℓ CE,opt,submax values as predicted by the model and as derived from experiments (see Table 2 in Kistemaker et al. [7]) over five stimulation levels = 0.55, 0.28, 0.22, 0.17, 0.08. Note that = applies in the isometric situation (see (2) and compare (3)). Further note that experimental data for muscle contractions at very low stimulation levels are missing in the literature so far: the lowest analysed level available for Kistemaker et al. [7] was = 0.08, that is, comparable to the second rows from top in Figures 2 and 3.
The optimisation results are summarised in Table 2. The higher the ] value, the smaller the optimisation error. The predicted width values WIDTH or Δ , respectively, decrease along with the error. We would yet tend to exclude the case ] = 4 because the predicted width values seem unrealistically low when compared to published values from other sources (e.g., WIDTH = 0.56 [9], Δ = 0.35 [25]). Furthermore, 0 decreases with ] using the parabola model for ℓ (ℓ CE ) whereas it saturates between ] = 3 and ] = 4 for the bell-shaped model. The bell-shaped model shows the most realistic Δ in the case ] = 3 (Δ = 0.32). Fitting the same model to other contraction modes of the muscle [25], a value of Δ = 0.32 had been found. In contrast, when using the parabola model, realistic WIDTH values between 0.5 and 0.6 are predicted by our optimisation for ] = 2.
When comparing the optimised parameter values across all start values of the ℓ (ℓ CE ) widths, across all ] values, and across both ℓ (ℓ CE ) model functions, we find that the resulting optimal parameter sets are more consistent for bellshaped ℓ (ℓ CE ) than for the parabola function. The bellshaped force-length relation gives generally a better fit. For each single ] value, the corresponding optimisation error is smaller when comparing realistic, published WIDTH and Δ values that may correspond to each other (WIDTH = 0.56 [9] and Δ = 0.35 [25]). Additionally, the error values from our optimisation are generally smaller than the corresponding value calculated from Table 2 in Kistemaker et al. [7] (0.23 mm).
In a nutshell, we would say that the most realistic model for the isometric force isom at submaximal activity levels is the combination of Hatze's approach for activation dynamics with ] = 3 and a bell-shaped curve for the force-length relation ℓ (ℓ CE ) with ] asc = 3. As a side effect, we predict that the parameter value 0 , being a weighting factor of the first addend in the compact formulation of Hatze's activation dynamics (5), should be reduced by about 40% ( 0 = 3.25 ⋅ 10 4 L/mol) as compared to the value originally published in Hatze [13] ( 0 = 5.27 ⋅ 10 4 L/mol).

A Generalised Method for Calculating
Parameter Sensitivities. The findings in the last section were initiated by thoroughly comparing two different biomechanical models of muscular activation using a systematic sensitivity analysis as introduced in Dickinson and Gelinas [17] and Lehman and Table 2: Parameters minimising the sum over five submaximal stimulation levels = = 0.55, 0.28, 0.22, 0.17, 0.08 of squared differences between shifts in optimal CE length Δℓ CE,opt,submax ( ) (Δ MA,opt by Roszek et al. (1994) [27] in the third column of Table 2 in Kistemaker et al. [7]) at these levels predicted by the model with the isometric force isom ( , ℓ CE ) = max ⋅ ( = , ℓ CE ) ⋅ ℓ (ℓ CE ) and by experiments; simulated data represent a rat gastrocnemius muscle with an optimal CE length ℓ CE,opt = 14.8 mm [26]; start value of 0 was 6.0⋅10 4 L/mol; the exponents of the bell-shaped force-length relations ℓ (ℓ CE ) were fixed according to Mörl et al. [25] (] asc = 3, ] des = 1.5); the corresponding width values in the ascending and descending branch were assumed to be equal: Δ asc = Δ des = Δ ; van Soest and Bobbert [9] and Kistemaker et al. [7] used a parabola for ℓ (ℓ CE ); for all other model parameters see Sections 7.3 and 2; optimisation was done by minsearch (Nelder-Mead algorithm) in MATLAB with error tolerances of 10 −8 ; error is the square root of the above-mentioned sum divided by five; corresponding error value given in Table 2 in Kistemaker et al. [7] was 0.23 mm.

]
Bell-shaped [11,25] P a r a b o l a [ 7,9]  Stark [2], respectively. Starting with the latter formulation, Scovil and Ronsky [1] calculated specific parameter sensitivities for muscular contractions. They applied three variants of this method. Method 1 applies to state variables that are explicitly known to the modeller as in, for example, an eye model [2], a musculoskeletal model for running that includes a Hilltype muscle model [1], or the activation models analysed in our study. Scovil and Ronsky [1] calculated the change in the value of a state variable averaged over time per a finite change in a parameter value, both normalised to each of their unperturbed values. They thus calculated just one (mean) sensitivity value for a finite time interval (e.g., a running cycle) rather than time-continuous sensitivity functions.
Method 2: whereas Dickinson and Gelinas [17] and Lehman and Stark [2] had introduced the full approach for calculating such sensitivity functions, Scovil and Ronsky [1] distorted this approach by suggesting that the partial derivative of the right hand side of an ODE, that is, of the rate of change of a state variable, with respect to a model parameter would be a "model sensitivity. " The distortion becomes explicitly obvious from our formulation: this partial derivative is just one of the two addends that contribute to the rate of change of the sensitivity function (10), rather than defining the sensitivity of the state variable itself (i.e., the solution of the ODE) with respect to a model parameter (8).
Method 3: Scovil and Ronsky [1] had also asked for calculating the influence of, for example, a parameter of the activation dynamics (like the time constant) on an arbitrary joint angle, that is, a variable that quantifies the overall output of a coupled dynamical system. Of course, the time constant does not explicitly appear in the mechanical differential equation for the acceleration of this very joint angle, which renders applicability of method 2 impossible. The conclusion in Scovil and Ronsky [1] was to apply method 1. Here, the potential of our formulation comes particularly to the fore. It enables calculating the time-continuous sensitivity of all components of the coupled solution, that is, any state variable ( ). This is because all effects of a parameter change are in principle reflected within any single state variable, and the time evolution of a sensitivity according to (10) takes this into account.
In this paper, we have further worked out the sensitivity function approach by Lehman and Stark [2], presenting the differential equations for sensitivity functions in more detail to those modellers who want to apply the method. Furthermore, we enhanced the approach by Lehman and Stark [2] to also calculating the sensitivities of the state variables with respect to their initial conditions (17). This should be helpful not only in biomechanics but also, for example, in meteorology when predicting the behaviour of storms [28]. Since initial conditions are often just known approximately but start with the relative sensitivity values of 1, their influence should be traced to verify how their uncertainty propagates during a simulation. In the case of muscle activation dynamics, the sensitivities̃, 0 and̃, 0 , respectively, decreased rapidly to zero: initial activity has no effect on the solution very early before steady state is reached.
Furthermore, we included a second order sensitivity analysis which is not only helpful for an enhanced understanding of the parameter influence but also part of mathematical optimisation techniques [29]. The values of̃could be interpreted either as the relative sensitivity of the sensitivitỹ with respect to another parameter (and vice versa: with respect to ) or as the curvature of the graph of the solution ( ) in the + -dimensional solutionparameter space. The latter may help to connect the results to the field of mathematical optimisation in which the second derivative (Hessian) of a function is often included in objective functions to find optimal parameter sets.

Insights into Global
Methods. Some additional conclusions can be drawn from global sensitivity analysis, in particular from comparing results in Section 6.3 to those based on local sensitivity analysis (Sections 6.1, 6.2, and 7.1).
For Zajac's activation dynamics, global analysis confirms local analysis in stating that there are no significant second or higher order sensitivities, with the slight exception of the phase of rapid change in activity after a step in stimulation. An experimenter who wants to measure the activation time constant can exclude influence from potentially slower deactivation processes ( < 1) by starting from high activity levels (Figure 2, bottom). It should yet be kept in mind that build-up of activity to the new level is not solely determined by but might be biased by other parameters than because TSI peaks during the build-up phase (Figure 4, right).
In Hatze's activation dynamics, the higher order sensitivities play a clearly more significant role, even in the nearsteady-state case ( Figure 5: stronger deviation from 1 of both VBS and TSI). When arguing in terms of controllability of the models in Section 7.1, we speculated that Zajac's dynamics might be easier to control than Hatze's dynamics. Notwithstanding, Figure 5 shows that the stimulation is the most important control factor with even a higher importance than in Zajac's formulation.
At first sight unapparent, another result is the importance of . From a strictly local point of view we concluded that this parameter should have the same sensitivity as since they both are formally equivalent multipliers in Hatze's ODE (see relative sensitivities in Figure 3). However, the importance of is significantly smaller than that of , in fact almost negligible. Their different global variabilities of values can give an explanation. The parameter in the product ⋅ ∈ [4; 11] × [0; 1] has a clearly lower relative variability than , measured in maximum percentage deviation from the respective mean value. The parameter thus acts as an amplifier for . Similarly, the parameter ] has a relatively small variability throughout the literature. So, although its differential sensitivity is quite large, ] is found to have a low importance for the model output. For the latter fact there is yet another reason. In Section 6.2, we have emphasised that ] has a very changeful influence on solutions, depending on activity level. Additionally, its influence is highly dependent on other parameters like length ℓ CE and (see end of Section 7.1). Its strong influence in some situations and configurations is thus hidden by global averaging.
This demonstrates that the findings of global sensitivity analysis must be treated with caution because the whole dynamics of a system is condensed to a single average function per whole parameter range. Without local analyses of the solution space as exemplified in Sections 6.1 and 6.2 crucial features of its topology might be lost when solely relying on global analysis. Activation frequency constant in Hatze [6]; value: range: 3.67, . . . , 11.25 (1/s); here: 10 (1/s) :