The impact of anatomical uncertainties on the predictions of a musculoskeletal hand model – a sensitivity study

Abstract

Outputs of musculoskeletal models should be considered probabilistic rather than deterministic as they are affected by inaccuracies and estimations associated with the development of the model. One of these uncertainties being critical for modeling arises from the determination of the muscles' line of action and the physiological cross-sectional area. Therefore, the aim of this study was to evaluate the outcome sensitivity of model predictions from a musculoskeletal hand model in comparison to the uncertainty of these input parameters. For this purpose, the kinematics and muscle activities of different hand movements (abduction of the fingers, abduction of the thumb, and flexion of the thumb) were recorded. One thousand simulations were calculated for each movement using the Latin hypercube sampling method with a corresponding variation of the muscle origin/insertion points and the cross-sectional area. Comparing the standard hand to simulations incorporating uncertainties of input parameters shows no major deviations in on- and off-set time point of muscle activities. About 60% of simulations are located within a ± 30% interval around the standard model concerning joint reaction forces. The comparison with the variation of the input data leads to the conclusion that the standard hand model is able to provide not over-scattered outcomes and, therefore, can be considered relatively stable. These results are of practical importance to the personalization of a musculoskeletal model with subject-specific bone geometries and hence changed muscle line of action.

Keywords:

• Musculoskeletal
• hand
• AnyBody
• Latin hypercube sampling
• sensitivity analysis

Introduction

The outcome of musculoskeletal models, e.g. joint moments and forces, is mostly reported as single values with a deterministic character. The incorporation of any uncertainties within the model is often left unexamined (Cook et al. 2014). When inaccuracies associated with model development are taken into considerations, simulation results have to be interpreted as a single case out of many possible outcomes. These uncertainties may arise either from simplifications in the model (Gagnon et al. 2011) or inaccuracies in underlying input variables as a function of data extraction (Żuk and Pezowicz 2016). The process of identifying such parameters is usually associated with complex measurements (Scheys et al. 2009; Taddei et al. 2012), which are based on an average error greater than zero (Scheys et al. 2008).

Critical variables for modeling are the muscle parameters, such as the line of action (Raikova and Prilutsky 2001), physiological cross-sectional area (PCSA), and thus the maximum force capacity of muscles (Herzog 1992). These are crucial parameters for determining muscle activities and joint reaction forces. Nussbaum et al. (1995) has already shown that such uncertainties in the muscle line of action result in a change of more than 100 N in the joint reaction force within the spine.

If the model to be analyzed is based only on a single anatomical data set, the inaccuracy in the determination of muscle parameters is reduced to measurement errors as the data is not diversified. This lack of diversity may also lead to a systematic error if the subject's anatomical preconditions in patient-specific modeling deviate significantly from the cadaver, which the model is based on. To achieve ideal population-based models, a more comprehensive anatomical data set is necessary, which can be gained by including a large number of different cadavers when determining muscle parameters. The advantage of including a larger number of cadavers increases the representativeness of the sample means and reduces their susceptibility to be skewed by outliers. However, scattering of the respective parameters also leads to considerable uncertainties in parameters. Another disadvantage depicts the fact that there is high variability within human performance measures within a population, and hence only a few people can be defined as an average person (Smith et al. 2014). However, generalized musculoskeletal models' overall purpose should not be the ability to represent only an average person but to be applicable to a large number of subjects with idiosyncratic characteristics through patient-specific scaling.

The recently published hand model by Engelhardt and Melzner et al. (2020) within the AnyBody™ Modeling System (AMS) (Anybody, Aalborg, Denmark) is based on a general data set, which includes data of 16 donors (Havelková et al. 2020). This model consists of 80 extrinsic and intrinsic hand muscles with 31 degrees of freedom and has been shown to correlate highly with moment arm studies from prior literature (Engelhardt and Melzner et al. 2020). Since this model's anatomical data has some unavoidable known indeterminacy, the outcome sensitivity of subject-specific model predictions (muscle activities and joint reaction forces) compared to the input parameters should be researched. Thus, this study aimed to investigate the sensitivity and stability of the model and at the same time determine whether the model can be applied to a broad population beyond the average hand. Due to its increasing importance in using various modern technologies, such as smartphones (Gustafsson et al. 2010), we particularly focused on the thumb.

