Sensitivity Analysis of Upper Limb Musculoskeletal Models During Isometric and Isokinetic Tasks

Personalized parameters can be used to quantitatively describe the musculoskeletal system. These parameters (e.g., height, weight, physiological cross-sectional area) are typically derived from subject-specific anatomy, which is known to vary due to size, age, sex, race, and lifestyle [1–5]. Researchers use personalized parameters to generate subject-specific musculoskeletal models and evaluate human mechanics through computer simulations. Subject-specific musculoskeletal models vary in detail from simply scaling a model by height and weight to deriving personalized parameters from three-dimensional anatomical reconstructions of medical images [6–9]. No matter which methods are used to develop subject-specific models, the personalized parameters do not exactly match the subject due to the inherent errors of measurement equipment. Understanding how measurement or parameter uncertainty influences musculoskeletal model predictions over a diverse range of healthy and pathologic subjects is needed for subject-specific models to be more clinically and scientifically relevant. For example, current methods for subject-specific modeling require extensive imaging and data processing, but can still result in high errors (e.g., [10,11]). However, parameter sensitivity can address these issues by identifying key parameters to measure, thereby narrowing the scope and experimental time needed for model for personalization. Sensitivity analyses can also provide error bounds on predictions, thereby improving confidence in and interpretability of results.

Toward understanding parameter uncertainty, sensitivity analyses have been used to evaluate how changes in musculoskeletal parameters influence lower limb simulations. For example, Scovil and Ronsky simulated isolated muscle contraction and locomotion, while individually altering 14 Hill-type model parameters by ±50% across all muscles; they concluded that simulations of movement (walking and running) were less sensitive to parameter perturbations than simulations of muscle contraction [12]. Scovil and Ronsky also concluded that each simulation was most sensitive to a different set of parameters [12]. Ackland et al. expanded this work by using a Monte Carlo approach to explore sensitivity to multiple simultaneous parameter changes during gait, including muscle moment arms, optimal fiber length, maximum isometric force, and tendon slack length [13]. Ackland et al. concluded that predicted muscle forces during gait were most sensitive to tendon slack length, but the magnitude of sensitivity varied between individual muscles [13]. De-Groote et al. similarly demonstrated a high sensitivity to tendon slack length for certain muscles in a study that applied sensitivity analyses to estimate subject-specific model parameters using dynamometry with isometric knee flexion/extension [14]. Interestingly, De-Groote et al. also noted that some muscles had considerable sensitivity to optimal fiber length, but sensitivity across all parameters varied based on the method used to calculate muscle forces [14].

Prior studies provide key insights into the relative importance of various musculoskeletal parameters. However, to our knowledge, prior sensitivity analyses have not explored the wide range of physiological diversity that is present in humans. Notably, in prior studies that calculate sensitivity coefficients by computing the slope between two points, the points have been chosen by changing muscle parameters by a range of ±1% to ±10% [11,13–17]. Scovil et al. is the only outlier at ±50% [12]. These small ranges are not representative of reported muscle and bone changes in humans. For example, maximum isometric force is known to decrease by approximately 30% due to age-related atrophy [2,4] and can increase by approximately 20% due to exercise-induced hypertrophy [3,5]. Furthermore, by defining sensitivity using a two-point slope, the computed sensitivity is only valid if the relationship between muscle parameters and simulated outputs is linear. However, the underlying curves of the Hill-muscle model are nonlinear, and do not always scale linearly across large parameter changes [18]. Previous studies also predominately focus on locomotion [11–16]. Carmichael and Liu started to address this in a sensitivity study of atrophy that evaluated maximum forces in 46 isometric arm positions [17]. Atrophy was simulated by reducing maximum isometric force by 25%, 50%, and 62.5%, and then changing four other muscle parameters by ±1%; their results indicated overall parameter sensitivity increased with simulated atrophy [17]. However, the isometric tasks did not capture the wide range of upper limb movement, which has over 15 degrees-of-freedom [19], 30 forms of grasping [20], and can perform identical tasks with power or precision [21].

