## Abstract

Patient-specific biomechanical models including patient-specific finite-element (FE) models are considered potentially important tools for providing personalized healthcare to patients with musculoskeletal diseases. A multi-step procedure is often needed to generate a patient-specific FE model. As all involved steps are associated with certain levels of uncertainty, it is important to study how the uncertainties of individual components propagate to final simulation results. In this study, we considered a specific case of this problem where the uncertainties of the involved steps were known and the aim was to determine the uncertainty of the predicted strain distribution. The effects of uncertainties of three important components of patient-specific models, including bone density, musculoskeletal loads and the parameters of the material mapping relationship on the predicted strain distributions, were studied. It was found that the number of uncertain components and the level of their uncertainty determine the uncertainty of simulation results. The ‘average’ uncertainty values were found to be relatively small even for high levels of uncertainty in the components of the model. The ‘maximum’ uncertainty values were, however, quite high and occurred in the areas of the scapula that are of the greatest clinical relevance. In addition, the uncertainty of the simulation result was found to be dependent on the type of movement analysed, with abduction movements presenting consistently lower uncertainty values than flexion movements.

## 1. Introduction

Patient-specific biomechanical modelling is considered a promising approach for providing personalized healthcare [1]. Musculoskeletal diseases are of particular interest and importance, because they are in many cases related to mechanical forces and the structural performance of tissues. Patient-specific finite-element (FE) models are among the most important biomechanical models that could be used for providing patients with personalized healthcare. For example, both generic and patient-specific FE models have been used for design of orthopaedic implants [2,3], predicting the longevity and performance of orthopaedic implants [4–6], predicting the risk of fracture in osteoporotic patients [7–10] and orthopaedic surgery planning [11–13].

A typical procedure for creating patient-specific FE models of the musculoskeletal system includes the following steps [14]: (1) acquisition of images using clinical image acquisition techniques such as computer tomography (CT) or magnetic resonance imaging (MRI), (2) segmentation of images and creation of FE meshes, (3) assignment of patient-specific mechanical properties based on density values and empirical correlations between bone mineral density and Young's modulus, and (4) application of patient-specific musculoskeletal loads and boundary conditions that are estimated from patient-specific musculoskeletal models [15]. Owing to the increased availability of imaging modalities, development of efficient image processing programmes, increased computational power and new approaches in image-based modelling such as the use of statistical shape and appearance models [16], patient-specific FE modelling of the skeletal system has become feasible during the last decade and a number of researchers have investigated the different aspects of patient-specific FE modelling for different applications [17–21].

Patient-specific FE models are assumed to be more accurate than generic FE models, because their elements are specific to the patient for whom the FE model is created. Nevertheless, not much information is available in the literature regarding the actual accuracy of patient-specific FE models of the musculoskeletal system. One important issue that may adversely affect the accuracy of patient-specific FE models is the accumulation of uncertainty through several steps of model generation. The four major steps that were mentioned in the previous paragraph are associated with certain degrees of uncertainty and error that may propagate to simulation results.

Patient-specific FE modelling and the accuracy of such models have received a lot of attention recently [22–26]. Some of the above-mentioned aspects of FE modelling, such as the effects of density values, density–modulus relationships [22,24] and their uncertainty [23,26] on simulations results, have been studied. However, other than in a few isolated studies (e.g. [27]), there is not much information available regarding the interplay of the different sources of uncertainty in patient-specific FE modelling and how that influences the ultimate results.

One specific case in this regard is when the uncertainties of all procedures involved in the generation of patient-specific FE models are known. In this paper, we consider this specific case and study the uncertainty of final simulation results. From the four above-mentioned steps involved in patient-specific FE modelling of bones, we consider the uncertainties associated with three major steps, namely densitometry (part of step 1), material mapping (step 3) and load estimation (step 4). The uncertainty associated with imprecise segmentation of anatomical shapes (part of step 1 and entire step 2) is not considered here partly owing to the complexities associated with considering variations in anatomical shapes and partly because the other three sources of uncertainty are considered to be more important. In general, the shape of the scapula and other bones can be quite accurately segmented from CT images, particularly given that the gradient of bone density could be used as auxiliary information for detecting bone boundaries. Moreover, uncertainties in detecting the boundary of bones do not significantly change the shape of bones and are therefore expected to have minimal effects on the accuracy of the predicted strain distribution.