Materials and methods

The PCSA and the muscles’ line of action are crucial parameters in determining acting muscle forces and hence the acting joint reaction forces. In order to evaluate the sensitivity of the recently published musculoskeletal hand model from Engelhardt and Melzner et al. (2020) to these parameters, the underlying anatomical data (Havelková et al. 2020) was varied according to their known measurement uncertainties – see Appendix A, supplementary material. Subsequently, the simulated on- and off-set time points of the muscle activities were compared with experimentally derived electromyographic activities. Further, the effects on the change of joint reaction forces within the hand were investigated.

Experiment set-up

The test subject (weight = 78.3 kg, hand length = 20.5 cm) was informed about the measurement procedure in advance and provided written consent confirming voluntary participation. Then, the test subject performed one trial with the abduction of all fingers, one trial with the abduction of the thumb, and one trial with the flexion of the thumb. Each trial consisted of five repetitive movements, where the hand was in pronation for the abduction of the fingers and neutral position otherwise.

The movements' kinematics were captured with a camera-based motion capture system (Vicon©, Vicon Motion Systems Ltd., UK) using a marker-set including the upper arm, forearm, and all fingers. The EMG activities of 10 extrinsic and intrinsic (see Table 1 and Figure 1) muscles were captured via surface EMG (Delsys Trigno IM and Delsys Trigno Mini, Delsys® Inc., MA, US.) following the recommendation of sensor placement according to Criswell and Cram (2011) and Barbero et al. (2012). Prior to the sensor placement, the skin was cleaned with an alcoholic solution to minimize the acquired signal's influences due to dead skin cells, salt, or grease. The muscles were selected according to the movements to be investigated (abduction of all fingers/thumb and flexion of the thumb), such that in addition to the primarly active intrinsic muscles, extrinsic muscles stabilizing the wrist were also covered. For normalization of the EMG data, maximum voluntary contraction (MVC) measurements following Kendall (2010) were performed beforehand.

Figure 1. Experimental set-up with marker and EMG sensor placement of the forearm and hand of the test subject.

Table 1. Measured EMG signal of intrinsic and extrinsic muscles.

Extrinsic muscles	Sensor	Protocol
Flexor carpi radialis (FCR)	T-IM	**
Flexor carpi ulnaris (FCU)	T-IM	*
Extensor carpi radialis (ECR)	T-IM	*
Extensor carpi ulnaris (ECU)	T-Min	*
Flexor digitorum superficialis (FDS)	T-Min	*
Extensor digitorum (ED)	T-IM	*
Intrinsic muscles	Sensor	Protocol
Abductor pollicis brevis (APB)	T-Min	**
Flexor pollicis brevis (FPB)	T-Min	**
First dorsal interosseous (FDI)	T-Min	*
Abductor digiti minimi (ADM)	T-Min	**

The placement of the sensors was done according to Criswell and Cram (2011) * and Barbero et al. (2012). ** For the measurement, Trigno IM (T-IM) and Trigno Mini (T-Min) sensors from Delsys were applied.

Simulation

The gathered kinematic data was fed into the musculoskeletal simulation software AnyBody™ (AnyBody Technology A/S, Aalborg, Denmark, V. 7. 3). Then, the model was scaled according to the anthropometric data of the test subject to allow for patient-specific simulation, and simple muscle models are applied. The change of input parameters and the simulation were performed externally using the python module AnyPyTools (Lund et al. 2019) to conduct the subsequent sensitivity study. The PCSA and the origin, as well as the insertion points of all extrinsic and intrinsic hand/wrist muscles were amended according to their measurement inaccuracies published by (Havelková et al. 2020) – see Appendix A, supplementary material.

A modification in the PCSA is directly incorporated into the cost function G, which has to be minimized by the algorithm $$G=\mathrm{}\sum _i{\left(\frac{f_i}{N_i}\right)}^p$$ whereby p is the polynomial coefficient symbolizing the muscle synergy, f_i is the applied force, and N_i is the strength of a muscle, which is linear proportional to the PCSA. Hence, the magnitude of the PCSA influences muscle recruitment and, thus, the joint reaction forces to be determined. The outputs are also influenced by changing the muscles' origin and insertion points, resulting in a resulting alteration in the moment arm.