In this context, the objectives of this study were twofold: (1) compare different methods (two-point sensitivity coefficient, linear regression sensitivity coefficient, and sensitivity function) for calculating parameter sensitivity across a large physiological range of musculoskeletal parameters in the upper limb and (2) analyze the sensitivity of upper limb models to maximum isometric force (⁠$Fmmax$⁠), optimal fiber length (⁠$Lmopt$⁠), pennation angle (α), and mass during isometric and isokinetic tasks at the hand, wrist, and elbow. We explicitly evaluated parameter ranges inclusive of both healthy and pathological populations during three isometric and two isokinetic tasks. We hypothesized that sensitivity functions would most accurately account for the variation in simulated muscle activations caused by measurement uncertainty, as it would take into account the local nonlinear region of the Hill-type parameters.

Sensitivity analysis was performed using 2005 computed muscle control simulations performed in opensim 4.1 [22]. The muscle and bone parameters of two generic upper limb models were varied stepwise using reported values (Table 1). Computed muscle control simulations were then batch processed using the opensim API in python 3.7 to predict the muscle activations required to perform five tasks. The muscle activations from each simulation were summarized using the area under the curve (AUC). Sensitivity coefficients were calculated and statistically compared for the AUCs of each parameter using three methods: two-point method, linear regression method, and sensitivity function.

We performed sensitivity analysis individually for four musculoskeletal parameters: $Fmmax$⁠, $Lmopt$⁠, α, and mass. While previous studies included tendon slack length, we chose not to study this parameter as there is no way to measure tendon slack length in vivo, and many models use it as a tuning variable.

Each parameter was altered stepwise over reported ranges for healthy and pathologic subjects. The generic opensim models [24–26], one validated for lateral pinch and another for elbow and wrist movements, were altered by applying a uniform percent change to the selected parameters across all modeled muscles and bones. The ranges for the percent changes varied for each parameter based on reported values for hypertrophy, atrophy, osteoporosis, and osteopetrosis (Table 1). While it is not known how parameters directly map to subject diversity due to age, sex, lifestyle, etc., we assumed that by using bounds associated with disease and/or physical training, we would cover a large range of this diversity. The generic models were changed 100 times for each parameter, with 50 steps each for the positive and negative directions. This approach created 401 parameter sets (100 changes × 4 parameters + baseline), referred to as synthetic subjects. Note, this approach was inclusive of 50 sets for each parameter that represented atrophy/osteoporosis (decreases in parameters) and 50 sets for each parameter that represented hypertrophy/osteopetrosis (increases in parameters).

Isometric and isokinetic tasks were tested to capture multiple upper limb functions. Isometric tasks are used to grasp and support objects, while isokinetic tasks are used to move the limb and position objects. Five tasks (two isokinetic and three isometric) were selected to test a portion of upper limb tasks resulting in 2005 simulations (401 synthetic subjects × 5 tasks).

All isokinetic tasks started with the models positioned with a neutral shoulder, a fully pronated forearm (90 deg), and a flexed elbow (90 deg) (Fig. 1). This initial posture was chosen to match experimental work [27–30], while also ensuring a consistent initial state across simulations. The isokinetic tasks were designed to simulate knocking motions to capture both flexion and extension. To perform each motion, either the elbow or wrist joint was rotated in the sagittal plane to an angle that was halfway between the initial posture and maximum angle. The decision to simulate motion to the halfway point was determined based on the time required to run a single simulation and initial pilot work showing changes in activation occurred primarily at the motion boundaries, regardless of where these boundaries were set. Upon reaching the target angle, the joint was rotated back to the initial posture. Isokinetic elbow flexion was simulated by flexing the elbow from 90 deg to 110 deg for 0.27 s, stopping for 0.01 s, and then returning to 90 deg in 0.27 s for a total time of 0.54 s. Likewise, isokinetic wrist extension was simulated by extending the wrist from neutral to 48.5 deg for 0.2 s, stopping for 0.01 s, and then returning to neutral in 0.2 s for a total time of 0.4 s. The kinematics were defined to mirror experimental tasks performed using isokinetic dynamometers, which include brief periods of rest, acceleration, and deceleration at the endpoints [31,32]. The flexion/extension rates correspond to approximately 80 deg/s (elbow) and 194 deg/s (wrist), which is within the range used to experimentally evaluate isokinetic flexion/extension [31,32].