We have chosen to model the scapula in this study because patient-specific FE modelling of the upper extremity is generally more challenging than that of the lower extremity. For example, uncertainty may be higher in measurement of the scapula density distribution, as certain parts of the scapula are extremely thin and practically ‘see-through’. Therefore, the partial volume effects may play an important role and result in inaccurate density values. In this study, we used high-resolution CT scanning to minimize the effects of thin parts on the segmented shape. In addition, information regarding the density gradient was used during to improve the sensitivity of the segmentation process to bone boundaries. It is important to realize that the thin parts of the scapula are away from the point of application of the joint reaction force that is by far the largest force. As most surgical procedures are performed close to the glenohumeral joint, the areas close to the joint are much more important than the thin areas that are far from the joint.

In addition, not many material mapping relationships are available for the upper extremity. One therefore needs to use the relationships that are developed for bones of other anatomical sites. As these relationships are known to be site specific [28], the material mapping procedure may be less accurate for the upper extremity bones. Finally, the movements of the upper extremity are more complex. It is therefore more challenging to obtain an accurate estimation of the musculoskeletal loads of the upper extremity.

We use a Monte Carlo approach for assessment of the accuracy of patient-specific FE models of the scapula. First, a ‘reference model’ is generated based on one cadaveric scapula for which a complete set of modelling parameters are measured for both patient-specific musculoskeletal and patient-specific FE modelling. Uncertainties are then introduced to the different components of the model based on Gaussian distributions with known average and standard deviations. Subsequently, Monte Carlo simulations are carried out to analyse the effects of uncertainties of individual components of the FE model on the accuracy of the resulting simulation results.

## 2. Material and methods

The study was carried out in two steps. In the first step (§2.1), a ‘reference’ patient-specific model of one scapula was created. In the second step (§2.2), the reference model was used to study the effects of uncertainty in individual components of the reference model on the accuracy of FE simulation results.

### 2.1. The reference models

A patient-specific model of one scapula was created as the reference model. The cadaveric scapula that was used for creating the reference model belongs to a donor (male, 57 years) for whom the parameters of the upper extremity soft-tissue structure such as muscle volume, optimum fibre length, tendon length and pennation angle were directly measured. The relationships between the soft-tissue structures and bones, such as muscle and ligament attachment sites, were also measured. The bony structure was scanned using a clinical CT scanner (Siemens, SOMATOM Definition Flash, construction diameter: 230 mm) with an isotropic resolution of ≈0.6 mm. This constitutes a unique dataset of consistent soft- and hard-tissue parameters that could provide a very high level of accuracy in patient-specific biomechanical modelling of the scapula. To the best of the authors’ knowledge, there is currently no second scapula for which this rich dataset of consistent soft- and hard-tissue parameters are available. This particular scapula is therefore an obvious choice for the reference patient-specific model of the scapula. A detailed description of soft- and hard-tissue parameters and their measurement techniques can be found in earlier studies [29–34].

Creating a patient-specific FE model of the scapula requires information about patient-specific musculoskeletal loads. Estimation of musculoskeletal loads is normally performed either using simple mass–spring–damper models [35–37] that only give the generic trends of loading or using large-scale musculoskeletal models [38–40] that provide the detailed loading conditions of joints and muscles. We used a patient-specific large-scale musculoskeletal model of the upper extremity, namely the Delft shoulder and elbow model (DSEM), for estimating the detailed musculoskeletal modelling of the scapula. This model has been validated against experimental measurements of the joint force using instrumented shoulder implants [29]. As previously mentioned, the entire set of soft- and hard-tissue parameters that are needed for musculoskeletal modelling of the upper extremity were measured for this particular donor. The DSEM includes 31 muscles and muscle parts and 139 muscle lines of action and is capable of estimating the glenohumeral joint reaction force as well as the magnitude and direction of all muscle forces for any given movement of the upper extremity. However, one needs to provide kinematic data regarding these movements so that the positions, velocities and acceleration of the different body parts can be calculated and related to forces through Newton's laws. The kinematic data were measured for a number of healthy subjects whose anthropometrical parameters were close to those of the donor. The average of those kinematic measurements was assumed to represent the likely kinematic data of the donor. The musculoskeletal loads were calculated for two specific movements, namely 90° abduction and 90° flexion. For every one of those movements, the joint reaction forces as well as all individual muscle forces were calculated.