For the PCSAs, origin, and insertion points, a normal distribution (with the expectation value µ being the value of the standard model and standard deviation (SD) of the anatomical data – see Appendix A, supplementary material) was assumed. Further, we made the following assumptions regarding origin and insertion points:

• The origin and insertion area was split up according to the number of muscle representatives within the model.

• Since the anatomical dataset provided only the area in a two-dimensional plane, the SD in the third dimension was assumed to be 1 mm.

• To not exceed the limits of the origin and insertion areas while randomly creating input parameters, the boundaries of the area were set to ±3SDs. This restriction assures that origin and insertion points stay at 99.7% within the anatomical area.

Based on this set-up, 1.000 simulations were calculated for each trial. We used the Latin-Hypercube sampling (LHS) technique (Helton and Davis 2003) to perform sensitivity analysis, which is a standard method in biomechanics as it provides improved sampling efficiency compared to the otherwise commonly used Monte Carlo analysis (Valero-Cuevas et al. 2003; Reinbolt et al. 2005).

Data processing

For the EMG signal, a root mean square with a window size of 50 ms was applied, the baseline offset was subtracted, and the signal was normalized according to the MVC measurements. To compare the numerical with the experimentally measured muscle activity, a compilation of corresponding representatives of a muscle was created whereby only the maximum value at each time step of the grouped muscle elements was taken into account. A detailed list of enveloped muscles can be found in the Appendix B, supplementary material.

Following prior research, which has used an evaluation of on- and off-set timing of muscle activities to validate musculoskeletal models (Pieri et al. 2018), the on- and off-set time points of the experiment, the standard model, and all simulations run by the LHS method were compared. This comparison is particularly interesting since a larger variation of the on- and off-set timing would imply a questionable experimental validation or instability of the model.

For automatization of on- and off-sets detection, change point detection (Killick et al. 2012) as implemented by the MatLab (The MathWorks Inc., Natick, USA) FindChangePts() function was used. Based on this output, the difference between each simulation and the experiment was evaluated. Next, the joint reaction forces acting on the wrist and all finger joints were calculated using the inverse dynamic approach of the AMS. This inverse dynamic approach provides a computation of the muscle forces and the resulting joint reaction forces, but it can not be used directly for motion prediction. Since a relatively small range of motion in the distal and proximal interphalangeal joint took place during the conducted movements, only the carpometacarpal (CMC) and the metacarpophalangeal (MCP) joint were considered in the investigation of the joint reaction forces.

For analysis of the model output and measured data, Rebba and Mahadevan (2008) showed that a metric r $$r=P(-\epsilon <m-c<\epsilon )$$ where m-c is the difference between experimental data and numerical simulation, and $\epsilon$ is the accuracy requirement, which is more suitable than the hypothesis testing using the p-value. Therefore, the r metric was used for comparing the experimental data, the standard model, and the outcome of all LHS simulations regarding the on- and off-set time points. This metric was also applied to comparing joint reaction forces, although the model cannot be checked against experimental data. Instead, the results of the LHS simulations were compared to the standard model.

Results

As the aim of this study was to test the sensitivity and stability of the musculoskeletal hand model, Figure 2 exemplifies the extent to which variation in input parameters can affect changes in muscle activity (Figure 2(a)) and joint reaction forces (Figure 2(b)). Figure 2 displays the activity of the FDI muscle and the joint reaction force of the thumb's CMC joint during the abduction of all fingers. In both cases, the standard model is within the range of one SD of the LHS simulations. It can be stated that the variation in the joint reaction force is higher than that of the individual muscle activities. All muscle activities and joint reaction forces of each motion can be found in the Appendix C–E, supplementary material.

Figure 2. (a) Muscle activity of the FDI muscle and (b) the joint reaction force of the CMC joint of the thumb during the abduction of the fingers. The red dashed line represents the outcome of the standard model, whereas the blue line indicates the mean of the LHS simulation with its SD as a shaded area. The force of (b) is given in percent body weight (BW).

