- Email: [email protected]
Subject-specific modeling of the scapula bone tissue adaptation
Journal of Biomechanics 46 (2013) 2434–2441
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Subject-specific modeling of the scapula bone tissue adaptation Gianni Campoli a, Harrie Weinans a,b, Frans van der Helm a, Amir A. Zadpoor a,n a Department of Biomechanical Engineering, Faculty of Mechanical, Maritime, and Materials Engineering, Delft University of Technology (TU Delft), Mekelweg 2, Delft 2628 CD, The Netherlands b Department of Orthopaedics, Erasmus University Medical Center, Rotterdam, The Netherlands
art ic l e i nf o
a b s t r a c t
Article history: Accepted 12 July 2013
Adaptation of the scapula bone tissue to mechanical loading is simulated in the current study using a subject-specific three-dimensional finite element model of an intact cadaveric scapula. The loads experienced by the scapula during different types of movements are determined using a subject-specific large-scale musculoskeletal model of the shoulder joint. The obtained density distributions are compared with the CT-measured density distribution of the same scapula. Furthermore, it is assumed that the CTmeasured density distribution can be estimated as a weighted linear combination of the density distributions calculated for different loads experienced during daily life. An optimization algorithm is used to determine the weighting factors of fourteen different loads such that the difference between the weighted linear combination of the calculated density distributions and the CT-measured density is minimal. It is shown that the weighted linear combination of the calculated densities matches the CTmeasured density distribution better than every one of the density distributions calculated for individual movements. The weighting factors of nine out of fourteen loads were estimated to be zero or very close to zero. The five loads that had larger weighting factors were associated with either one of the following categories: (1) small-load small-angle abduction or flexion movements that occur frequently during our daily lives or (2) large-load large-angle abduction or flexion movements that occur infrequently during our daily lives. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Shoulder joint Scapula Bone tissue adaptation Finite element Patient-specific Optimization process
1. Introduction The morphology of bone tissue is partially determined by the mechanical loads it experiences during daily activities (Frost, 1988; Ruff et al., 2006; Turner, 1998). Using adaptation theories (Carter et al., 1987; Cowin and Hegedus, 1976; Huiskes et al., 1987), one could theoretically connect the musculoskeletal loads caused by daily activities to the measured morphology of bone tissue. However, there are two major difficulties in establishing this connection. First, there is no easy way for non-invasive measurement of musculoskeletal loads. Second, bones are loaded differently in the various movements that we carry out in our daily living, and it is not clear how every movement contributes to the measured morphology of bone. These two difficulties exist in the study of both lower- and upper-extremity bones. For example, many researchers have studied bone tissue adaptation of lower extremity bones in both physiological (Campoli et al., 2012; Turner et al., 2005; Weinans et al., 1992; Weinans and Prendergast, 1996) and post-operative conditions (Bernd-Arno et al., 2009; Huiskes et al., 1987; Lengsfeld et al., 2002). The adaptation of the trabecular
n
Corresponding author. Tel.: +31 15 2786794; fax: +31 15 2784717. E-mail addresses: [email protected], [email protected] (A.A. Zadpoor). 0021-9290/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jbiomech.2013.07.024
bone is also extensively studied (Biewener et al., 1996; Jang and Kim, 2008; Shefelbine et al., 2005; Tanck et al., 2006; Tsubota et al., 2009). However, the two above-mentioned difficulties have not been addressed before. In most cases, researchers have simply assumed generic and highly simplified loading conditions. As for upper-extremity bones, even studies that use generic and simplified loading conditions are rare. Indeed, no numerical model of glenoid bone remodeling was available until recently (Sharma et al., 2009, 2010; Sharma and Robertson, 2013). In this study, we present a methodology for establishing the connection between the movements of daily living and the morphology of the scapula and overcoming both above-mentioned difficulties. Once the above-mentioned difficulties are overcome, we could study the effects of individual movements on bone tissue adaptation and separate the effects of different movements from each other. Scapula was used for this study, because distinguishing between different movements is easier in the scapula. That is due to the fact that upper-extremity loading is more varied as compared to the lower-extremity loading that is dominated by gait. Even though currently available models (Sharma and Robertson, 2013) could predict the density distribution of the scapula quite accurately, they do not differentiate between the effects of different types of movements on bone tissue adaptation. We use a validated (Nikooyan et al., 2008, 2010) large-scale musculoskeletal model of the shoulder and elbow, namely the Delft
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
2435
difference between weighted linear summation of the predicted density distributions and the CT-measured density distribution is minimal. The proposed methodology is used for relating the measured density distribution of a scapula to the musculoskeletal loads estimated for the different movements of the same individual using the musculoskeletal model.
Shoulder and Elbow Model (DSEM) (van der Helm, 1994a,b; van der Helm and Pronk, 1995; van der Helm and Veenbaas, 1991). For every movement, this model can provide accurate prediction of detailed musculoskeletal loads including joint reaction and muscle forces and their point of application. In order to overcome the second obstacle, we propose to use an optimization procedure for linking the loads estimated by the musculoskeletal model and their resulting density distributions to the density distribution measured using computed tomography (CT). In this scheme, the CT-measured density distribution is assumed to be explained by a weighted linear combination of the density distributions caused by a number of movements. The optimization process finds the weighting factors of different movements such that the
2. Materials and methods 2.1. Computed tomography One cadaveric scapula of a male donor (57 years) was CT-scanned using a clinical scanner (Siemens, SOMATOMs Definition Flash, construction diameter: 230 mm)
c 44 45
79 58 7 7
78
114
113 52
116
8 47 5 40 51 49
C 1L
8
0 5 69
111
61
126 8 9197257 62 96
53 54 03 1 80
83
JGH
115
lcs TS_ 63 RTS
82
LT
46 112
0 JAC898 789 85 84 6 s 8 lc AA_
G 11072 1
75
106 129 125
98
124
118 74
101
119
122 5 10
100 120
64
109
65
121 6969
73
104
123 110 107
67
89 70 66
131 0 132
108 I lcs RAAI_
Fig. 1. The FE mesh used in simulations (a), bony landmarks and boundary conditions (b) as well as an example of applied joint and muscle loads (c). The inferior–superior and anterior–superior distances of the glenoid fossa, see von Schroeder et al. (2001), were respectively 34 and 27 mm.
2436
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
with an isotropic resolution of 0.6 mm. The CT images were imported into Materialise Mimics and were semi-automatically segmented. Segmentation of the scapula is relatively difficult due to the very thin ‘see-through’ parts of the scapula. After pre-selecting bone areas using a threshold-based technique, a more accurate segmentation was obtained by manually ‘wrapping’ the masked image of the scapula blade cross-section to fill the gaps and voids between disconnected pixels, followed by ‘erosion’ of the pixels that were too far away from the actual image of the thinnest areas of the scapula blade. The segmented scapula was subsequently meshed (Fig. 1a) using 4-node linear tetrahedron elements (edge length¼ 1.1 mm, 349,283 elements).
2.2. Loads and boundary conditions A schematic drawing of the complete methodology is presented in Fig. 2. The details of every step are described below. DSEM was used for estimation of the detailed loading of the scapula while performing two different types of movements: abduction and flexion. These two movements were chosen because of their high frequency in daily living and because many other more complex movements can be related to these two movements. The parameters of the DSEM such as the muscle attachment sites, muscle volume, anatomical parameters and landmarks, shape of the thorax, and several other parameters were directly measured for the same donor.
Kinematic data, movement 1 Kinematic data, movement 2 ...
Musculoskeletal model
Kinematic data, movement n ... Forces, movement 1
Forces, movement 2
Forces, movement n
Finite element model ...
Density distribution, movement 1
Density distribution, movement 2
Density distribution, movement n ...
CT-measured density distribution
w2
wn
New values of weighting factors
w1
Objective function value
No Optimization algorithm
Stopping criteria met?
Stopping criteria of optimization are met, for example, when the objective function does not drop further (within a specified tolerance) for a certain number of iterations.
Yes
End. Fig. 2. A flowchart of the methods used in the current study for establishing the connection between different movements and the morphology of the bone.
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441 These parameters and their definitions can be found in references (Nikooyan et al., 2010, 2011; van der Helm, 1994a,b; van der Helm and Pronk, 1995; van der Helm and Veenbaas, 1991). DSEM is a large-scale musculoskeletal model of the upper extremity and includes thorax, clavicle, scapula, humerus, ulna, and radius. The parameters of the model such as muscle mass, muscle attachments sites, physiological cross-section area of muscle, force-length parameters, optimum fiber lengths, and tendon lengths were measured through an extensive cadaveric study. The model includes 31 muscles and muscle parts and 139 muscle lines of action. Using these directly measured subject-specific parameters, the glenohumeral joint reaction forces as well as muscle forces were estimated with the DSEM for seven different angles of abduction and flexion, namely 201, 401, 601, 801, 1001, 1201 and 1401. The kinematic data needed for solution of the musculoskeletal model were measured in a population of healthy subjects whose anatomies were similar to the donor's anatomy. The average of the measured kinematics data (position of body segments) was used as input to the inverse dynamics musculoskeletal model. By differentiating the vector of position changes with respect to time, one can calculate the velocities and accelerations of different body segments. The mass and inertia properties of the body segments are also known and one could use Newton's equations for calculating internal forces. A physiological muscle load-sharing criterion is used for determining the contribution of every muscle to the calculated forces. The forces estimated with DSEM were applied on the surface of the generated mesh (Fig. 1c). Transferring the musculoskeletal loads to the finite element (FE) mesh required slightly modified points of application of the forces (0.3–2.0 mm distance of muscle forces from the scapula geometry), because the geometry used in the DSEM is a simplified geometry and had slight deviations from the CT-measured geometry. The force and moment equilibrium were re-established after mapping the points of application of musculoskeletal loads onto the FE mesh. The distance between the points of force application given by the DSEM with all the nodes of the FE mesh was calculated. For every load, the node of the mesh that had the least distance from the point of force application was used for applying the force. Similar to other studies (Suaŕ ez et al., 2009, 2012), an optimization algorithm was used to calculate the smallest changes in the magnitude of musculoskeletal forces that will re-establish the equilibrium. Loads were applied as traction surface loads distributed over the area of one single element. The naming convention of bones and bony landmarks that is used here is in agreement with the International Society of Biomechanics (ISB) recommendations (Wu et al., 2005). Three nodes corresponding to three standard bony landmarks of the scapula, namely AA (Anglus Acromial, most posterior), TS (Trigonum Scapula, most medial), and AI (Angulus Inferior, most inferior), were partially or totally constrained. For AA, all degrees of freedom were constrained. The superior–inferior and anterior–posterior displacements (respectively y and z directions in Fig. 1b) were constrained for TS. The displacements of AI were only constrained in the anterior–posterior direction (z direction in Fig. 1b).
2.3. Bone tissue adaptation model The mechanical properties were based on the following correlation between the apparent density and elastic modulus (Morgan et al., 2003): E ¼ aρ ; b
a ¼ 6950; b ¼ 1:49
ð1Þ
Poisson's ratio was fixed at 0:3. The density values were calculated from CTmeasured gray values using a phantom-calibrated linear relationship. A strain-based tissue adaptation model (Campoli et al., 2012; Weinans et al., 1992) was used. In this model, the change in bone density, ρ, is linked to the first principal strain, ε1 : Δρ ¼ B
ε1 k Δt ρ
ð2Þ
The strain threshold to trigger bone remodeling is assumed to be around 1000 micro-strain (Frost, 1998). The maximum density value measured by CT is 2.2 g/cm3. The reference stimulus, k, is calculated as the ratio of the strain threshold to the maximum density value: 1000 10 6/2.2¼0.0005 cm3/g. The parameter B regulates the speed of remodeling. The results of the bone tissue adaptation model were found to be largely independent from the choice of the parameter B. To speed up calculations, a B value of 1000 g2/cm6 s was used in the current study. There are two additional parameters: ρmax and ρmin . Similar to other studies (Sharma et al., 2009, 2010), the density values that are calculated to be large than ρmax were registered as ρmax whereas the ones that were smaller than ρmin were registered as ρmin . These enforced boundaries ensure that the density values remain within the physiological range. One can determine the minimum and maximum density values, ρmax and ρmin , based on the values observed in the CT images (see values in Fig. 4a). An implicit FE solver (ABAQUS) was used. The simulation started from a uniform density value (¼ 0.8 g/cm3) for all elements of the FE model. In every iteration, the first principal strain and the current density of all elements were sent to a subroutine that calculated the change in the density of all elements. The density and, thus, mechanical properties of the elements (Eq. (1)) were accordingly updated (Campoli et al., 2012; Zadpoor et al., 2013). The tissue adaptation model was solved for 25 iterations to ensure that the mean absolute difference between
2437
Fig. 3. The mean of absolute difference between the density distributions of two consecutive iterations divided by the first density distribution plotted vs. iteration number. The results are presented only for the movement (1201 flexion) that had the greatest difference between the density distributions calculated for two first consecutive iterations. two consecutive density distributions divided by the first density distribution was less than 0.1% (Fig. 3). In total, we carried out 14 simulations: one independent simulation for every one of the above-mentioned movements (Fig. 2). 2.4. Linear combination of the predicted density distributions The measured density distribution has resulted from adaptation of the bone tissue in response to a variety of mechanical loads that occur in daily living. We assumed that the measured density distribution is a weighted linear combination of the density distributions that result from the individual loads experienced during daily living: ρCT ¼ ∑i wi ρi
ð3Þ
where ρCT is the CT-measured vector of density distribution, ρi is the vector of the density distribution that results from the load i, and wi is the weight associated with the load i. The weight wi may be a function of the frequency and/or intensity of the load. The weight factors associated with the considered loading cases were determined through a gradient-based nonlinear optimization algorithm with the following objective function: f ðxÞ ¼ jjρCT ∑i wi ρi jj
ð4Þ
The algorithm evaluated the objective function at the end of all iterations (Fig. 2). Even though it is obviously not possible to capture all the loads experienced during our daily activities, one has to include the most important loads that are most likely to influence bone tissue adaptation. The loads that were considered in the current study were the fourteen sets of loads that were estimated by the DSEM for seven different angles of flexion and seven different angles of abduction (see Section 2.2). As previously explained, flexion and abduction are very important movements due to their high frequency of occurrence in daily living and because many other movements can be related to flexion and abduction.
3. Results The FE results showed that no single abduction or flexion movement could accurately explain the CT-measured density distribution (Fig. 4a). In order to make a quantitative comparison between the predicted, ρFE , and measured, ρCT , density distributions, the difference between the two distributions was calculated for all the elements of the FE mesh. For every element j, the prediction error was calculated as ρ ρ j;FE j;CT ð5Þ Errorj ¼ 100 ρj;CT The mean value of the prediction error is, indeed, the mean percent error of the absolute difference between the predicted and CT densities. Two movements, i.e. 801 flexion (Fig. 4b) and 801 abduction (Fig. 4c) are among the movements that best predict the
2438
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
Fig. 4. The CT-measured density distribution (g/cm3) (a) as well as density distribution (g/cm3) calculated for 801 abduction (b), 801 flexion (c), and combined loading case (d) displayed in six different planes. The distance between two consecutive planes is 20 mm.
Table 1 The mean and standard deviation of the absolute and percentage difference between the predicted and CT-measured density distributions for the fourteen considered movements as well as for the linear combination of the fourteen density distributions. Abduction
Mean error (%) SD, error (%) Mean error (g/cm3) SD, error (g/cm3)
Flexion
Combi
201
401
601
801
1001
1201
1401
201
401
601
801
1001
1201
1401
—
36.7 12.3 0.396 0.201
32.7 13.4 0.351 0.195
31.4 13.7 0.337 0.194
29.2 14.0 0.312 0.191
28.8 13.9 0.309 0.188
33.4 13.4 0.361 0.200
36.6 13.0 0.398 0.212
36.9 12.5 0.400 0.206
35.2 13.0 0.381 0.203
33.4 13.4 0.360 0.200
32.6 13.8 0.351 0.202
33.3 13.9 0.360 0.206
34.1 13.6 0.370 0.206
37.4 12.7 0.407 0.212
15.9 13.4 0.170 0.150
Abbreviations: Combi: weighted linear combination of the considered density distributions, SD: standard deviation.
CT-measured density distribution (Table 1). In general, the mean values of the calculated prediction errors are between 29% and 38%.
The weighted linear combination of density distributions better matches the measured density distribution (compare Fig. 4d and a). The results of the applied optimization procedure (Table 2) suggest
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
2439
Table 2 Optimized values of the weighting factors of the fourteen considered movements. These weight factors are the ones for which the minimum value of the objective function is calculated by the optimization algorithm. Abduction
Flexion
201
401
601
801
1001
1201
1401
201
401
601
801
1001
1201
1401
w1 0.21
w2 0.20
w3 0
w4 0
w5 0.48
w6 0
w7 0
w8 0.18
w9 0
w10 0.07
w11 0
w12 0
w13 0
w14 0.38
that the contributions of certain movements are much higher than those of the other movements. In particular, five movements were found to dominate the linear combination due to their much larger weight factors (Table 1, wi 4 0), while nine movements did not contribute much to the calculated density distribution (Table 1, wi ≈0). The mean and median prediction errors were respectively 15.9% and 12.8% for the optimized linear combination of density distributions. The mean prediction error is about half of the smallest mean prediction error calculated for single movements (Table 1). The median of the prediction error (12.8%) is also much smaller than the values calculated for single movements (30.2–38.1%).
4. Discussion It was observed that no single movement could successfully explain all measured density variations. Nevertheless, the density distribution is better predicted in certain areas of the scapula even for single movements. These areas are particularly high-density areas that experience high strain values in various movements. Because of their highly loaded condition, the density distribution in these areas may be better explained by a mechanical bone adaptation theory as compared to the areas that are often less stressed and their density distribution is therefore to a lesser extent driven by mechanical stimulus. Moreover, some movements predict the density distribution better than the other movements (e.g. 80–1001 abduction and 801 flexion, Table 1). Some of these movements (e.g. 8–1001 abduction) are the ones that cause high levels of strain and therefore have important effects on bone adaptation. Since largest density values are found in different areas depending on the movement, one expects the bone adaptation stimulus to be a combination of the stimuli caused by different movements. The linear combination of the different density distributions can combine different loading cases. The weight of every loading case,wi , represents the contribution of that loading case to the overall density distribution. Interestingly, the contribution of eight considered loads is calculated to be zero. The contribution of another loading case (w10 ¼0.07) is close to zero. The most important contributions are therefore from the five remaining load cases. Three of these come from small-angle abduction and flexion movements (w1 ; w2 and w8 ). The small-angle movements are more likely to happen in our daily living as compared to very large-angle movements. However, the loads associated with smallangle abduction and flexion movements are small (Nikooyan et al., 2010). The other important contributions are due to a very largeangle flexion movement (1401, w14 ) and one relatively large-angle abduction movement (1001, w5 ). In vivo measurement of the glenohumeral joint forces has shown that, both in flexion and abduction movements, the joint reaction force (JRF) increases more or less monotonically as the elevation angle increases (Nikooyan et al., 2010). Since JRF is the largest force applied on the scapula, it influences the calculated micro-strains to a large extent. Despite their high loads, large-angle movements occur less frequently in our daily lives than small-angle movements. In short, the movements that contribute the most to the predicted density
Fig. 5. Interrelationship between loading cycles and bone adaptation. Higher levels of micro-strains are needed to induce bone formation, if the number of loading cycles is low. In contrast, loads that generate low micro-strain levels could induce bone formation, if applied at large enough numbers. Adapted by permission from Macmillan Publishers Ltd: Nature Reviews Rheumatology, Ozcivici et al. (2010), copyright 2010.
distribution (Fig. 4d) are either frequent small-angle movements or infrequent large-angle movements. This is in agreement with the general consensus in the field about the relationship between loading cycles and bone adaptation, nicely illustrated by Ozcivici et al. (2010) (Fig. 5). The intermediate loads that do not generate very high levels of micro-strains nor occur very frequently appear to have less contribution to the measured density distribution. Even when movements are combined, there is 15.9% difference between the predicted and measured density distributions (Table 1). This difference can be explained as detailed below. First, we are trying to explain the density distribution only in terms of mechanical loading. It is well known that the bone density distribution also depends on genetic and biochemical factors (Akhter et al., 2000; Robling and Turner, 2002; Ruff et al., 2006). Furthermore, the subject-specific aspects mentioned here mostly refer to bone density and mechanical loading. There are also biological subjectspecific aspects that are not considered in this study. Second, measurement of density distribution using computed tomography is associated with some levels of uncertainty. The accuracy of computed tomography in measurement of density distribution is reported to be up to 10% (Prevrhal et al., 1999) (Table 1). Moreover, the partial volume effect plays a role in the densitometry of the scapula due to the presence of certain thin (see-through) areas within the scapula. When the voxel is too large compared to the dimension of the area within which bone density largely varies, the same voxel includes areas with high and low bone densities. The average density is then too far from both the peak and the valley of actual density values. Despite the relatively high CT resolution (0.6 mm), the voxel size is not small enough to capture the density distribution very accurately in those areas. A voxel size of 100 μm or less would result in more accurate estimation of apparent density values.
2440
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
The results reported here could be compared with those of a recent study (Sharma and Robertson, 2013) in which a 3D bone adaptation simulation of a scapula cadaver specimen using literature-based boundary conditions was performed. The prediction error was found to be less than 0.4 g/cm3 for 53% of the elements of the whole scapula and 66% of the elements in the glenoid area (Sharma and Robertson, 2013). A difference of 0.4 g/cm3 between predicted and measured density values is more than one standard deviation (¼0.150 g/cm3) away from the mean value of the prediction error of our combined loading case (¼0.170 g/cm3). The prediction accuracy of the model presented by Sharma and Robertson (2013) is remarkable considering the fact that they did not use a musculoskeletal model to directly relate musculoskeletal movements to muscle and joint reaction forces experienced by the scapula. The advantage of using literature-based boundary conditions is that much less data is required for running the model. Since the aim of the current study was to relate musculoskeletal movements to bone remodeling results and distinguish between the effects of different types of movements, it was not possible to use the methodology proposed by Sharma and Robertson (2013). An important limitation of the current research is the limited number of subjects studied. It would be interesting to study the prediction accuracy of the proposed method in a large statistical population with different activity patterns (e.g. volleyball players vs. office workers) to see whether the proposed methodology could distinguish between the different types of activities. It should be, however, noted that obtaining the required data for subject-specific modeling of the shoulder joint is difficult and may require extensive image acquisition (both CT and MRI), image processing, and functional testing (Poelert et al., 2013; Zadpoor, in press). Some other parameters of the musculoskeletal model can only be obtained in the dissection room. Therefore, it is a data collection as well as a modeling challenge to do the same type of analysis for a large statistical population. The method proposed in this study has the potential to be used for validation of the loads estimated by musculoskeletal models. In general, two different types of mechanical models are used for estimation of musculoskeletal loads, namely simple mass–spring– damper models (Nikooyan and Zadpoor, 2011a, 2011b, 2012) and large-scale musculoskeletal models (Arnold et al., 2001; Pandy and Anderson, 2000; Pandy et al., 1998; Scheys et al., 2008; Shelburne and Pandy, 1997; van der Helm, 1994b; White et al., 1989). Mass– spring–damper models are simple and require few parameters. However, they cannot provide detailed musculoskeletal loading. Musculoskeletal models need large number of parameters, but are capable of providing the detailed musculoskeletal loading that is needed for FE modeling of bones. However, validation of musculoskeletal models is extremely difficult, because there is no simple way for non-invasive in vivo measurement of musculoskeletal loads. Even though forces measured using instrumented implants (Nikooyan et al., 2010) could be also used for corroboration of musculoskeletal models, instrumented implants cannot measure muscle forces. The method proposed in this study can be used for further corroboration of musculoskeletal models, ideally, when the movements of the subjects are known for a long enough period of time. Much more research is needed regarding the validity of the proposed method and its capability in distinguishing between different types of movements before such an ambitious application can be realized. In conclusion, no single movement could accurately predict the CT-measured density distribution, while a weighted linear combination of the density distributions resulting from several movements could provide a more accurate prediction. The movements that contribute the most to the predicted density distribution are either small-angle (20–401) abduction and flexion movements that generate low forces but occur frequently or very-large angle
flexion (1401) and relatively large-angle (1001) abduction movements that generate large forces but occur infrequently.
Conflict of interest None. References Akhter, M., Iwaniec, U., Covey, M., Cullen, D., Kimmel, D., Recker, R., 2000. Genetic variations in bone density, histomorphometry, and strength in mice. Calcified Tissue International 67, 337–344. Arnold, A.S., Blemker, S.S., Delp, S.L., 2001. Evaluation of a deformable musculoskeletal model for estimating muscle–tendon lengths during crouch gait. Annals of Biomedical Engineering 29, 263–274. Bernd-Arno, B., Ingo, N., Patrick, W., Anas, B., 2009. Numerical investigations on the strain-adaptive bone remodelling in the periprosthetic femur: influence of the boundary conditions. BioMedical Engineering OnLine 8, 1–9. Biewener, A., Fazzalari, N.L., Konieczynski, D., Baudinette, R.V., 1996. Adaptive changes in trabecular architecture in relation to functional strain patterns and disuse. Bone 19, 1–8. Campoli, G., Weinans, H., Zadpoor, A.A., 2012. Computational load estimation of the femur. Journal of the Mechanical Behavior of Biomedical Materials 10, 108–119. Carter, D., Fyhrie, D., Whalen, R., 1987. Trabecular bone density and loading history: regulation of connective tissue biology by mechanical energy. Journal of Biomechanics 20, 785–787. (789–794). Cowin, S., Hegedus, D., 1976. Bone remodeling I: theory of adaptive elasticity. Journal of Elasticity 6, 313–326. Frost, H., 1988. Vital biomechanics: proposed general concepts for skeletal adaptations to mechanical usage. Calcified Tissue International 42, 145–156. Frost, H.M., 1998. A brief review for orthopedic surgeons: fatigue damage (microdamage) in bone (its determinants and clinical implications). Journal of Orthopaedic Science 3, 272–281. Huiskes, R., Weinans, H., Grootenboer, H., Dalstra, M., Fudala, B., Slooff, T., 1987. Adaptive bone-remodeling theory applied to prosthetic-design analysis. Journal of Biomechanics 20, 1135–1150. Jang, I.G., Kim, I.Y., 2008. Computational study of Wolff's law with trabecular architecture in the human proximal femur using topology optimization. Journal of Biomechanics 41, 2353–2361. Lengsfeld, M., G¸nther, D., Pressel, T., Leppek, R., Schmitt, J., Griss, P., 2002. Validation data for periprosthetic bone remodelling theories. Journal of Biomechanics 35, 1553–1564. Morgan, E.F., Bayraktar, H.H., Keaveny, T.M., 2003. Trabecular bone modulus– density relationships depend on anatomic site. Journal of Biomechanics 36, 897–904. Nikooyan, A., Veeger, H., van der Helm, F., Westerhoff, P., Graichen, F., Bergmann, G., 2008. Comparing model-predicted gh-joint contact forces by in-vivo measured forces. Journal of biomechanics 41. (S145). Nikooyan, A., Veeger, H., Westerhoff, P., Graichen, F., Bergmann, G., Van der Helm, F., 2010. Validation of the Delft Shoulder and Elbow Model using in-vivo glenohumeral joint contact forces. Journal of Biomechanics 43, 3007–3014. Nikooyan, A.A., Veeger, H.E.J., Chadwick, E.K.J., Praagman, M., van der Helm, F.C.T., 2011. Development of a comprehensive musculoskeletal model of the shoulder and elbow. Medical and Biological Engineering and Computing 49, 1425–1435. Nikooyan, A.A., Zadpoor, A.A., 2011a. An improved cost function for modeling of muscle activity during running. Journal of Biomechanics 44, 984–987. Nikooyan, A.A., Zadpoor, A.A., 2011b. Mass–spring–damper modelling of the human body to study running and hopping—an overview. Journal of Engineering in Medicine 225, 1121–1135. Nikooyan, A.A., Zadpoor, A.A., 2012. Effects of muscle fatigue on the ground reaction force and soft-tissue vibrations during running: a model study. IEEE Transactions on Biomedical Engineering 59, 797–804. Ozcivici, E., Luu, Y.K., Adler, B., Qin, Y.X., Rubin, J., Judex, S., Rubin, C.T., 2010. Mechanical signals as anabolic agents in bone. Nature Reviews Rheumatology 6, 50–59. Pandy, M.G., Anderson, F.C., 2000. Dynamic simulation of human movement using large-scale models of the body. Phonetica 57, 219–228. Pandy, M.G., Sasaki, K., Kim, S., 1998. A three-dimensional musculoskeletal model of the human knee joint. Part 1: theoretical construct. Computer Methods in Biomechanics and Biomedical Engineering 1, 87. Poelert, S.L., Valstar, E., Weinans, H., Zadpoor, A.A., 2013. Patient-specific finite element modeling of bones. Journal of Engineering in Medicine 227, 464–478. Prevrhal, S., Engelke, K., Kalender, W.A., 1999. Accuracy limits for the determination of cortical width and density: the influence of object size and CT imaging parameters. Physics in Medicine and Biology 44, 751. Robling, A., Turner, C., 2002. Mechanotransduction in bone: genetic effects on mechanosensitivity in mice. Bone 31, 562–569. Ruff, C., Holt, B., Trinkaus, E., 2006. Who's afraid of the big bad Wolff?: “Wolff's law” and bone functional adaptation. American Journal of Physical Anthropology 129, 484–498. Scheys, L., Van Campenhout, A., Spaepen, A., Suetens, P., Jonkers, I., 2008. Personalized MR-based musculoskeletal models compared to rescaled generic
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
models in the presence of increased femoral anteversion:effect on hip moment arm lengths. Gait and Posture 28, 358–365. Sharma, G.B., Debski, R.E., McMahon, P.J., Robertson, D.D., 2009. Adaptive glenoid bone remodeling simulation. Journal of Biomechanics 42, 1460–1468. Sharma, G.B., Debski, R.E., McMahon, P.J., Robertson, D.D., 2010. Effect of glenoid prosthesis design on glenoid bone remodeling: adaptive finite element based simulation. Journal of Biomechanics 43, 1653–1659. Sharma, G.B., Robertson, D.D., 2013. Adaptive scapula bone remodeling computational simulation: relevance to regenerative medicine. Journal of Computational Physics 244, 312–320. Shefelbine, S.J., Augat, P., Claes, L., Simon, U., 2005. Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic. Journal of Biomechanics 38, 2440–2450. Shelburne, K.B., Pandy, M.G., 1997. A musculoskeletal model of the knee for evaluating ligament forces during isometric contractions. Journal of Biomechanics 30, 163–176. ́ Suarez, D.R., Van Der Linden, J.C., Valstar, E.R., Broomans, P., Poort, G., Rozing, P.M., van Keulen, F., 2009. Influence of the positioning of a cementless glenoid prosthesis on its interface micromotions. Proceedings of the Institution of Mechanical Engineer, Part H: Journal of Engineering in Medicines 223, 795–804. ́ Suarez, D.R., Weinans, H., van Keulen, F., 2012. Bone remodelling around a cementless glenoid component. Biomechanics and Modeling in Mechanobiology 11 (6), 903–913. Tanck, E., Ruimerman, R., Huiskes, R., 2006. Trabecular architecture can remain intact for both disuse and overload enhanced resorption characteristics. Journal of Biomechanics 39, 2631–2637. Tsubota, K., Suzuki, Y., Yamada, T., Hojo, M., Makinouchi, A., Adachi, T., 2009. Computer simulation of trabecular remodeling in human proximal femur using large-scale voxel FE models: approach to understanding Wolff's law. Journal of Biomechanics 42, 1088–1094. Turner, A., Gillies, R., Sekel, R., Morris, P., Bruce, W., Walsh, W., 2005. Computational bone remodelling simulations and comparisons with DEXA results. Journal of Orthopaedic Research 23, 705–712.
2441
Turner, C., 1998. Three rules for bone adaptation to mechanical stimuli. Bone 23, 399–407. van der Helm, F.C.T., 1994a. Analysis of the kinematic and dynamic behavior of the shoulder mechanism. Journal of Biomechanics 27, 527–550. van der Helm, F.C.T., 1994b. A finite element musculoskeletal model of the shoulder mechanism. Journal of Biomechanics 27. (551-553–555-569). van der Helm, F.C.T., Pronk, G.M., 1995. Three-dimensional recording and description of motions of the shoulder mechanism. Journal of Biomechanical Engineering 117, 27. van der Helm, F.C.T., Veenbaas, R., 1991. Modelling the mechanical effect of muscles with large attachment sites: application to the shoulder mechanism. Journal of Biomechanics 24, 1151–1163. von Schroeder, H.P., Kuiper, S.D., Botte, M.J., 2001. Osseous anatomy of the scapula. Clinical Orthopedics and Related Research 383, 131–139. Weinans, H., Huiskes, R., Grootenboer, H., 1992. The behavior of adaptive boneremodeling simulation models. Journal of Biomechanics 25, 1425–1441. Weinans, H., Prendergast, P., 1996. Tissue adaptation as a dynamical process far from equilibrium. Bone 19, 143–150. White, S.C., Yack, H.J., Winter, D.A., 1989. A three-dimensional musculoskeletal model for gait analysis. Anatomical variability estimates. Journal of Biomechanics 22, 885–893. Wu, G., van der Helm, F.C.T., Veeger, H.E.J., Makhsous, M., Van Roy, P., Anglin, C., Nagels, J., Karduna, A.R., McQuade, K., Wang, X., Werner, F.W., Buchholz, B., 2005. ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion–Part II: shoulder, elbow, wrist and hand. Journal of Biomechanics 38, 981–992. Zadpoor, A.A., Open forward and inverse problems in theoretical modeling of bone tissue adaptation.Journal of the Mechanical Behavior of Biomedical Materials, http://dx.doi.org/10.1016/j.jmbbm.2013.05.017, in press. Zadpoor, A.A., Campoli, G., Weinans, H., 2013. Neural network prediction of load from the morphology of trabecular bone. Applied Mathematical Modelling 37, 5260–5276.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Subject-specific modeling of the scapula bone tissue adaptation Gianni Campoli a, Harrie Weinans a,b, Frans van der Helm a, Amir A. Zadpoor a,n a Department of Biomechanical Engineering, Faculty of Mechanical, Maritime, and Materials Engineering, Delft University of Technology (TU Delft), Mekelweg 2, Delft 2628 CD, The Netherlands b Department of Orthopaedics, Erasmus University Medical Center, Rotterdam, The Netherlands
art ic l e i nf o
a b s t r a c t
Article history: Accepted 12 July 2013
Adaptation of the scapula bone tissue to mechanical loading is simulated in the current study using a subject-specific three-dimensional finite element model of an intact cadaveric scapula. The loads experienced by the scapula during different types of movements are determined using a subject-specific large-scale musculoskeletal model of the shoulder joint. The obtained density distributions are compared with the CT-measured density distribution of the same scapula. Furthermore, it is assumed that the CTmeasured density distribution can be estimated as a weighted linear combination of the density distributions calculated for different loads experienced during daily life. An optimization algorithm is used to determine the weighting factors of fourteen different loads such that the difference between the weighted linear combination of the calculated density distributions and the CT-measured density is minimal. It is shown that the weighted linear combination of the calculated densities matches the CTmeasured density distribution better than every one of the density distributions calculated for individual movements. The weighting factors of nine out of fourteen loads were estimated to be zero or very close to zero. The five loads that had larger weighting factors were associated with either one of the following categories: (1) small-load small-angle abduction or flexion movements that occur frequently during our daily lives or (2) large-load large-angle abduction or flexion movements that occur infrequently during our daily lives. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Shoulder joint Scapula Bone tissue adaptation Finite element Patient-specific Optimization process
1. Introduction The morphology of bone tissue is partially determined by the mechanical loads it experiences during daily activities (Frost, 1988; Ruff et al., 2006; Turner, 1998). Using adaptation theories (Carter et al., 1987; Cowin and Hegedus, 1976; Huiskes et al., 1987), one could theoretically connect the musculoskeletal loads caused by daily activities to the measured morphology of bone tissue. However, there are two major difficulties in establishing this connection. First, there is no easy way for non-invasive measurement of musculoskeletal loads. Second, bones are loaded differently in the various movements that we carry out in our daily living, and it is not clear how every movement contributes to the measured morphology of bone. These two difficulties exist in the study of both lower- and upper-extremity bones. For example, many researchers have studied bone tissue adaptation of lower extremity bones in both physiological (Campoli et al., 2012; Turner et al., 2005; Weinans et al., 1992; Weinans and Prendergast, 1996) and post-operative conditions (Bernd-Arno et al., 2009; Huiskes et al., 1987; Lengsfeld et al., 2002). The adaptation of the trabecular
n
Corresponding author. Tel.: +31 15 2786794; fax: +31 15 2784717. E-mail addresses: [email protected], [email protected] (A.A. Zadpoor). 0021-9290/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jbiomech.2013.07.024
bone is also extensively studied (Biewener et al., 1996; Jang and Kim, 2008; Shefelbine et al., 2005; Tanck et al., 2006; Tsubota et al., 2009). However, the two above-mentioned difficulties have not been addressed before. In most cases, researchers have simply assumed generic and highly simplified loading conditions. As for upper-extremity bones, even studies that use generic and simplified loading conditions are rare. Indeed, no numerical model of glenoid bone remodeling was available until recently (Sharma et al., 2009, 2010; Sharma and Robertson, 2013). In this study, we present a methodology for establishing the connection between the movements of daily living and the morphology of the scapula and overcoming both above-mentioned difficulties. Once the above-mentioned difficulties are overcome, we could study the effects of individual movements on bone tissue adaptation and separate the effects of different movements from each other. Scapula was used for this study, because distinguishing between different movements is easier in the scapula. That is due to the fact that upper-extremity loading is more varied as compared to the lower-extremity loading that is dominated by gait. Even though currently available models (Sharma and Robertson, 2013) could predict the density distribution of the scapula quite accurately, they do not differentiate between the effects of different types of movements on bone tissue adaptation. We use a validated (Nikooyan et al., 2008, 2010) large-scale musculoskeletal model of the shoulder and elbow, namely the Delft
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
2435
difference between weighted linear summation of the predicted density distributions and the CT-measured density distribution is minimal. The proposed methodology is used for relating the measured density distribution of a scapula to the musculoskeletal loads estimated for the different movements of the same individual using the musculoskeletal model.
Shoulder and Elbow Model (DSEM) (van der Helm, 1994a,b; van der Helm and Pronk, 1995; van der Helm and Veenbaas, 1991). For every movement, this model can provide accurate prediction of detailed musculoskeletal loads including joint reaction and muscle forces and their point of application. In order to overcome the second obstacle, we propose to use an optimization procedure for linking the loads estimated by the musculoskeletal model and their resulting density distributions to the density distribution measured using computed tomography (CT). In this scheme, the CT-measured density distribution is assumed to be explained by a weighted linear combination of the density distributions caused by a number of movements. The optimization process finds the weighting factors of different movements such that the
2. Materials and methods 2.1. Computed tomography One cadaveric scapula of a male donor (57 years) was CT-scanned using a clinical scanner (Siemens, SOMATOMs Definition Flash, construction diameter: 230 mm)
c 44 45
79 58 7 7
78
114
113 52
116
8 47 5 40 51 49
C 1L
8
0 5 69
111
61
126 8 9197257 62 96
53 54 03 1 80
83
JGH
115
lcs TS_ 63 RTS
82
LT
46 112
0 JAC898 789 85 84 6 s 8 lc AA_
G 11072 1
75
106 129 125
98
124
118 74
101
119
122 5 10
100 120
64
109
65
121 6969
73
104
123 110 107
67
89 70 66
131 0 132
108 I lcs RAAI_
Fig. 1. The FE mesh used in simulations (a), bony landmarks and boundary conditions (b) as well as an example of applied joint and muscle loads (c). The inferior–superior and anterior–superior distances of the glenoid fossa, see von Schroeder et al. (2001), were respectively 34 and 27 mm.
2436
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
with an isotropic resolution of 0.6 mm. The CT images were imported into Materialise Mimics and were semi-automatically segmented. Segmentation of the scapula is relatively difficult due to the very thin ‘see-through’ parts of the scapula. After pre-selecting bone areas using a threshold-based technique, a more accurate segmentation was obtained by manually ‘wrapping’ the masked image of the scapula blade cross-section to fill the gaps and voids between disconnected pixels, followed by ‘erosion’ of the pixels that were too far away from the actual image of the thinnest areas of the scapula blade. The segmented scapula was subsequently meshed (Fig. 1a) using 4-node linear tetrahedron elements (edge length¼ 1.1 mm, 349,283 elements).
2.2. Loads and boundary conditions A schematic drawing of the complete methodology is presented in Fig. 2. The details of every step are described below. DSEM was used for estimation of the detailed loading of the scapula while performing two different types of movements: abduction and flexion. These two movements were chosen because of their high frequency in daily living and because many other more complex movements can be related to these two movements. The parameters of the DSEM such as the muscle attachment sites, muscle volume, anatomical parameters and landmarks, shape of the thorax, and several other parameters were directly measured for the same donor.
Kinematic data, movement 1 Kinematic data, movement 2 ...
Musculoskeletal model
Kinematic data, movement n ... Forces, movement 1
Forces, movement 2
Forces, movement n
Finite element model ...
Density distribution, movement 1
Density distribution, movement 2
Density distribution, movement n ...
CT-measured density distribution
w2
wn
New values of weighting factors
w1
Objective function value
No Optimization algorithm
Stopping criteria met?
Stopping criteria of optimization are met, for example, when the objective function does not drop further (within a specified tolerance) for a certain number of iterations.
Yes
End. Fig. 2. A flowchart of the methods used in the current study for establishing the connection between different movements and the morphology of the bone.
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441 These parameters and their definitions can be found in references (Nikooyan et al., 2010, 2011; van der Helm, 1994a,b; van der Helm and Pronk, 1995; van der Helm and Veenbaas, 1991). DSEM is a large-scale musculoskeletal model of the upper extremity and includes thorax, clavicle, scapula, humerus, ulna, and radius. The parameters of the model such as muscle mass, muscle attachments sites, physiological cross-section area of muscle, force-length parameters, optimum fiber lengths, and tendon lengths were measured through an extensive cadaveric study. The model includes 31 muscles and muscle parts and 139 muscle lines of action. Using these directly measured subject-specific parameters, the glenohumeral joint reaction forces as well as muscle forces were estimated with the DSEM for seven different angles of abduction and flexion, namely 201, 401, 601, 801, 1001, 1201 and 1401. The kinematic data needed for solution of the musculoskeletal model were measured in a population of healthy subjects whose anatomies were similar to the donor's anatomy. The average of the measured kinematics data (position of body segments) was used as input to the inverse dynamics musculoskeletal model. By differentiating the vector of position changes with respect to time, one can calculate the velocities and accelerations of different body segments. The mass and inertia properties of the body segments are also known and one could use Newton's equations for calculating internal forces. A physiological muscle load-sharing criterion is used for determining the contribution of every muscle to the calculated forces. The forces estimated with DSEM were applied on the surface of the generated mesh (Fig. 1c). Transferring the musculoskeletal loads to the finite element (FE) mesh required slightly modified points of application of the forces (0.3–2.0 mm distance of muscle forces from the scapula geometry), because the geometry used in the DSEM is a simplified geometry and had slight deviations from the CT-measured geometry. The force and moment equilibrium were re-established after mapping the points of application of musculoskeletal loads onto the FE mesh. The distance between the points of force application given by the DSEM with all the nodes of the FE mesh was calculated. For every load, the node of the mesh that had the least distance from the point of force application was used for applying the force. Similar to other studies (Suaŕ ez et al., 2009, 2012), an optimization algorithm was used to calculate the smallest changes in the magnitude of musculoskeletal forces that will re-establish the equilibrium. Loads were applied as traction surface loads distributed over the area of one single element. The naming convention of bones and bony landmarks that is used here is in agreement with the International Society of Biomechanics (ISB) recommendations (Wu et al., 2005). Three nodes corresponding to three standard bony landmarks of the scapula, namely AA (Anglus Acromial, most posterior), TS (Trigonum Scapula, most medial), and AI (Angulus Inferior, most inferior), were partially or totally constrained. For AA, all degrees of freedom were constrained. The superior–inferior and anterior–posterior displacements (respectively y and z directions in Fig. 1b) were constrained for TS. The displacements of AI were only constrained in the anterior–posterior direction (z direction in Fig. 1b).
2.3. Bone tissue adaptation model The mechanical properties were based on the following correlation between the apparent density and elastic modulus (Morgan et al., 2003): E ¼ aρ ; b
a ¼ 6950; b ¼ 1:49
ð1Þ
Poisson's ratio was fixed at 0:3. The density values were calculated from CTmeasured gray values using a phantom-calibrated linear relationship. A strain-based tissue adaptation model (Campoli et al., 2012; Weinans et al., 1992) was used. In this model, the change in bone density, ρ, is linked to the first principal strain, ε1 : Δρ ¼ B
ε1 k Δt ρ
ð2Þ
The strain threshold to trigger bone remodeling is assumed to be around 1000 micro-strain (Frost, 1998). The maximum density value measured by CT is 2.2 g/cm3. The reference stimulus, k, is calculated as the ratio of the strain threshold to the maximum density value: 1000 10 6/2.2¼0.0005 cm3/g. The parameter B regulates the speed of remodeling. The results of the bone tissue adaptation model were found to be largely independent from the choice of the parameter B. To speed up calculations, a B value of 1000 g2/cm6 s was used in the current study. There are two additional parameters: ρmax and ρmin . Similar to other studies (Sharma et al., 2009, 2010), the density values that are calculated to be large than ρmax were registered as ρmax whereas the ones that were smaller than ρmin were registered as ρmin . These enforced boundaries ensure that the density values remain within the physiological range. One can determine the minimum and maximum density values, ρmax and ρmin , based on the values observed in the CT images (see values in Fig. 4a). An implicit FE solver (ABAQUS) was used. The simulation started from a uniform density value (¼ 0.8 g/cm3) for all elements of the FE model. In every iteration, the first principal strain and the current density of all elements were sent to a subroutine that calculated the change in the density of all elements. The density and, thus, mechanical properties of the elements (Eq. (1)) were accordingly updated (Campoli et al., 2012; Zadpoor et al., 2013). The tissue adaptation model was solved for 25 iterations to ensure that the mean absolute difference between
2437
Fig. 3. The mean of absolute difference between the density distributions of two consecutive iterations divided by the first density distribution plotted vs. iteration number. The results are presented only for the movement (1201 flexion) that had the greatest difference between the density distributions calculated for two first consecutive iterations. two consecutive density distributions divided by the first density distribution was less than 0.1% (Fig. 3). In total, we carried out 14 simulations: one independent simulation for every one of the above-mentioned movements (Fig. 2). 2.4. Linear combination of the predicted density distributions The measured density distribution has resulted from adaptation of the bone tissue in response to a variety of mechanical loads that occur in daily living. We assumed that the measured density distribution is a weighted linear combination of the density distributions that result from the individual loads experienced during daily living: ρCT ¼ ∑i wi ρi
ð3Þ
where ρCT is the CT-measured vector of density distribution, ρi is the vector of the density distribution that results from the load i, and wi is the weight associated with the load i. The weight wi may be a function of the frequency and/or intensity of the load. The weight factors associated with the considered loading cases were determined through a gradient-based nonlinear optimization algorithm with the following objective function: f ðxÞ ¼ jjρCT ∑i wi ρi jj
ð4Þ
The algorithm evaluated the objective function at the end of all iterations (Fig. 2). Even though it is obviously not possible to capture all the loads experienced during our daily activities, one has to include the most important loads that are most likely to influence bone tissue adaptation. The loads that were considered in the current study were the fourteen sets of loads that were estimated by the DSEM for seven different angles of flexion and seven different angles of abduction (see Section 2.2). As previously explained, flexion and abduction are very important movements due to their high frequency of occurrence in daily living and because many other movements can be related to flexion and abduction.
3. Results The FE results showed that no single abduction or flexion movement could accurately explain the CT-measured density distribution (Fig. 4a). In order to make a quantitative comparison between the predicted, ρFE , and measured, ρCT , density distributions, the difference between the two distributions was calculated for all the elements of the FE mesh. For every element j, the prediction error was calculated as ρ ρ j;FE j;CT ð5Þ Errorj ¼ 100 ρj;CT The mean value of the prediction error is, indeed, the mean percent error of the absolute difference between the predicted and CT densities. Two movements, i.e. 801 flexion (Fig. 4b) and 801 abduction (Fig. 4c) are among the movements that best predict the
2438
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
Fig. 4. The CT-measured density distribution (g/cm3) (a) as well as density distribution (g/cm3) calculated for 801 abduction (b), 801 flexion (c), and combined loading case (d) displayed in six different planes. The distance between two consecutive planes is 20 mm.
Table 1 The mean and standard deviation of the absolute and percentage difference between the predicted and CT-measured density distributions for the fourteen considered movements as well as for the linear combination of the fourteen density distributions. Abduction
Mean error (%) SD, error (%) Mean error (g/cm3) SD, error (g/cm3)
Flexion
Combi
201
401
601
801
1001
1201
1401
201
401
601
801
1001
1201
1401
—
36.7 12.3 0.396 0.201
32.7 13.4 0.351 0.195
31.4 13.7 0.337 0.194
29.2 14.0 0.312 0.191
28.8 13.9 0.309 0.188
33.4 13.4 0.361 0.200
36.6 13.0 0.398 0.212
36.9 12.5 0.400 0.206
35.2 13.0 0.381 0.203
33.4 13.4 0.360 0.200
32.6 13.8 0.351 0.202
33.3 13.9 0.360 0.206
34.1 13.6 0.370 0.206
37.4 12.7 0.407 0.212
15.9 13.4 0.170 0.150
Abbreviations: Combi: weighted linear combination of the considered density distributions, SD: standard deviation.
CT-measured density distribution (Table 1). In general, the mean values of the calculated prediction errors are between 29% and 38%.
The weighted linear combination of density distributions better matches the measured density distribution (compare Fig. 4d and a). The results of the applied optimization procedure (Table 2) suggest
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
2439
Table 2 Optimized values of the weighting factors of the fourteen considered movements. These weight factors are the ones for which the minimum value of the objective function is calculated by the optimization algorithm. Abduction
Flexion
201
401
601
801
1001
1201
1401
201
401
601
801
1001
1201
1401
w1 0.21
w2 0.20
w3 0
w4 0
w5 0.48
w6 0
w7 0
w8 0.18
w9 0
w10 0.07
w11 0
w12 0
w13 0
w14 0.38
that the contributions of certain movements are much higher than those of the other movements. In particular, five movements were found to dominate the linear combination due to their much larger weight factors (Table 1, wi 4 0), while nine movements did not contribute much to the calculated density distribution (Table 1, wi ≈0). The mean and median prediction errors were respectively 15.9% and 12.8% for the optimized linear combination of density distributions. The mean prediction error is about half of the smallest mean prediction error calculated for single movements (Table 1). The median of the prediction error (12.8%) is also much smaller than the values calculated for single movements (30.2–38.1%).
4. Discussion It was observed that no single movement could successfully explain all measured density variations. Nevertheless, the density distribution is better predicted in certain areas of the scapula even for single movements. These areas are particularly high-density areas that experience high strain values in various movements. Because of their highly loaded condition, the density distribution in these areas may be better explained by a mechanical bone adaptation theory as compared to the areas that are often less stressed and their density distribution is therefore to a lesser extent driven by mechanical stimulus. Moreover, some movements predict the density distribution better than the other movements (e.g. 80–1001 abduction and 801 flexion, Table 1). Some of these movements (e.g. 8–1001 abduction) are the ones that cause high levels of strain and therefore have important effects on bone adaptation. Since largest density values are found in different areas depending on the movement, one expects the bone adaptation stimulus to be a combination of the stimuli caused by different movements. The linear combination of the different density distributions can combine different loading cases. The weight of every loading case,wi , represents the contribution of that loading case to the overall density distribution. Interestingly, the contribution of eight considered loads is calculated to be zero. The contribution of another loading case (w10 ¼0.07) is close to zero. The most important contributions are therefore from the five remaining load cases. Three of these come from small-angle abduction and flexion movements (w1 ; w2 and w8 ). The small-angle movements are more likely to happen in our daily living as compared to very large-angle movements. However, the loads associated with smallangle abduction and flexion movements are small (Nikooyan et al., 2010). The other important contributions are due to a very largeangle flexion movement (1401, w14 ) and one relatively large-angle abduction movement (1001, w5 ). In vivo measurement of the glenohumeral joint forces has shown that, both in flexion and abduction movements, the joint reaction force (JRF) increases more or less monotonically as the elevation angle increases (Nikooyan et al., 2010). Since JRF is the largest force applied on the scapula, it influences the calculated micro-strains to a large extent. Despite their high loads, large-angle movements occur less frequently in our daily lives than small-angle movements. In short, the movements that contribute the most to the predicted density
Fig. 5. Interrelationship between loading cycles and bone adaptation. Higher levels of micro-strains are needed to induce bone formation, if the number of loading cycles is low. In contrast, loads that generate low micro-strain levels could induce bone formation, if applied at large enough numbers. Adapted by permission from Macmillan Publishers Ltd: Nature Reviews Rheumatology, Ozcivici et al. (2010), copyright 2010.
distribution (Fig. 4d) are either frequent small-angle movements or infrequent large-angle movements. This is in agreement with the general consensus in the field about the relationship between loading cycles and bone adaptation, nicely illustrated by Ozcivici et al. (2010) (Fig. 5). The intermediate loads that do not generate very high levels of micro-strains nor occur very frequently appear to have less contribution to the measured density distribution. Even when movements are combined, there is 15.9% difference between the predicted and measured density distributions (Table 1). This difference can be explained as detailed below. First, we are trying to explain the density distribution only in terms of mechanical loading. It is well known that the bone density distribution also depends on genetic and biochemical factors (Akhter et al., 2000; Robling and Turner, 2002; Ruff et al., 2006). Furthermore, the subject-specific aspects mentioned here mostly refer to bone density and mechanical loading. There are also biological subjectspecific aspects that are not considered in this study. Second, measurement of density distribution using computed tomography is associated with some levels of uncertainty. The accuracy of computed tomography in measurement of density distribution is reported to be up to 10% (Prevrhal et al., 1999) (Table 1). Moreover, the partial volume effect plays a role in the densitometry of the scapula due to the presence of certain thin (see-through) areas within the scapula. When the voxel is too large compared to the dimension of the area within which bone density largely varies, the same voxel includes areas with high and low bone densities. The average density is then too far from both the peak and the valley of actual density values. Despite the relatively high CT resolution (0.6 mm), the voxel size is not small enough to capture the density distribution very accurately in those areas. A voxel size of 100 μm or less would result in more accurate estimation of apparent density values.
2440
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
The results reported here could be compared with those of a recent study (Sharma and Robertson, 2013) in which a 3D bone adaptation simulation of a scapula cadaver specimen using literature-based boundary conditions was performed. The prediction error was found to be less than 0.4 g/cm3 for 53% of the elements of the whole scapula and 66% of the elements in the glenoid area (Sharma and Robertson, 2013). A difference of 0.4 g/cm3 between predicted and measured density values is more than one standard deviation (¼0.150 g/cm3) away from the mean value of the prediction error of our combined loading case (¼0.170 g/cm3). The prediction accuracy of the model presented by Sharma and Robertson (2013) is remarkable considering the fact that they did not use a musculoskeletal model to directly relate musculoskeletal movements to muscle and joint reaction forces experienced by the scapula. The advantage of using literature-based boundary conditions is that much less data is required for running the model. Since the aim of the current study was to relate musculoskeletal movements to bone remodeling results and distinguish between the effects of different types of movements, it was not possible to use the methodology proposed by Sharma and Robertson (2013). An important limitation of the current research is the limited number of subjects studied. It would be interesting to study the prediction accuracy of the proposed method in a large statistical population with different activity patterns (e.g. volleyball players vs. office workers) to see whether the proposed methodology could distinguish between the different types of activities. It should be, however, noted that obtaining the required data for subject-specific modeling of the shoulder joint is difficult and may require extensive image acquisition (both CT and MRI), image processing, and functional testing (Poelert et al., 2013; Zadpoor, in press). Some other parameters of the musculoskeletal model can only be obtained in the dissection room. Therefore, it is a data collection as well as a modeling challenge to do the same type of analysis for a large statistical population. The method proposed in this study has the potential to be used for validation of the loads estimated by musculoskeletal models. In general, two different types of mechanical models are used for estimation of musculoskeletal loads, namely simple mass–spring– damper models (Nikooyan and Zadpoor, 2011a, 2011b, 2012) and large-scale musculoskeletal models (Arnold et al., 2001; Pandy and Anderson, 2000; Pandy et al., 1998; Scheys et al., 2008; Shelburne and Pandy, 1997; van der Helm, 1994b; White et al., 1989). Mass– spring–damper models are simple and require few parameters. However, they cannot provide detailed musculoskeletal loading. Musculoskeletal models need large number of parameters, but are capable of providing the detailed musculoskeletal loading that is needed for FE modeling of bones. However, validation of musculoskeletal models is extremely difficult, because there is no simple way for non-invasive in vivo measurement of musculoskeletal loads. Even though forces measured using instrumented implants (Nikooyan et al., 2010) could be also used for corroboration of musculoskeletal models, instrumented implants cannot measure muscle forces. The method proposed in this study can be used for further corroboration of musculoskeletal models, ideally, when the movements of the subjects are known for a long enough period of time. Much more research is needed regarding the validity of the proposed method and its capability in distinguishing between different types of movements before such an ambitious application can be realized. In conclusion, no single movement could accurately predict the CT-measured density distribution, while a weighted linear combination of the density distributions resulting from several movements could provide a more accurate prediction. The movements that contribute the most to the predicted density distribution are either small-angle (20–401) abduction and flexion movements that generate low forces but occur frequently or very-large angle
flexion (1401) and relatively large-angle (1001) abduction movements that generate large forces but occur infrequently.
Conflict of interest None. References Akhter, M., Iwaniec, U., Covey, M., Cullen, D., Kimmel, D., Recker, R., 2000. Genetic variations in bone density, histomorphometry, and strength in mice. Calcified Tissue International 67, 337–344. Arnold, A.S., Blemker, S.S., Delp, S.L., 2001. Evaluation of a deformable musculoskeletal model for estimating muscle–tendon lengths during crouch gait. Annals of Biomedical Engineering 29, 263–274. Bernd-Arno, B., Ingo, N., Patrick, W., Anas, B., 2009. Numerical investigations on the strain-adaptive bone remodelling in the periprosthetic femur: influence of the boundary conditions. BioMedical Engineering OnLine 8, 1–9. Biewener, A., Fazzalari, N.L., Konieczynski, D., Baudinette, R.V., 1996. Adaptive changes in trabecular architecture in relation to functional strain patterns and disuse. Bone 19, 1–8. Campoli, G., Weinans, H., Zadpoor, A.A., 2012. Computational load estimation of the femur. Journal of the Mechanical Behavior of Biomedical Materials 10, 108–119. Carter, D., Fyhrie, D., Whalen, R., 1987. Trabecular bone density and loading history: regulation of connective tissue biology by mechanical energy. Journal of Biomechanics 20, 785–787. (789–794). Cowin, S., Hegedus, D., 1976. Bone remodeling I: theory of adaptive elasticity. Journal of Elasticity 6, 313–326. Frost, H., 1988. Vital biomechanics: proposed general concepts for skeletal adaptations to mechanical usage. Calcified Tissue International 42, 145–156. Frost, H.M., 1998. A brief review for orthopedic surgeons: fatigue damage (microdamage) in bone (its determinants and clinical implications). Journal of Orthopaedic Science 3, 272–281. Huiskes, R., Weinans, H., Grootenboer, H., Dalstra, M., Fudala, B., Slooff, T., 1987. Adaptive bone-remodeling theory applied to prosthetic-design analysis. Journal of Biomechanics 20, 1135–1150. Jang, I.G., Kim, I.Y., 2008. Computational study of Wolff's law with trabecular architecture in the human proximal femur using topology optimization. Journal of Biomechanics 41, 2353–2361. Lengsfeld, M., G¸nther, D., Pressel, T., Leppek, R., Schmitt, J., Griss, P., 2002. Validation data for periprosthetic bone remodelling theories. Journal of Biomechanics 35, 1553–1564. Morgan, E.F., Bayraktar, H.H., Keaveny, T.M., 2003. Trabecular bone modulus– density relationships depend on anatomic site. Journal of Biomechanics 36, 897–904. Nikooyan, A., Veeger, H., van der Helm, F., Westerhoff, P., Graichen, F., Bergmann, G., 2008. Comparing model-predicted gh-joint contact forces by in-vivo measured forces. Journal of biomechanics 41. (S145). Nikooyan, A., Veeger, H., Westerhoff, P., Graichen, F., Bergmann, G., Van der Helm, F., 2010. Validation of the Delft Shoulder and Elbow Model using in-vivo glenohumeral joint contact forces. Journal of Biomechanics 43, 3007–3014. Nikooyan, A.A., Veeger, H.E.J., Chadwick, E.K.J., Praagman, M., van der Helm, F.C.T., 2011. Development of a comprehensive musculoskeletal model of the shoulder and elbow. Medical and Biological Engineering and Computing 49, 1425–1435. Nikooyan, A.A., Zadpoor, A.A., 2011a. An improved cost function for modeling of muscle activity during running. Journal of Biomechanics 44, 984–987. Nikooyan, A.A., Zadpoor, A.A., 2011b. Mass–spring–damper modelling of the human body to study running and hopping—an overview. Journal of Engineering in Medicine 225, 1121–1135. Nikooyan, A.A., Zadpoor, A.A., 2012. Effects of muscle fatigue on the ground reaction force and soft-tissue vibrations during running: a model study. IEEE Transactions on Biomedical Engineering 59, 797–804. Ozcivici, E., Luu, Y.K., Adler, B., Qin, Y.X., Rubin, J., Judex, S., Rubin, C.T., 2010. Mechanical signals as anabolic agents in bone. Nature Reviews Rheumatology 6, 50–59. Pandy, M.G., Anderson, F.C., 2000. Dynamic simulation of human movement using large-scale models of the body. Phonetica 57, 219–228. Pandy, M.G., Sasaki, K., Kim, S., 1998. A three-dimensional musculoskeletal model of the human knee joint. Part 1: theoretical construct. Computer Methods in Biomechanics and Biomedical Engineering 1, 87. Poelert, S.L., Valstar, E., Weinans, H., Zadpoor, A.A., 2013. Patient-specific finite element modeling of bones. Journal of Engineering in Medicine 227, 464–478. Prevrhal, S., Engelke, K., Kalender, W.A., 1999. Accuracy limits for the determination of cortical width and density: the influence of object size and CT imaging parameters. Physics in Medicine and Biology 44, 751. Robling, A., Turner, C., 2002. Mechanotransduction in bone: genetic effects on mechanosensitivity in mice. Bone 31, 562–569. Ruff, C., Holt, B., Trinkaus, E., 2006. Who's afraid of the big bad Wolff?: “Wolff's law” and bone functional adaptation. American Journal of Physical Anthropology 129, 484–498. Scheys, L., Van Campenhout, A., Spaepen, A., Suetens, P., Jonkers, I., 2008. Personalized MR-based musculoskeletal models compared to rescaled generic
G. Campoli et al. / Journal of Biomechanics 46 (2013) 2434–2441
models in the presence of increased femoral anteversion:effect on hip moment arm lengths. Gait and Posture 28, 358–365. Sharma, G.B., Debski, R.E., McMahon, P.J., Robertson, D.D., 2009. Adaptive glenoid bone remodeling simulation. Journal of Biomechanics 42, 1460–1468. Sharma, G.B., Debski, R.E., McMahon, P.J., Robertson, D.D., 2010. Effect of glenoid prosthesis design on glenoid bone remodeling: adaptive finite element based simulation. Journal of Biomechanics 43, 1653–1659. Sharma, G.B., Robertson, D.D., 2013. Adaptive scapula bone remodeling computational simulation: relevance to regenerative medicine. Journal of Computational Physics 244, 312–320. Shefelbine, S.J., Augat, P., Claes, L., Simon, U., 2005. Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic. Journal of Biomechanics 38, 2440–2450. Shelburne, K.B., Pandy, M.G., 1997. A musculoskeletal model of the knee for evaluating ligament forces during isometric contractions. Journal of Biomechanics 30, 163–176. ́ Suarez, D.R., Van Der Linden, J.C., Valstar, E.R., Broomans, P., Poort, G., Rozing, P.M., van Keulen, F., 2009. Influence of the positioning of a cementless glenoid prosthesis on its interface micromotions. Proceedings of the Institution of Mechanical Engineer, Part H: Journal of Engineering in Medicines 223, 795–804. ́ Suarez, D.R., Weinans, H., van Keulen, F., 2012. Bone remodelling around a cementless glenoid component. Biomechanics and Modeling in Mechanobiology 11 (6), 903–913. Tanck, E., Ruimerman, R., Huiskes, R., 2006. Trabecular architecture can remain intact for both disuse and overload enhanced resorption characteristics. Journal of Biomechanics 39, 2631–2637. Tsubota, K., Suzuki, Y., Yamada, T., Hojo, M., Makinouchi, A., Adachi, T., 2009. Computer simulation of trabecular remodeling in human proximal femur using large-scale voxel FE models: approach to understanding Wolff's law. Journal of Biomechanics 42, 1088–1094. Turner, A., Gillies, R., Sekel, R., Morris, P., Bruce, W., Walsh, W., 2005. Computational bone remodelling simulations and comparisons with DEXA results. Journal of Orthopaedic Research 23, 705–712.
2441
Turner, C., 1998. Three rules for bone adaptation to mechanical stimuli. Bone 23, 399–407. van der Helm, F.C.T., 1994a. Analysis of the kinematic and dynamic behavior of the shoulder mechanism. Journal of Biomechanics 27, 527–550. van der Helm, F.C.T., 1994b. A finite element musculoskeletal model of the shoulder mechanism. Journal of Biomechanics 27. (551-553–555-569). van der Helm, F.C.T., Pronk, G.M., 1995. Three-dimensional recording and description of motions of the shoulder mechanism. Journal of Biomechanical Engineering 117, 27. van der Helm, F.C.T., Veenbaas, R., 1991. Modelling the mechanical effect of muscles with large attachment sites: application to the shoulder mechanism. Journal of Biomechanics 24, 1151–1163. von Schroeder, H.P., Kuiper, S.D., Botte, M.J., 2001. Osseous anatomy of the scapula. Clinical Orthopedics and Related Research 383, 131–139. Weinans, H., Huiskes, R., Grootenboer, H., 1992. The behavior of adaptive boneremodeling simulation models. Journal of Biomechanics 25, 1425–1441. Weinans, H., Prendergast, P., 1996. Tissue adaptation as a dynamical process far from equilibrium. Bone 19, 143–150. White, S.C., Yack, H.J., Winter, D.A., 1989. A three-dimensional musculoskeletal model for gait analysis. Anatomical variability estimates. Journal of Biomechanics 22, 885–893. Wu, G., van der Helm, F.C.T., Veeger, H.E.J., Makhsous, M., Van Roy, P., Anglin, C., Nagels, J., Karduna, A.R., McQuade, K., Wang, X., Werner, F.W., Buchholz, B., 2005. ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion–Part II: shoulder, elbow, wrist and hand. Journal of Biomechanics 38, 981–992. Zadpoor, A.A., Open forward and inverse problems in theoretical modeling of bone tissue adaptation.Journal of the Mechanical Behavior of Biomedical Materials, http://dx.doi.org/10.1016/j.jmbbm.2013.05.017, in press. Zadpoor, A.A., Campoli, G., Weinans, H., 2013. Neural network prediction of load from the morphology of trabecular bone. Applied Mathematical Modelling 37, 5260–5276.