The geometry used in FE modelling was obtained by segmentation of the scapula from the acquired CT images using Mimics 14.0 (Materialise, Leuven, Belgium). Subsequently, the obtained geometry was meshed (figure 1*a*) with 4-node linear tetrahedron elements (edge length = 1.2 mm, 349 283 elements). The number of elements was chosen based on a convergence study. The number of elements used for discretization of the scapula geometry was increased, in several discretization steps, to 1.84 million elements. The maximum strain values and the strain energy of the FE model were analysed. It was observed that the FE solution converged within 1% of the calculated values once around 220 000 elements were used. The heterogeneous material properties of the scapula were calculated based on bone density values. The following phantom-calibrated relationship between the Hounsfield units (HU) of the CT scanner and apparent density (*ρ*) was used:
2.1

Young's modulus was then related to calculated values of apparent density using Morgan's relationship [41] 2.2

Poisson's ratio was fixed at 0.3.

The musculoskeletal loads calculated using the DSEM for both considered movements, namely 90° abduction and 90° flexion, were applied on the geometry of the scapula as traction forces spread over a few elements, and two FE models were created for those two movements. We therefore obtained two ‘reference models’, each corresponding to one of the considered movements. The boundary conditions were applied at three points within the bone geometry, namely points AA (acromial angle, most posterior), TS (trigonum scapula, most medial) and AI (inferior angle, most inferior) as specified in figure 1*b,c*. Boundary conditions were specified to constrain rigid body motions of the geometry and facilitate convergence of the FE solution. For AA, all degrees of freedom were constrained. The displacements in the *y-* and *z*-directions were constrained for TS. The displacements of AI were only constrained in the *z*-direction (figure 1*b,c*). Both FE models corresponding to two considered movements were created in Abaqus v. 6.10 (Simulia, Providence, RI) and were solved using an implicit FE solver (Abaqus Standard). The stress and strain distributions resulting from the applied musculoskeletal loads (figure 1*c*) were calculated. In a recent study [42], we demonstrated the accuracy of the reference model in a relatively indirect way. In that study, the reference model was coupled with a bone tissue adaptation model to predict the density distribution of the scapula. It was shown that, when coupled with the bone tissue adaptation model, the reference model predicts the CT-measured density distribution quite accurately.

### 2.2. Monte Carlo study

As already discussed, the uncertainty in patient-specific FE models may originate from different steps involved in generation of the models. Three sources of uncertainty were considered in this study, namely uncertainty in the estimation of musculoskeletal loads, uncertainty in the measurement of bone density and uncertainty in the parameters of the material mapping relationship.

We assumed that the uncertainty in the modelling parameters is normally distributed with an average value equal to the value used in the reference models and a standard deviation that is a certain percentage of the average values. For all three sources of uncertainty, four levels of uncertainty were considered corresponding to 5%, 10%, 20% and 30% of the average parameter values.

In this context, the parameters of equation (2.2), which were used in the Monte Carlo study, were picked from two Gaussian distributions with average values equal to the values presented in equation (2.2) for parameters *a* and *b*, and standard deviations equal to 5%, 10%, 20% or 30% of the average values.

As for bone density, the density of every element in the FE models was picked from a Gaussian distribution with an average value equal to the density of that element in the reference models and a standard deviation equal to 5%, 10%, 20% or 30% of that average value.

The musculoskeletal loads were modified in a similar manner. The directions of all forces were kept the same as in the reference model. However, the magnitude of all musculoskeletal loads including the glenohumeral joint reaction force and all individual muscle forces were modified. The new value of every load was drawn from a Gaussian distribution with an average value equal to the magnitude of that force in the reference models and a standard deviation equal to 5%, 10%, 20% or 30% of that average value. As there is no guarantee that the scapula is in equilibrium once musculoskeletal loads are subjected to such perturbations, the equilibrium of the scapula was re-established using an optimization procedure that minimized the sum of the absolute values of modifications needed for individual loads to re-establish force and moment equilibriums.

For every case of uncertainty analysis, the simulations were carried out 50 times (50 realizations) with 50 different values drawn from all involved Gaussian distributions. It was then possible to calculate the average and standard deviation of the uncertainty values caused by imposed perturbations.