The mean values of the difference between all calculated on and off time points for each movement are displayed in Table 2. In contrast to the fingers' abduction, during movements of the thumb, mainly FPB and APB muscles were addressed. Hence, only these muscles were taken into account for calculating the on and off times. Comparing the mean values of the LHS simulations and the standard model, only a slight shift of the on-off times can be observed.

Table 2. The average time difference of the on- and off-set time points of the muscles between the experimental measured and simulated predictions during (a) the abduction of all fingers, (b) the flexion of the thumb, and (c) the abduction of the thumb.

(a) On-/Off-Sets	ED	ECR	FCR	FCU	FDI	ADM	APB	FPB	FDS	ECU
Standard model	160	185	181	136	190	178	90	106	283	100
Mean – LHS	167	258	209	170	194	201	88	121	276	100
SD – LHS	39	180	173	177	169	213	27	63	245	25
(b) On-/Off-Sets							APB	FPB
Standard model							309	258
Mean – LHS							311	269
SD – LHS							26	61
(c) On-/Off-Sets							APB	FPB
Standard model							204	224
Mean – LHS							203	225
SD – LHS							7	22

The standard model is compared to the mean value of the LHS simulations with modified input parameters and its SD. The on- and the off-sets time difference is given in millisecond.

The joint reaction forces (normalized to the body weight) acting on the wrist and finger joints for each movement are shown in Table 3. When determining the joint reaction forces, only times when the muscles of the standard model were active were taken into account and averaged. The highest acting forces were found in the wrist joint, whereby the CMC joint of the index finger was also exposed to high loads during thumb movements. Considering the joint reaction forces, the standard model also stays within one SD of the LHS simulations except for the index finger's MCP joint during the abduction of the thumb. The coefficient of variation (CV) – which is defined as the ratio of SD and mean - varies between 9% (MCP_thumb – abduction thumb) and 56% (CMC_finger4 – abduction thumb) for individual joints.

Table 3. The average joint reaction forces of the wrist, CMC, and MCP finger joints during the (a) the abduction of all fingers, (b) the flexion of the thumb, and (c) the abduction of the thumb.

(a)
	Wrist	CMCthumb	CMCfinger2	CMCfinger3	CMCfinger4	CMCfinger5	MCPthumb	MCPfinger2	MCPfinger3	MCPfinger4	MCPfinger5
Standard model	35.4	6.1	14.0	5.2	2.8	14.0	1.2	3.5	5.0	4.5	9.3
Mean – LHS	31.8	5.7	11.4	5.3	2.7	13.1	1.3	3.4	4.6	4.1	9.2
SD – LHS	9.4	2.4	4.8	2.2	0.5	2.5	0.3	1.2	1.4	0.6	1.7
CV – LHS	0.29	0.42	0.42	0.41	0.19	0.19	0.23	0.34	0.3	0.14	0.18
(b)
	Wrist	CMCthumb	CMCfinger2	CMCfinger3	CMCfinger4	CMCfinger5	MCPthumb	MCPfinger2	MCPfinger3	MCPfinger4	MCPfinger5
Standard model	31.6	4.1	14.6	6.4	4.1	7.8	2.5	4.0	5.4	6.8	7.8
Mean – LHS	30.2	4.0	11.6	8.4	4.4	8.0	2.5	4.0	5.0	5.8	8.0
SD – LHS	10.2	1.6	3.9	3.3	1.6	2.9	0.4	1.0	1.7	2.4	2.9
CV – LHS	0.34	0.39	0.33	0.39	0.37	0.36	0.15	0.24	0.35	0.42	0.37
(c)
	Wrist	CMCthumb	CMCfinger2	CMCfinger3	CMCfinger4	CMCfinger5	MCPthumb	MCPfinger2	MCPfinger3	MCPfinger4	MCPfinger5
Standard model	56.8	12.6	33.6	7.3	1.7	6.0	11.5	14.9	7.4	2.8	5.8
Mean – LHS	57.9	13.7	33.6	7.7	2.8	5.1	11.7	18.6	7.8	3.8	4.9
SD – LHS	8.9	2.3	4.2	2.3	1.6	2.1	1.0	2.8	2.3	1.6	2.1
CV – LHS	0.15	0.17	0.12	0.3	0.56	0.42	0.09	0.15	0.29	0.43	0.42