Three isometric tasks were simulated: isometric elbow flexion, isometric wrist extension, and lateral pinch (Fig. 1). Elbow and wrist tasks were simulated, as these joints are required for stably positioning the hand [32]. Shoulder tasks were excluded due to the increased complexity of modeling the scapulothoracic joint. Lateral pinch was selected as the hand task because it is a commonly used hand grip [33,34].

The isometric tasks used the same initial posture as the isokinetic tasks. However, as the tasks were isometric, the kinematics were held constant. An external force profile was applied to the model to simulate isometric force against a known load. The force profiles consisted of a linear ramp to the target force, followed by a plateau in which the force was held. The force ramp is included to provide an initialization time period needed for computational stability and to limit simulation failure. The target force equaled 50% of the force required for the baseline model predict a single muscle activation of 1 throughout the task; this selection prevented ceiling effects in individual muscle activations during atrophy that would impede calculation of the sensitivity coefficient.

Isometric elbow flexion was simulated by applying the force profile to the distal radius; the force profile included a 1 s ramp to a target force of 200 N followed by holding the force for 0.5 s. Isometric wrist extension was simulated by applying the force profile to the third metacarpal; the force profile included a 1 s ramp followed by holding 100 N for 0.5 s. These tasks used the same musculoskeletal model as the isokinetic tasks. To ensure physiologically reasonable solutions and mimic experimental studies that constrained subjects to a given posture using ergometers, all modeled joint were constrained, except for either the elbow or wrist [35,36]. This constraint was implemented by locking the joint, which effectively sets the joint angle to a given value and prevents deviation from that value [37,38].

A previously validated model was used to simulate lateral pinch, as the low masses and inertias at the hand require unique considerations (compared to the wrist or elbow) to accurately simulate force producing tasks at the hand [24,26,39,40]. The model was positioned as described in previous works with a neutral wrist, thumb positioned in a pinch posture, and the force profile applied to the thumb-tip. The force profile consisted of a 2 s linear ramp from 10 N to 40 N, followed by sustaining 40 N for 1 s.

Sensitivity coefficients are calculated by identifying how known changes in parameters influence predicted metrics. Herein, each sensitivity coefficient describes the change in parameter versus the change in predicted muscle activations.

To calculate sensitivity coefficients, a single value output metric is required regardless of the method being tested. Similar to prior studies [12,15], to reduce the output time-series muscle activations and force profiles (Fig. 2), the AUC was calculated. Choosing the AUC as the summary metric allows for sensitivity coefficients from tasks of different durations and/or magnitudes to be compared. Herein, a low-pass 4th order Butterworth filter with a cutoff frequency of 6 Hz was used to filter the predicted muscle activations [7] (Fig. 2). The AUC was calculated using NumPy's trapezoidal method in python 3.7 for the entire simulation time for isokinetic tasks, but only the simulation time after the force ramp for isometric tasks (i.e., time after initialization during which isometric task is performed). The AUCs were task normalized using the theoretical maximum and minimum AUC. The theoretical maximum and minimum respectively represent a muscle maintaining activation of 1 (fully active) or 0 (no activity) during the simulation.

In sensitivity analysis, the sensitivity coefficient (ε) is defined as the ratio of change in model behavior to changes in parameter behavior [41]. The coefficient is calculated by dividing the difference of the normalized metric (⁠$ΔM/M0$⁠) by the difference of the normalized parameter value (⁠$ΔP/P0$⁠), where $P0$ and $M0$ are the original values. The sensitivity coefficient can then be used to determine how uncertainty in parameter measurements influences simulated predictions. This is done by multiplying ε by the parameter percent change to calculate the impact of a parameter on a predicted metric. In this paper, ε of 0.01 means that if the parameter were changed by 1%, the normalized AUC would change by 0.01. Since the AUC is normalized, a change of 0.01 corresponds to a 1% change in AUC.