The simulations were carried out for the cases where only one of the above-mentioned components of the patient-specific modelling process was uncertain, for the cases where all three components were uncertain, and for the cases where two out of three components were uncertain.

In order to analyse the deviations of the perturbed models from the reference models, we had to restrict ourselves to a few results of the FE analysis. The maximum (major) and minimum (minor) principal strains were chosen for uncertainty analysis, because they are important variables and present a fairly complete picture of strain and stress distributions.

The deviations of the perturbed models from the reference models were quantified using three uncertainty metrics. For every perturbed model, the maximum principal strains, *E*_{max}, of all elements were compared with those of the elements of the reference model. The elements were first ranked according to the absolute difference between the value of their maximum principal strain and that of their corresponding element in the reference model. The top 1% of the elements was then selected for analysis. Percentage error, Error_{1,1%}, was calculated between the maximum principal strain values of elements *i*, *E*_{max,i,perturb}, and those of the corresponding elements of the reference model, *E*_{max,i,ref},
2.3where *i* = *1*, …, *n* and *n* (i.e. *N*/100) is 1% of the total number of elements in the FE models (*N*). Similar error values were calculated for minimum principal strain, *E*_{min}. The error values, Error_{1,1%}, calculated in this way quantify the maximum deviation of the maximum tensile and compressive strains of perturbed models from those of the reference model. Similar values were calculated for the top 0.1% and 0.01% of the elements (see the electronic supplementary material).

In order to also quantify the *average* difference between the perturbed and reference models, the average value of the percentage error, Error_{2}, was also calculated
2.4Finally, the root mean squared error (RMSE) between the perturbed and reference models was also calculated
2.5Similar to Error_{1,1%}, the values of Error_{2} and RMSE were calculated also for minimum (minor) principal strain. All three quantifiers of uncertainty were calculated for both considered movements and, thus, for both reference models. As every perturbed model was simulated 50 times (50 different realizations), 50 different values were obtained. The average and standard deviation of the uncertainty metrics were then calculated. It is important to remember that, despite their names, Error_{1,1%}, Error_{2} and RMSE are actually ‘uncertainty’ values and not actual error values of the model.

When discussing the ‘range of uncertainty’, we refer to the average value of uncertainty metrics plus its standard deviation. Error_{1,1%} and Error_{2} values show the percentage difference between the simulation results of perturbed models and those of reference models. Using Error_{1,1%} and Error_{2} values, one could relate a certain percentage of uncertainty in model inputs to the resulting percentage of uncertainty in simulation results.

The Error_{1,1%} values calculated for different levels of uncertainty (5%, 10%, 20% and 30%) were compared with each other using the analysis of variance (ANOVA) followed by Tukey–Kramer post-hoc analysis. The Error_{1,1%} values calculated for flexion and abduction were compared with each other using Student's *t*-test. A significance threshold of *p* < 0.05 was used for both types of statistical analysis.

## 3. Results

In general, the overall pattern of strain distribution did not drastically deviate from that of the reference models (figure 2). The most notable changes were observed in the extreme strain values (figure 2). An important observation was that the maximum deviations of perturbed models from the reference model, quantified by Error_{1,1%}, were comparable to the percentage of uncertainty in the inputs of the perturbed models (figures 3–7 and table 1). However, the average values of uncertainty over the entire geometry of the scapula, quantified by Error_{2}, were in general much lower than the percentage of uncertainty induced in the components of the perturbed models (tables 2 and 3). In general, the range of Error_{1,1%} values calculated for both considered movements had considerable overlap (figures 3–7). However, statistical analysis showed that the uncertainty values were, in a relatively large number of cases, significantly different between both types of movements (table 4). Whenever there was a significant difference between the Error_{1,1%} values calculated for flexion and abduction, the Error_{1,1%} values of abduction were lower than those of flexion (table 4).