The reaction force of the standard model is compared to the mean value of the LHS simulations with modified input parameters, its SD, and corresponding CV. The joint reaction forces are given in percent body weight (BW).

The results for applying the metric proposed by Rebba and Mahadevan (2008) to the outcome of the on- and off-set time difference and the joint reaction forces, respectively, are shown in Figure 3–5 for each movement. Figures 3(a)–5(a) show the probability of finding an outcome of the predicted joint reaction forces of the LHS simulations within a specific interval around the calculated value of the standard model (all forces were normalized to the standard model). Figures 3(b)–5(b) show the probability that the time difference between the experimental and numerically calculated data is within a specified time interval for the standard model and the LHS simulations. Considering the joint reaction forces, almost 70% of all simulations stay within a 30% interval around the standard model for the abduction of all fingers. Notable exceptions are the CMC joint of the thumb and index finger, for which only 60% of results occur within this interval (Figure 3(a)).

Figure 3. Probability for the abduction of all fingers (a) to find an outcome of the predicted joint reaction forces of the LHS simulations within a specific interval around the calculated value of the standard model (all forces are normalized to the standard model) (b) that the time difference between the experimental and numerical calculated data is within a specific time interval for the standard model and the LHS simulations.

Figure 4. Probability for the flexion of the thumb (a) to find an outcome of the predicted joint reaction forces of the LHS simulations within a specific interval around the calculated value of the standard model (all forces are normalized to the standard model) (b) that the time difference between the experimental and numerical calculated data is within a specific time interval for the standard model and the LHS simulations. Only the APB and FPB muscles were considered.

Figure 5. Probability for the abduction of the thumb (a) to find an outcome of the predicted joint reaction forces of the LHS simulations within a specific interval around the calculated value of the standard model (all forces are normalized to the standard model) (b) that the time difference between the experimental and numerical calculated data is within a specific time interval for the standard model and the LHS simulations. Only the APB and FPB muscles were considered.

Regarding the thumb flexion, approximately 60% of simulation outcomes stay within a 30% interval around the standard model (Figure 4(a)); whereas for the thumb abduction, it is about 70% except for the joints of digitus IV and V (Figure 5(a)). Regarding the probabilities of finding an on- or off-set point in a specific time interval, the results of the LHS simulations follow the trend of the standard model (see Figure 3(b)–5(b)).

Discussion

The aim of this study was to examine the sensitivity and stability of the musculoskeletal hand model published by Engelhardt and Melzner et al. (2020) and thus to determine whether the model can be applied to a broad population outside the average hand for patient-specific modeling. Examining the comparison between experimental data, the standard model, and LHS simulations suggests that the hand model is able to provide convergence for a broad range of input data and, therefore, can be considered relatively stable.

Combining all joints' CV values results in an average scattering of 30% around the expected value. The probability that the calculated joint reaction forces of the LHS simulations are located within an interval of ±30% around the value of the standard model is about 60% for the three investigated movements. Even when comparing the experimental data with the numerical models, the highest deviation between the prediction of the standard model and the mean value of the LHS simulations is 73 ms (abduction of all fingers - ECR muscle) and the probability curves of Figures 3(b)–5(b) show no substantial deviations. Further, all SDs are smaller than equal to approximately 250 ms, which is within an acceptable range once electromechanical delays (which are about 75 ms, according to Seth and Pandy (2007)) are taken into account.

Investigating joint reaction forces, the level of particular forces within the hand during the different movements is of particular interest. Especially in the abduction of the thumb, the wrist, and the index finger's CMC joint seem to be exposed to high forces (Table 3(c)). These increased forces are especially interesting, as this movement is mainly used in the operation of modern multimedia devices. However, because there are no in vivo measured force data of these joints - compared to other joints like the shoulder (Bergmann 2008) - these results cannot be validated quantitatively.