Three different methods for calculating ε were evaluated (Fig. 3): two-point sensitivity coefficient, linear regression sensitivity coefficient, and sensitivity function. For the two-point sensitivity coefficient (⁠$εtp$⁠), the original parameter (⁠$P0$⁠) is changed by a single positive (⁠$P+$⁠) and negative (⁠$P−$⁠) parameter change. Herein, we used parameter changes of ±1% to calculate the $εtp$⁠, this is similar to previous works [17]. The linear regression sensitivity coefficient (⁠$εlrs$⁠) computes ε as the linear regression slope. This approach has been used in studies that evaluate sensitivity across multiple time points [12,41]. These studies develop ε by averaging $εtp$ across all points, and since ε represents the slope of the line between the two points, the average sensitivity is equivalent to the linear regression slope [14]. Finally, a sensitivity function can be created using a series of parameter adjustments, and a sensitivity coefficient (⁠$εder$⁠) can be computed as the local numerical derivative for each point [41]. We computed the numerical derivative using a sliding window. The window was set so that the synthetic subjects for the positive and negative parameter were five steps away from the centered subject. The step size was determined based on reported errors in calculating subject-specific parameters from magnetic resonance images [42]. Note, instead of producing a single sensitivity coefficient, this method created a sensitivity function that was equivalent to the numerical derivative of the AUC versus parameter change. $εtp$ and $εlrs$ assumed a uniform parameter sensitivity across possible parameter changes, while $εder$ calculated a localized form of $εtp$ at individual parameter values. For linear functions and small sampling regions $εtp$⁠,$εlrs$⁠, and $εder$ would all approximate each other.

To determine which method for calculating ε best accounts for measurement uncertainty across the population, the accuracy of predicted AUCs was compared to simulated ground truths. A sliding window was used to predict the AUC of a synthetic subject five steps away in both the positive and negative direction. This represents using ε to determine how a 5% measurement error affected the predicted muscle activations. The 5% measurement error was selected based on Holzbauer et al. who reported a maximum variation in reconstructed muscle volumes of 4.4% [42]; however, this number was rounded to 5% as all ranges were divisible by 5. The sliding window also allowed the simulated AUC of the synthetic subjects five steps away to serve as the ground truth for each predicted AUC. The root mean squared error (RMSE) of each sensitivity method was then calculated for each set of predicted AUCs. Statistical differences between the three sensitivity methods were evaluated using analysis of variance for each muscle within each task. A Bonferroni posthoc was used for muscles that were found to have significant differences (p < 0.05).

Sensitivity coefficients calculated using the sensitivity function (⁠$εder$⁠) had lower RMSEs compared to coefficients calculated from linear regression slopes (⁠$εlrs$⁠) or two-point approximations (⁠$εtp$⁠). $εder$ generally had lower RMSE regardless of muscle or task (Fig. 4) and had significant differences for muscles in the isometric wrist extension, isometric elbow flexion, wrist knock, and elbow knock tasks when compared with AUCs predicted with $εlrs$ and $εtp$⁠. For isometric wrist extension, muscle activation AUCs computed with $εder$ were found to be significantly lower than $εtp$ for extensor carpi radialis longus (ECRL) (p < 0.01), extensor carpi radialis brevis (ECRB) (p < 0.05), and extensor pollicis longus (EPL) (p < 0.01). AUCs calculated with $εder$ were also significantly lower than AUCs calculated than $εtp$ for flexor carpi ulnaris (FCU) (p < 0.01), flexor carpi radialis (FCR) (p < 0.01), and flexor pollicis brevis (FPL) (p < 0.01) during wrist knock.

Parameter sensitivity was found to vary based on task (Supplemental Figure S1 available in the Supplemental Materials on ASME Digital Collection). Isometric tasks had higher sensitivities compared to isokinetic tasks. When examining muscle parameters, all tasks were sensitive to $Fmmax$ and $Lmopt$⁠, but insensitive to α (Fig. 5), while bone mass was primarily influential during isokinetic tasks. The insensitivity to changes in α and bone mass can be seen graphically as horizontal lines while plotting the AUCs versus parameter changes (Fig. 6). Isokinetic tasks showed increased noise in AUC values when changing parameter values, which was most likely a result of variations in timing and recruitment during a motion task (Fig. 6).