Increasing the level of uncertainty did not necessarily result in significantly higher levels of uncertainty in simulation results. When the different sources of uncertainty were independently considered (figure 3), it was clear that the level of uncertainty caused by the various sources of uncertainty is different (figure 3). When only musculoskeletal loads were uncertain (figure 3*a*), the range of Error_{1,1%} values were more or less similar to or slightly smaller than the percentage of uncertainty present in the applied muscle and joint reaction forces. Moreover, there was a significant difference between Error_{1,1%} values calculated for the four different levels of uncertainty in musculoskeletal loading (table 5). When the parameters of the material mapping relationship (equation (2.2)) were inaccurate (figure 3*b,c*), the maximum percentage of uncertainty in calculated strains, i.e. Error_{1,1%} values, was sometimes markedly smaller than the percentage of uncertainty present in the parameters of equation (2.2). When both parameters of equation (2.2) were compared (cf. figure 3*b,c*), it was observed that uncertainties in parameter *b* generally translate to smaller uncertainty values in calculated strains (figure 3*c*) when compared with uncertainties in parameter *a*. Roughly speaking, the percentage of uncertainty in calculated strain values was less than half of the percentage of uncertainty in the values of parameter *b* (figure 3*c*). In addition, Error_{1,1%} values calculated for 5% uncertain *a* values were not significantly different from the Error_{1,1%} values calculated for 10% uncertain *a* values (table 5). In the case of parameter *b*, Error_{1,1%} values were not significantly different between 5% and 10% and between 20% and 30% of uncertainty (table 5). Uncertainty in bone density had a larger impact on the simulation results than uncertainty in the parameters of the material mapping relationship. The impact of bone density was, however, somewhat smaller than that of musculoskeletal loads (figure 3*d*).

When all three components of perturbed models were uncertain, the percentage of uncertainty in the calculated strain values was much larger than the percentage of uncertainty in individual components of the perturbed models (figure 4). In the four cases considered here, i.e. 5%, 10%, 20% and 30% uncertainty, the percentage of uncertainty in calculated strains was around 1.5 times larger than the percentage of uncertainty in the individual components of the perturbed models (figure 4). Moreover, the Error_{1,1%} values calculated for all four cases (5%, 10%, 20% and 30%) were significantly different from each other (table 5).

When bone density and material mapping parameters were uncertain, the additional uncertainty caused by uncertainty in estimation of musculoskeletal loads was moderate in most cases (figure 5). The gradual increase in the level of uncertainty of musculoskeletal loads from 0% to 30% increased the maximum uncertainty value in the calculation of strain values from 25.50 ± 10.78% (figure 5*a*) to 31.36 ± 12.07% (figure 5*e*). For smaller levels of uncertainty in estimation of musculoskeletal loads, i.e. 5–10%, and smaller levels of uncertainty in other components of the model, e.g. 5%, the range of uncertainty values in the calculation of strain hardly changed (compare figure 5*a* with figure 5*b,c*).

By assuming the musculoskeletal loads to be perfectly accurate (figures 6 and 7), one could study the isolated effects of uncertainty in bone material properties on the uncertainty of calculated strain values. When both parameters of the material mapping relationship were uncertain but density values were assumed to be perfectly accurate, the range of uncertainty in the calculation of strain values was more or less equal to the range of uncertainty in the parameters of the material mapping equation (figure 6*a*). Similar to the case of musculoskeletal loads, when the parameters of the material mapping equation were uncertain, the additional uncertainty because of uncertain density values was moderate (figures 6 and 7).

## 4. Discussion

The results of this study show that, as expected, the accuracy of patient-specific FE models depends on the accuracy of the steps that need to be followed for generation of the models. However, there is a major difference between the average and maximum values of uncertainty. The average uncertainty values stay relatively low (≈20 ± 9%) even when all three components considered here are up to 30% inaccurate (tables 1 and 2). However, the maximum uncertainty values were found to be more than 1.5 times larger than the percentage of uncertainty in the individual components of the model (figure 4). This value (1.5) is dependent on the percentage of elements selected for uncertainty analysis. If the top 0.1% or 0.01% of the elements had been used for uncertainty analysis instead of the top 1%, this value would have been higher and around 2% (see the results presented in the electronic supplementary material).

Another important result is that the level of uncertainty may also be dependent on the type of movement for which the analysis is performed. In many cases including the case where all sources of uncertainty were simultaneously present, the Error_{1,1%} values calculated for abduction were consistently lower than those of flexion. This trend was consistently observed for all cases where there was a significant difference between the Error_{1,1%} values of abduction and those of flexion. As different movements of the upper extremity result in different patterns of muscle recruitment and activity, the magnitude and orientation of muscle forces change with the type of movement. That could result in different modes of mechanical loading and, subsequently, different levels of uncertainty in simulation results.