Taking a closer look at the outliers in the probability metrics of Figures 3–5, it is evident that outliers are mainly finger joints that were rarely involved in the respective movement, such as the CMC and MCP joint of fingers 3–5 for the abduction of the thumb or the CMC joint of finger 3 for the flexion of the thumb. Only the MPC joint of the index finger (during the abduction of the thumb) and the CMC joints of the thumb and index finger (during the abduction of all fingers) seemed to be sensitive to changes in input parameters, although the movements addressed these joints. One possible reason for these joints' increased sensitivity could be the large number of intrinsic muscles running through the area between the index finger and the thumb.

When evaluating the scattering of the model’s outcome, the distribution of input parameters needs to be considered. The scattering magnitude is difficult to quantify in terms of the muscles’ insertion/origin area, but for PCSA, this is quite manageable. The anatomical data set from (Havelková et al. 2020) can be used to recalculate the average SD of PCSA values for all hand muscles, which is 38% of the corresponding mean. Weighting the individual muscles according to their size results only in a slightly smaller weighted average SD of 36.5%. Assuming a normal distribution of the underlying PCSA values, roughly 68% of the values are within an interval of ±36.5%. Comparing this measure of dispersion to that of the numerical outcomes shows that the model predictions do not scatter more than the input parameters (with the limitation of reference to the PCSA).

The LHS method used in this study has improved sampling efficiency compared to the Monte Carlo analysis, which has the disadvantage of requiring a large number of simulation runs. Accordingly, the number of 1.000 simulations seems to be sufficient compared to the Monte Carlo studies of Valente et al. (2014) and Myers et al. (2015) with 500 and 3.000 simulations, respectively.

Although the results show that the model seems to be relatively stable concerning changes in input parameters, this stability depends on the type of movement performed, which should also be considered in future analyses. Therefore, an attempt was made to depict a correspondingly comprehensive range of relevant movements being essential for modern everyday life.

This study has important implications for patient-specific scaling of musculoskeletal models. A common approach to patient-specific modeling is to scale the model’s force-specifically, which is analogous to changing the PCSA. Alternatively, it is also possible in the AMS to morph patient-specific bone geometries, which results in a shift of the muscle attachment points. The effects of both methods were illustrated in this study.

However, some limitations of this research need to be noted. The comparison of experimentally gained and numerically determined muscle activities introduces crosstalk, which is associated with the measurement of EMG signals (Konrad 2005). The kinematic inaccuracies that inevitably occur during the measurement can also impact model predictions (Myers et al. 2015). In combination with the inverse dynamic, these kinematic inaccuracies can cause unrealistic joint torques approach according to Bailly et al. (2021). In addition, due to the time-independent nature of the static optimization, an inverse dynamic calculation prevents the activation dynamics from being taken into account. Both limitations can be circumvented by the forward dynamic approach. Another limitation of the used inverse dynamics approach is the impossibility of any direct kind of motion prediction.

Further, some assumptions about the distribution of input parameters had to be made when the anatomical data set did not provide exact data (e.g. for the SD of the third dimension of the origin/insertion area). In addition, it must be taken into account that this study did not investigate the direct effect of a change in muscle moment arms and thus the resulting variation in joint reaction forces. Rather, a change in the muscle line of action indirectly caused changes in the muscle moment arms. Finally, it is challenging to derive generally valid data from the measurement of only three different movements of a single person. However, we argue that these results provide valuable preliminary information, which can be useful to future investigations.

Conclusion

The study presented a sensitivity and stability analysis of the musculoskeletal hand model by Engelhardt and Melzner et al. (2020) by changing the input parameters of the muscles' line of action. The results indicate the model's corresponding stability, which provides a reasonable basis for the patient-specific simulation of a broad population regarding the hand.

Disclosure statement

All authors of this manuscript have no conflict of interest.

Additional information

Funding

This work was supported by the project n. 182 Obstetrics 2.0 – virtual models for injury prevention during childbirth, realized within the frame of the Program INTERREG V-A: Cross-border cooperation between the Czech Republic and the Federal State of Germany Bavaria, Aim European Cross-border cooperation 2014–2020. This work is further supported by the Bavarian Academic Forum (BayWISS) Doctoral Consortium ‘Health Research’ and is funded by the Bavarian State Ministry of Science and Arts.