All tasks were most sensitive to $Lmopt$⁠, with the exception of lateral pinch. For example, isometric wrist knock was most sensitive to $Lmopt$ with ECRL's AUC changing by 0.95% for a 1% change in $Lmopt$ (Fig. 6). This $εtp$ indicates that the ECRL's AUC changed by 28.5% over the tested range. In contrast, lateral pinch was most sensitive to changes in $Fmmax$ with a 1% change in $Fmmax$ decreasing the adductor pollicis oblique head's (ADPo) AUC by 0.8%. This change indicates that the ADPo's AUC decreased by 40% over the 50% change in $Fmmax$⁠. Isokinetic elbow knock saw the highest sensitivity to α with a 1% change decreasing the triceps brachii long head's (TRIlong) AUC by 0.02%. This means that over the range of tested percent changes, a 40% change, the TRIlong's AUC only changed by 0.8%.

Our study demonstrates that sensitivity coefficients computed using a sensitivity function better account for parameter uncertainty than sensitivity coefficients computed using the two-point method and linear regression slope. Sensitivity coefficients pulled from a sensitivity function, $εder$⁠, had lower RMSE across all tasks and muscles. This is likely a result of $εder$ looking at the local slopes of individual percent changes, instead of using a single sample, $εlrs$⁠, or averaged sample, $εtp$⁠. By examining the local trend, $εder$ accounts for nonlinear regions created by the Hill-type muscle curves or other nonlinearities of the simulation pipeline, such as the nonlinear cost function used by computed muscle control. This is illustrated by changes in how $Lmopt$ influences muscle activations during isometric elbow flexion (Fig. 3). In contrast, $εlrs$ and $εtp$ result in higher errors the more a region diverges from the assumed linear behavior. For this reason $εlrs$ and $εtp$were appropriate in previous studies that focused on small linear regions around the baseline value, but may not be appropriate when studying large parameter changes at the population-level. This interpretation is supported by previous work that found linear scaling of Hill-type parameters is most accurate for small parameter changes [43]. While calculating $εder$ is more computationally costly, it provides the benefit of (i) identifying how the Hill-type parameters influence model predictions across the population and (ii) develops sensitivity as a function.

This study also demonstrates that parameter sensitivity is task-dependent for greater physiological ranges for hand, wrist, and elbow tasks. Similar to previous studies [13,14,16], we found sensitivity was both task and muscle dependent. For example, ECRL and ERCB were primarily sensitive to both isometric and isokinetic wrist tasks. However, ECRL and ECRB were twice as sensitive to $Lmopt$ and $Fmmax$ during isometric wrist extension compared to isokinetic wrist extension/flexion. While previous studies at the lower [12–14,16] and upper [17] limb found that simulations were always more sensitive to $Lmopt$ than $Fmmax$⁠, this trend was reversed for muscles during the lateral pinch task. In direct contrast to Redl et al. [16], who stated the sensitivity to $Fmmax$ was negligible during gait, lateral pinch was most sensitive to $Fmmax$⁠, highlighting that sensitivities at the lower limb cannot be extrapolated to the upper limb. These deviations highlight reported sensitivity trends do not hold for all tasks, thereby suggesting that only task-specific muscles may require personalization.

Pennation angle may not need to be personalized in subject-specific upper limb models. Varying α within physiological ranges did not change predicted muscle activations regardless of the task. This finding is consistent with Carmichael and Liu who found models within physiological bounds were insensitive to α; although they also found that models beyond physiological bound (i.e., 35% change in $Fmmax$⁠) were sensitive to changes in α [17]. In upper limb models, this insensitivity likely occurs because of the small angle approximation. The small magnitudes of α within upper limb muscles means that even extreme percentage changes result in similar, small values. Thus, when using the cosine of α to convert muscle force to tendon force, small changes in α result in near zero changes in the calculation.

Our findings are affected by some limitations. We chose to only change a single parameter within each model. This decision was made so that computational power could be used to have smaller parameter step sizes compared to the two points used by previous studies. Parameter interactions could change a task's overall sensitivity, but we expect sensitivity to still be task dependent. We also limited our tasks to the hand, wrist, and elbow, so individual parameter trends may not be valid at the shoulder. Likewise, we could not feasibly consider all postures, tasks, and movements at the hand, wrist, and elbow. Instead, we focused on representative tasks to understand how sensitivity changes across different joints using fundamental motions.