One needs to note that the areas close to the glenohumeral joint are of particular importance for the planning of joint replacement surgery and analysis of implant loosening. It is therefore important to have an accurate estimation of stress and strain distributions close to the glenohumeral joint. As the joint reaction force is by far the largest force in the set of musculoskeletal loads applied on the scapula [29], the largest stress and strain values are generally observed close to the glenohumeral joint. Given that the largest deviations from the reference models occur for the areas with maximum and minimum strain values (figure 2), the largest uncertainty values usually occur close to the glenohumeral joint. One could therefore conclude that maximum uncertainty values, i.e. Error_{1,1%} values, are more relevant for evaluation of the accuracy of patient-specific FE models than the Error_{2} values that give estimates of the average uncertainty values within the entire scapula.

### 4.1. Uncertainty levels and relevance for patient-specific finite-element modelling

In this study, the uncertainty levels of all three components of patient-specific models were assumed to be between 5% and 30%. We tried to relate those levels of uncertainty to the typical levels of uncertainty found in the literature. As for bone density, it has been shown that the uncertainty associated with bone density measurement using a clinical CT scanner could be as high as 10% [43]. This is within the range of variations (5–30%) considered here. As for musculoskeletal loads, the joint reaction force predicted with the DSEM was validated against measurements of the joint reaction force using an instrumented shoulder prosthesis [29]. The validation study showed that when the model was made (partially) patient specific, the prediction of the joint reaction force changed by about 10% for abduction movements, suggesting an uncertainty level of around 10% [29]. These values are also within the range of uncertainties (5–30%) considered in this study.

In a recent study [44], specimen-specific values of the parameters *a* and *b* were determined for different cadaveric femora using an optimization process. The values found for *a* varied between 9307 and 15 673, while *b* varied between 0.87 and 1.40. Such large deviations in the values of parameters *a* and *b* are in line with the values considered here, i.e. up to 30%. Although *a* and *b* values could be quite accurately determined for any *in vitro* specimen, there is no such possibility when dealing with patient-specific FE models used in clinical settings. One therefore has to resort to some sort of representative density–modulus relationships that could be used for different individuals. Owing to large variations between different individuals, any representative density–modulus relationship will result in relatively large errors for the individuals for whom the coefficients of the density–modulus relationship are far from the representative values. Moreover, the modulus–density relationships vary from one anatomical site to another and using the same relationships for bones from different anatomical locations could result in additional errors. This is an important point given that not all bones are equally well studied. For example, we are not aware of any studies that determine the density–modulus relationship of the scapula and the other bones of the upper extremity. Therefore, FE modelling of the scapula is currently performed using the density–modulus relationships developed for other bones, and that could cause additional inaccuracy in patient-specific FE models.

The range of variability considered in this study refers to patient-specific FE models created for individual patients. Subject- or specimen-specific FE models could be created very accurately for a few subjects or specimens in certain ideal conditions. For example, one could accurately measure the density distribution of cadaveric bone specimens *in vitro*, because there are no effects of soft tissues surrounding the bone, and the bone can be scanned using a high radiation dose for a long time. However, such practices are not possible when creating ‘patient’-specific FE models. Additionally, estimation of musculoskeletal loads is much less accurate when dealing with patients when compared with the accuracy that can be achieved for subjects for whom the musculoskeletal models have been developed. The range of uncertainties encountered in creating patient-specific FE models could therefore be much larger than those encountered in some of the specimen- and subject-specific FE models. As the ultimate goal of all such studies is to translate patient-specific models to clinical settings, it is important to consider the actual conditions that are likely to be encountered in clinical settings. That is why relatively large levels of uncertainty were also examined here.

If the above-mentioned values for the accuracy of density measurements and musculoskeletal load estimation could be considered ‘typical accuracy values’ and if the uncertainty of the material mapping relationship could be considered to be between 5% and 30%, one could estimate the ‘typical’ range of uncertainty that is caused by the three components of the models considered here. The results of this study suggest that such a ‘typical’ range of uncertainty would be between ≈10% and ≈35% (figure 5*c*).

### 4.2. Discussion on modelling assumptions

One of the most important assumptions used in this study is the assumption that the uncertainty associated with the different components of patient-specific models are normally distributed. The normal distribution assumption is a common assumption in statistical uncertainty analysis where no additional information is available regarding the distribution of uncertainty values. However, it is important to remember that the uncertainties associated with different components of patient-specific models (1) are not necessarily normally distributed and may conform to other statistical distributions and (2) may include a systematic component. The Monte Carlo analysis performed in this study cannot capture the effects of systematic uncertainties on the results of FE simulations and is limited to random uncertainties that are normally distributed. The presented results are, nevertheless, useful for analysis of the accuracy of patient-specific FE models, because there is currently no information available regarding the effects of uncertainty of individual components of models on the accuracy of patient-specific simulation results. The results presented here can therefore serve as a first approximation of uncertainties associated with patient-specific FE modelling of the scapula.

### 4.3. Limitations of the study

The limitations of this study are either related to generation of reference models or to the methods used in the uncertainty analysis. The limitations of the methods used for generation of reference models are very similar to those of other FE models of bones discussed elsewhere [14,45,46]. There are, however, a few remaining points that need to be taken into account. The most important point is regarding the material mapping relationships. As previously discussed, there is currently not much data available in the literature regarding the relationship between the apparent density and mechanical properties of the scapula. For example, a recent study that reviewed such relationships [28] found no study dedicated to the upper extremity bones. We therefore had to use a material mapping relationship that has been developed using lower extremity bones. As the relationship between the apparent density and the mechanical properties of bone is site dependent [41], this could introduce some errors in the reference models. It is therefore necessary to develop material mapping relationships that are specific to the scapula. Full-field strain measurement techniques such as digital image correlation that can be used for characterization of engineering materials [47–49] and biological tissues are good candidates for this purpose. If the strain distribution during mechanical testing of a number of scapulae is measured and those scapulae are also modelled using FE models, one could optimize the parameters of the material mapping relationships such that the difference between the measured and predicted strain distributions is minimized.

The limitations regarding the uncertainty analysis are twofold. As previously mentioned, the first limitation is owing to the assumption that uncertainties are normally distributed. Second, only three components of patient-specific FE models were considered in this study. The uncertainties associated with the image segmentation and mesh generation stage were not considered. It should, however, be noted that the uncertainties related to the mesh generation stage are relatively easier to control as one could use finer meshes for discretization of the bone geometry until calculated strains converge within sufficiently small tolerances; see the details of the convergence study presented above. As for the uncertainties related to the geometry of the model, it is not expected that small (i.e. sub-millimetre) uncertainties in geometry segmentation could result in large uncertainties in prediction of density distribution. The most important source of uncertainty related to the geometry of the scapula may be due to possible modulations between the overall scapula geometry and FE modelling errors. In that case, the effects of accuracy of individual components on the accuracy of patient-specific FE models may be shape specific. At present, it is not clear whether any strong modulation exists between the scapula shape and FE modelling results.

### 4.4. Implications for other areas of research

The uncertainty analysis performed here was focused on patient-specific FE modelling of the scapula. There are, however, many other applications that are closely related to the study presented here. Patient-specific FE modelling of the upper extremity in pathological conditions such as modelling of an implanted shoulder joint is an important example. Many of the steps involved in FE modelling of the pathological cases are similar to the steps studied here. The results presented here may therefore be at least partially valid for pathological cases. Moreover, patient-specific models of the lower extremity such as hip and knee joints are also somewhat similar to those of the upper extremity. The presented methodology and results may be therefore helpful in understanding the uncertainties involved in patient-specific FE modelling of the lower extremity in both physiological and pathological conditions.

## 5. Conclusion

The effects of uncertainty in three important components of patient-specific FE models of the scapula on the uncertainty of predicted strain predictions were studied using Monte Carlo simulations. The average uncertainty within the entire geometry of the scapula was found to be lower than the uncertainty of individual components of the model even when several components were highly (i.e. up to 30%) uncertain. In addition, the maximum uncertainty values within the geometry of the scapula were found to be generally much higher than the average values and occurred in the areas that are more relevant for analysis of pathological conditions. The percentage of uncertainty in calculated strain values was found to be dependent on the number of inaccurate components of the model and the level of uncertainty of individual components. It was also observed that the uncertainty values depend on the type of movement for which musculoskeletal loads are calculated and applied. In all significantly different cases, uncertainty values calculated for abduction were consistently lower than those of flexion.

- Received December 8, 2013.
- Accepted January 21, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.