This study presents a comparison of (1) methods to compute task-specific parameter sensitivity for subject-specific modeling and (2) task-specific parameter sensitivities for simulated tasks at the hand, wrist, and elbow. The methods described can be used to account for how inaccuracies in subject-specific measurements propagate through experimental analyses using musculoskeletal computer models and simulations. Performing an a priori task sensitivity analysis could inform which model parameters to prioritize when constructing personalized models. Future work can continue to develop the method of obtaining sensitivity coefficients using numerical derivatives by exploring how step size, parameter range, and sample size may influential prediction errors, thereby informing standardization of this method in biomechanics.

• The National Institute of Biomedical Imaging and Bio-engineering (NIBIB) of the National Institutes of Health (NIH) under a Trailblazer Award (Award No. R21 EB030068; Funder ID: 10.13039/100000070).

• National Science Foundation Graduate Research Fellowship Program (NSF GRFP) (Award No. 1315138; Funder ID: 10.13039/100000001).

The data and information that support the findings of this article are freely available at online².

$Fmmax$ =

maximum isometric force

$Lmopt$ =

optimal fiber length

$α$ =

pennation angle

$ε$ =

sensitivity coefficient

$εtp$ =

sensitivity coefficient calculated using two points

$εlrs$ =

sensitivity coefficient calculated using linear regression

$εder$ =

sensitivity coefficient calculated using local derivatives

$M0$ =

initial value of measured output

$ΔM$ =

change in value of measured output

$P0$ =

initial value of model parameter

$P−$ =

negative value change of model parameter

$P+$ =

positive value change of model parameter

$ΔP$ =

change in value of model parameter

Abbreviations

DELT1 =

anterior deltoid

DELT2 =

medial deltoid

DELT3 =

posterior deltoid

SUPSP =

supraspinatus

INFSP =

infraspinatus

SUBSC =

subscapularis

TMIN =

teres minor

TMAJ =

teres major

PECM1 =

pectoralis major clavicular

PECM2 =

pectoralis major medial

PECM3 =

pectoralis major inferior

LAT1 =

latissimus dorsi superior

LAT2 =

latissimus dorsi medial

LAT3 =

latissimus dorsi inferior

CORB =

coracobrachialis

TRIlong =

triceps brachii long head

TRIlat =

triceps brachii lateral head

TRImed =

triceps brachii medial head

ANC =

anconeus

SUP =

supinator

BIClong =

biceps brachii long head

BICshort =

biceps brachii short head

BRA =

brachialis

BRD =

brachioradialis

ECRL =

extensor carpi radialis longus

ECRB =

extensor carpi radialis brevis

ECU =

extensor carpi ulnaris

FCR =

flexor carpi radialis

FCU =

flexor carpi ulnaris

PL =

palmaris longus

PT =

pronator teres

PQ =

pronator quadratus

FDSL =

flexor digitorum superficialis digit 5

FDSR =

flexor digitorum superficialis digit 4

FDSM =

flexor digitorum superficialis digit 3

FDSI =

flexor digitorum superficialis digit 2

FDPL =

flexor digitorum profundus digit 5

FDPR =

flexor digitorum profundus digit 4

FDPM =

flexor digitorum profundus digit 3

FDPI =

flexor digitorum profundus digit 2

EDCL =

extensor digitorum communis digit 5

EDCR =

extensor digitorum communis digit 4

EDCM =

extensor digitorum communis digit 3

EDCI =

extensor digitorum communis digit 2

EDM =

extensor digiti minimi

EIP =

extensor indicis proprius

EPL =

extensor pollicis longus

EPB =

extensor pollicis brevis

FPL =

flexor pollicis longus

APL =

abductor pollicis longus

APB =

abductor pollicis brevis

FPB =

flexor pollicis brevis

OPP =

opponens pollicis

ADPt =

adductor pollicis tranverse head

ADPo =

adductor pollicis oblique head