- Email: [email protected]
Computational simulation of trabecular adaptation progress in human proximal femur during growth
ARTICLE IN PRESS Journal of Biomechanics 42 (2009) 573–580
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Computational simulation of trabecular adaptation progress in human proximal femur during growth In Gwun Jang a, Il Yong Kim b, a b
Department of Mechanical and Materials Engineering, Jackson Hall 213, 35 Fifth Field Company Lane, Queen’s University, Kingston, Ontario, Canada K7L 3N6 Department of Mechanical and Materials Engineering, McLaughlin Hall 305, 130 Stuart Street, Queen’s University, Kingston, Ontario, Canada K7L 3N6
a r t i c l e in fo
abstract
Article history: Accepted 18 December 2008
There are a large number of clinical and experimental studies that analyzed trabecular architecture as a result of bone adaptation. However, only a limited amount of quantitative data is currently available on the progress of trabecular adaptation during growth. In this paper, we proposed a two-step numerical simulation method that predicts trabecular adaptation progress during growth using a recently developed topology optimization algorithm, design space optimization (DSO), under the hypothesis that the mechanisms of DSO are functionally equivalent to those of bone adaptation. We applied the proposed scheme to trabecular adaptation simulation in human proximal femur. For the simulation, the full trabecular architecture in human proximal femur was represented by a two-dimensional mFE model with 50 mm resolution. In Step 1, we determined a reference value that regulates trabecular adaptation in human proximal femur. In Step 2, we simulated trabecular adaptation in human proximal femur during growth with the reference value derived in Step 1. We analyzed the architectural and mechanical properties of trabecular patterns through iterations. From the comparison with experimental data in the literature, we showed that in the early growth stage trabecular adaptation was achieved mainly by increasing bone volume fraction (or trabecular thickness), while in the later stage of the development the trabecular architecture gained higher structural efficiency by increasing structural anisotropy with a relatively low level of bone volume fraction (or trabecular thickness). We demonstrated that the proposed numerical framework predicted the growing progress of trabecular bone that has a close correlation with experimental data. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Trabecular adaptation Bone growth Human proximal femur Topology optimization Design space optimization
1. Introduction Although there are a number of clinical and experimental studies on trabecular adaptation, only few quantitative data are currently available on the architectural changes in trabecular adaptation during growth. Korstjens et al. (1995) found that fine trabecular patterns in young children develop into coarser ones as they age. Nafei et al. (2000) studied the change of morphological properties in growing trabecular ovine bone. Based on the experiment with five different age groups of female pigs, Tanck et al. (2001) showed that in the early phase of growth the bone density changes significantly, whereas in the later stage of the development mainly the trabecular architecture undergoes adaptation with a small change in overall density. In an attempt to numerically simulate trabecular adaptation progress, phenomenological update algorithms were used with the aid of the finite element method; these studies sought to
Corresponding author. Tel.: +1 613 533 3077; fax: +1 613 533 6489.
E-mail addresses: [email protected] (I.G. Jang), [email protected] (I.Y. Kim). 0021-9290/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2008.12.009
quantitatively describe the mechanical stimulus (cause) and the bone response (effect). In the phenomenological models, a variety of mechanical stimuli has been utilized: stress (Adachi et al., 2001; Pauwels, 1965), strain (Carter et al., 1981; Cowin and Hegedus, 1976), strain-energy density (Fyhrie and Carter, 1986; Huiskes et al., 1987; Weinans et al., 1992), and strain-energy density rate (Huiskes et al., 2000; Ruimerman et al., 2005). The aforementioned mechanical stimuli used in the different models have been shown to be consistent with one another in boneremodeling simulation. Recently developed phenomenological models with so-called mFE models enable us to simulate bone remodeling or adaptation at the tissue level even in threedimensions (Ruimerman et al., 2005; Tsubota et al., 2003). In these studies, mFE models were built to represent trabecular architecture with the resolution of tens of mm, and in situ mechanical stimuli were used for the regulation of the trabecular adaptation. As a different numerical approach, topology optimization methods (Bendsøe, 1989; Bendsøe and Kikuchi, 1988) have been successfully applied and produced a variety of meaningful results in bone-remodeling simulation (Bagge, 2000; Fernandes and Rodrigues, 1999; Hollister et al., 1994). Recently, the authors developed a new topology optimization scheme that implements
ARTICLE IN PRESS 574
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
an extra perimeter constraint for obtaining trabecular architecture in human proximal femur (Jang and Kim, 2008). The numerical results showed a close correlation with clinical data in the literature. It is interesting to note that phenomenological methods share the common underlying assumption with topology optimizationbased methods, in that the aim of bone remodeling is to use the bone material in an ‘‘optimal’’ way. Their similarities can also be revealed from the view point of fully stressed design formulation, where each member of the structure is fully stressed under the multiple loading conditions for the optimum design (Haftka and Gu¨rdal, 1992). Jang and Kim (2008) recently showed that the trabecular architecture in human proximal femur is constructed as a fully stressed structure to support multi-directional and multi-axial load cases in daily activities. Huiskes (2000) summarized this observation: ‘‘It has been shown that computer algorithms for material design optimization, when given the opportunity, show the tendency to produce trabecular patterns, similar to bone.’’ In addition to this qualitative comparison, the phenomenological and topology optimization-based approaches have been explicitly compared (Harrigan and Hamilton, 1994; Subbarayan and Bartel, 2000; Tovar, 2004), and it has been shown that their structural analogy is remarkable (Nowak, 2006). Especially Jang et al. (2009) derived the mathematical relations among the parameters in topology optimization and phenomenological approach. However, as Hart (2001) pointed out, a critical deficiency of topology optimization is that the numerical approach fails to obtain physiological ‘‘intermediate’’ adaptation process and is, therefore, unable to provide an understanding of the ‘‘entire’’ bone-adaptation process. The aforementioned limitation of topology optimization stems from the fact that the design domain does not change during the optimization process; actual trabecular adaptation involves the sequential change of the trabecular surface via osteoclast resorption and osteoblast deposition (Parfitt, 1994). Due to this limitation of conventional topology optimization, we could obtain only ‘‘final’’ trabecular architecture without any information on ‘‘intermediate’’ trabecular adaptation in our previous paper (Jang and Kim, 2008). Compared to the fixed design domain in conventional topology optimization, a recently developed design space optimization (DSO) incorporates an evolutionary feature in order to adaptively change the design domain (Jang and Kwak, 2006; Kim and Kwak, 2002). The objectives of this paper were to numerically predict the entire trabecular adaptation progress in human proximal femur during growth by using DSO for the first time in the field, and to compare the numerical results to the experimental data in the literature. We developed a two-step simulation procedure: determining a reference value for trabecular adaptation (Step 1) and then simulating trabecular adaptation (Step 2). We showed that in the early growth stage where the trabecular architecture was less optimal, trabecular adaptation was achieved mainly by increasing bone volume fraction (or trabecular thickness), whereas in later stage of the development the trabecular architecture achieved higher structural efficiency by increasing structural anisotropy with a relatively low level of bone volume fraction (or trabecular thickness).
2. Methods 2.1. Design space optimization and its application to trabecular adaptation Design space optimization consists of double loops: the inner loop which is the same as the conventional topology optimization and the outer loop which controls design domain change, as shown in Fig. 1. When an optimization process in the inner loop with the current design domain reaches convergence status, the outer loop seeks out a better design domain by expanding the design space (at the start of the inner loop) and by shrinking the design space (at the end of the inner loop).
Fig. 1. Flow chart of design space optimization (DSO).
Fig. 2. Schematic diagram for trabecular adaptation and design space optimization. This evolutionary nature of DSO allowed us to hypothesize that the mechanisms of DSO are functionally equivalent to those of the trabecular adaptation, as illustrated in Fig. 2. Osteoblast bone formation in trabecular adaptation may be considered to be equivalent to design space expansion in DSO, in which process one layer of new elements with very low density is added along the boundary of the current whole design domain. Osteoclast bone resorption is conceptually similar with design space reduction in the sense that inefficiently used elements are deleted from the current design domain. In this paper, starting from initial trabecular bone (Fig. 3(a)), we added new elements with the lowest density (0.05 in this paper) on all over the boundary of the bone structure (Fig. 3(b)). This means we did not use any specific criterion—such as stress, strain, or strain-energy density requirement—to add elements. After expanding the trabecular bone domain this way, we determined new, balanced morphology using ‘‘optimization’’. Then, inefficiently used elements which have bone densities below 0.1 were removed at the end of inner loop (Fig. 3(c)). The steps shown in Fig. 3(a)–(c) constitute a single iteration in our numerical method. As this numerical routine progressed, trabecular morphology gradually evolved into a more optimized one. 2.2. Two-step simulation for trabecular adaptation 2.2.1. Step 1 The main objective of Step 1 was to determine the reference value k used in the bone-remodeling update equation (Huiskes et al., 2000; Ruimerman et al., 2005;
ARTICLE IN PRESS I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
575
Fig. 3. Application of design space optimization to bone adaptation. (a) Initial trabecular morphology; (b) expanded trabecular morphology right after design space expansion; and (c) resultant trabecular morphology right after design space reduction. element, h1 and h2 are lengths of a finite element (both 50 mm in this paper), n1 and n2 are divisions of the domain along x1 and x2 directions, and ri,j represents the bone density located in the coordinate of (i, j) in the mesh (Petersson, 1999). The mass and perimeter constants in Eq. (2), M0 and P0, were calculated from the initial trabecular structure. Sensitivity equations for Eq. (2) are presented in our previous work (Jang and Kim, 2008). Because the reference value k for bone remodeling in human proximal femur is not known, we approximated the value using the simulation result by Eq. (2) as PN k ffi kapprox ¼
rÞ
i¼1 ðU a;i =
N
,
(3)
which is the average value of Ua,i/r over the trabecular architecture determined at the end of Step 1. It was shown in our previous work that Eq. (2) produced physiological results (Jang and Kim, 2008), and therefore kapprox obtained by Eq. (3) may be considered as an effective estimate for subsequent trabecular adaptation simulation, Step 2. 2.2.2. Step 2 We reformulated Eq. (2) to simulate trabecular adaptation during growth as 3 N X X 1 T cj r i vi uj Kuj þ a 2 j¼1 i¼1
Minimize f ðqÞ ¼ Fig. 4. FE model for an initial trabecular architecture in human proximal femur.
Subject to gðqÞ ¼ P ¼
n2 X n1 X
h2 jri;j ri1;j j
j¼0 i¼1
Weinans et al., 1992)
þ
U a;i dri ¼B k dt ri
Subject to
f ðqÞ ¼
h1 jri;j ri;j1 jXP 0
j¼1 i¼0
(1)
where ri represents the bone density of the ith element, B the constant, and Ua,i the apparent strain-energy density of the ith element. First, the trabecular domain was discretized into finite elements with the size of 50 mm 50 mm. In order to represent initial trabecular architecture, we randomly generated hollow circular patterns with an external diameter of 1400 mm and an internal diameter of 1200 mm (Fig. 4). The circular pattern was chosen to construct initial trabecular architecture with thin and near-isotropic trabecular patterns based on the observation by Nafei et al. (2000). The total number of finite elements needed for the definition of the initial trabecular architecture was 0.75 million. Check more detailed information on FE modeling and analysis in Appendix A. Then, total strain energy was to be minimized under three load cases subject to mass (M) and perimeter (P) constraints. The design variables were the bone densities (r) in the FE elements with the total number N. The bone density vector, r, is composed of each bone density ri of the ith element as q ¼ ½ r1 r2 . . . ri . . . rN T . The perimeter of a structure is the sum of the lengths of all inner and outer boundaries. In this study, a minimum allowable perimeter constraint was applied to effectively represent porosity of the trabecular bone and produce physiological numerical results (Jang and Kim, 2008). Using the information, we formulated the optimization problem as Minimize
n2 X n1 X
3 X 1 T uj Kuj cj 2 j¼1
g 2 ðqÞ ¼ P ¼
N X
i vi pM 0 i¼1 n n 2 1 XX
g 1 ðqÞ ¼ M ¼
r
n2 X n1 X
h1 jri;j ri;j1 jXP 0
j¼1 i¼0
0:05pri p1:0
(4)
a ¼ gk ffi gkapprox .
(5)
The reference value k with Eq. (5) represents the balance between two P3 T competing measures in Eq. (4): total strain energy j¼1 c j ðð1=2Þuj Kuj Þ and total PN bone mass i¼1 ri vi . As k increases, more emphasis is placed on the minimization of the second term (total bone mass), and the optimum solution would have a relatively lower total bone mass and higher total strain energy. It can also be seen from Eq. (1) that a large k will slow down the increase in bone density (accordingly bone mass) over the iterations. Sensitivity equations for Eq. (4) were derived as 3 X @f ðqÞ 1 T @K uj ¼ cj uj þ ave @re 2 @re j¼1
h2 jri;j ri1;j j
j¼0 i¼1
þ
0:05pri p1:0.
where a is to be determined on the basis of kapprox obtained in Step 1. Harrigan and Hamilton (1994) mathematically verified that a global function in the form of Eq. (4) is exactly equivalent to a strain-energy density-based boneremodeling equation, Eq. (1). In terms of bone-mass change during simulation, when using Eq. (2), we have ‘‘constant’’ total bone mass during simulation, which is not physiological in trabecular adaption during growth. With Eq. (4), however, we can vary the total bone mass during simulation while balancing the bone mass and total strain energy so as to simulate physiological trabecular adaptation progress. In our previous study (Jang et al., 2009), we mathematically derived the relationships among the reference value k in Eq. (1), the penalty parameter a in Eq. (4), and the exponent g in the relationship between bone density r and Young’s Modulus E of bone material (E ¼ E0rg) as (check more detailed mathematical derivation for Eq. (5) in Appendix B)
(2)
where cj is the normalized weighting factor for three load cases (c1 ¼ 0.6 for onelegged stance, c2 ¼ 0.2 for abduction, and c3 ¼ 0.2 for adduction), uj is the nodal displacement vector of the jth load case, K stiffness matrix, vi the volume of the ith
n2 X n1 X @ri;j @gðqÞ ¼ h2 sgnðri;j ri1;j Þ @re @re j¼0 i¼1
þ
n2 X n1 X j¼1 i¼0
h1
@ri;j @re
sgnðri;j ri;j1 Þ.
(6)
The overall procedure for this two-step simulation is briefly illustrated in Fig. 5.
ARTICLE IN PRESS 576
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Fig. 5. Schematic diagram for two-step simulation for trabecular adaptation during growth.
Trabecular adaptation in Step 2 resulted in the changes of morphological properties of trabecular patterns: bone volume fraction ratio (BV/TV), trabecular thickness (Tb.Th), structural anisotropy, and strain-energy density. In order to observe the histories of the four indices during growth, four volume of interest (VOIs) with the size of 10 mm 10 mm were selected from the epiphysis of proximal femur, the metaphysis of proximal femur, Ward’s triangle, and intertrochanteric trabecular group. The structural anisotropy was determined from the mean intercept length (MILmax/MILmin) (Whitehouse and Dyson, 1974).
3. Results In Step 1, we obtained kapprox ¼ 1.54 102 J g1, and this value was used to calculate a in Eq. (4) in Step 2. If kapprox obtained from the final trabecular architecture in Step 2 is close to the original kapprox in Step 1, we can verify that the original kapprox was an effective estimate for trabecular adaptation during growth in human proximal femur. Indeed, the two ks should be ideally the same. However, they may have some differences due to numerical issues including round-off errors and error accumulation during simulation. When we approximate the reference value k from the result in Step 1 and put it into Eq. (4) with the form of a in Step 2, it implicitly has the meaning that in Step 2 we are expected to obtain bone architecture with similar bone mass and the reference value k calculated in Step 1, because the derived value a with Eq. (5) represents the regulating measure for the total strain energy and bone mass in trabecular architecture in Eq. (4). We obtained kapprox of 1.56 102 J g1 in Step 2; the difference was only 1.3%. In Step 2, the initial trabecular structure in a human proximal femur adapted to its surrounding mechanical environment through 20 iterations as shown in Fig. 6; the final result presented in Fig. 6(d) is very similar to the actual patterns observed in the
literature as shown in Fig. 7 (Gray, 1918; Skedros and Baucom, 2007). Fig. 8 shows the complete histories of the four characteristics during iterations. The bone volume fraction increased rapidly in the early growth phase, reached its maximum at around iteration 5, decreased slightly, and stabilized from iteration 10 for all VOIs. Trabecular thickness in Fig. 8(b) showed a similar trend, which indicates a strong correlation between bone volume fraction and trabecular thickness. On average the early structural anisotropy was relatively low, and there was only a small change until iteration 5. After this point, anisotropy increased gradually until iteration 15, and then it flattened out (Fig. 8(c)). In Fig. 8(d), we investigated the history of strain-energy density, which is a measure of the efficiency of the structural architecture. Strain-energy density in all VOIs reached their minimum values at iteration 5 and then stayed at the same level.
4. Discussion The trabecular adaptation process may be explained from the view point of bone stiffness. Bone stiffness is determined by three parameters: tissue modulus, bone volume fraction, and bone architecture (Kabel et al., 1999; van Rietbergen et al., 1998). In this study, the tissue modulus was assumed to be constant, and therefore the bone stiffness developed only through the combined effect of bone volume fraction and anisotropy. This suggests that, in the early growth stage (iterations 1–5) when trabecular adaptation was achieved mainly by increasing bone volume fraction (or trabecular thickness) (Fig. 8(a) and (b)), the trabecular architecture was less optimal; on the other hand the trabecular architecture gained higher structural efficiency in the later development stage (iterations 5–15) by increasing structural
ARTICLE IN PRESS I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Fig. 6. Trabecular architecture in human proximal femur during growth.
577
anisotropy (Fig. 8(c)) even with a relatively low level of bone volume fraction (or thickness). General trends of the architectural and mechanical properties during growth in our simulation (Fig. 8) can be compared to the experimental data in the literature (Kneissel et al., 1997; Nafei et al., 2000; Tanck et al., 2001). In Fig. 9, bone volume fraction and anisotropy in the vertebra and proximal tibia of female pigs during 6–230 weeks of age are presented (Tanck et al., 2001). Examining the growth plate in the proximal tibia, Tanck et al. (2001) estimated that 104 weeks in pigs is roughly equivalent to late adolescence (about 15 years) for humans. Their experimental data in Fig. 9(a) shows a close correlation with our numerical data in Fig. 8(a), and therefore we estimated 1 iteration in this study to be 12.5 weeks for pigs or 1.8 years for human. In this study, bone volume fraction and trabecular thickness reached their maximum values at around iteration 5 (approximately 9 years for human), and a similar trend was observed in medieval Nubian skeletons by Kneissel et al. (1997): the highest bone volume fraction and trabecular thickness occurred during adolescence. Nafei et al. (2000) also found that architectural properties of trabecular bone in sheep change significantly with skeletal maturity. The bone volume fraction of the skeletally mature group was 22% higher than that of the skeletally immature group, and the structural anisotropy was 25% higher in the skeletally mature than in the skeletally immature sheep. As Huiskes et al. (2000) stated, homeostasis in a living organism is the property of a control system which regulates its internal parameters so as to maintain a stable, constant condition over possibly wide environmental variations. According to the classical control theory (Franklin et al., 2002), a general timeresponse of the system to any input or disturbance consists of rise time, peak overshoot, and others such as settling time, quarterdecay, etc. It is interesting to observe that trabecular adaptation history obtained in this paper (Fig. 8) shows common features of a control system, and it suggests that the proposed method plays the role of a regulator for bone adaptation. There are very few numerical studies that determined the trabecular adaptation progress of full trabecular architecture, because of computational difficulties. Specifically, with Euler forward integration in phenomenological approaches, the numerical stability theory typically requires that information does not propagate across more than one element per time step (Cook et al., 1989). Hence, the time step for a mFE model must be very small, and usually hundreds of iterations are needed whose computing time is very long (Ruimerman et al., 2005). On the other hand, topology optimization-based methods show better performance in convergence and computing time despite of their non-physiological intermediate progresses (Jang et al., 2009). In this paper, we successfully predicted physiological intermediate trabecular adaptation progress by virtue of evolutionary nature of DSO while improving the overall computational efficiency. We should address the limitation and challenge of this research. Due to a very limited amount of clinical data, it was very difficult to directly compare our numerical results to the literature data in terms of species (human in Kneissel et al.’s paper; sheep in Nafei et al.’s paper; pig in Tanck et al.’s paper), age span (0–60 years in Kneissel et al.’s paper; 3–80 months in Nafei et al.’s paper; 6–230 weeks in Tanck et al.’s paper), and the anatomical site of sampling (vertebra in Kneissel et al.’s paper; proximal tibia in Nafei et al.’s paper; proximal tibia and vertebra in Tanck et al.’s paper). Nevertheless, the common feature on these different samples is that they are major weight-bearing bones in the skeletal systems. The sites of sampling were also chosen from parts which are loaded predominantly in axial compression.
ARTICLE IN PRESS 578
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Fig. 7. Photograph and radiograph of human proximal femur (by courtesy of John G. Skedros; Skedros and Baucom, 2007) (Reprinted with permission from Elsevier).
Fig. 8. Architectural and mechanical properties in human proximal femur during growth. (a) Bone volume fraction (BV/TV); (b) Trabecular thickness (Tb. Th, mm); (c) Structural anisotropy; and (d) Strain energy density (J cm3).
Another limitation is two-dimensional aspects of geometry and loads in the FE model. Our previous paper (Jang and Kim, 2008) showed that final trabecular architecture from topology optimization agreed well with the real trabecular architecture even with two-dimensional models. The good correlation in terms of Ward’s classification (Gray, 1918; Ward, 1838; Whitehouse and Dyson, 1974) and intersection angles of trabecular groups is mainly because the hip joint contact force in the analysis is the greatest load on the mediolateral plane of a proximal femur, and moreover its magnitude is at least an order of magnitude higher
than other loads. Even though we obtained trabecular architecture which is structurally similar to photograph and radiograph images of human proximal femur, no quantitative comparisons could be made between the numerically predicted values of bone volume fraction, trabecular thickness, or structural anisotropy and those measured on bone specimens experimentally. Another consequence of using a two-dimensional model is that the mechanical properties would be more sensitive to reductions in the thickness and number of trabeculae than would be expected for a threedimensional model of trabecular bone, because two-dimensional
ARTICLE IN PRESS I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
579
Acknowledgements The authors would like to thank Dr. Krister Svanberg for providing the MMA code and Dr. John G. Skedros for providing the photograph and radiograph images of human proximal femur. The authors would also like to thank Ryan Willing at Queen’s University (Kingston, Canada) for helpful discussion.
Appendix. Supporting Information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.jbiomech.2008.12.009.
References
Fig. 9. History of trabecular architecture of female pig during growth (Tanck et al., 2001) (Reprinted with permission from Elsevier). (a) Bone volume fraction of trabecular architecture in a female pig through age and (b) anisotropy of trabecular architecture in a female pig through age.
networks have fewer connections, on average, at each node than three-dimensional networks (Silva and Gibson, 1997). The proposed method in this paper is mainly based on topology optimization and is enhanced by adding the evolutionary feature of design space optimization to mimic the bone formation and resorption activities in trabecular adaptation. We performed two-step simulation to obtain the intermediate trabecular architecture throughout trabecular adaptation as well as its final architecture in human proximal femur during growth. From the comparison with clinical and experimental data, it is shown that the proposed numerical framework predicted physiological trabecular adaptation progress during growth. It means that, as other researchers already found, topology optimizationbased methods are compatible with the phenomenological methods in bone-remodeling simulation. Moreover, the proposed method can be used to predict the progress of trabecular adaptation unlike other conventional topology optimizationbased approaches. The trends of the architectural and mechanical properties obtained by the simulation would improve the insight into trabecular adaptation during growth. However, with the twodimensional aspects of the model in this paper, we could not quantitatively compare morphometrical indices from our simulation with those from experiments. For thorough and exact comparison, it is required to perform trabecular adaptation simulation with the corresponding three-dimensional model, which would be our future work.
Conflict of interest It is stated that the authors have no conflict of interest including any financial and personal relationships with other people or organizations that could inappropriately influence their work.
Adachi, T., Tsubota, K., Tomita, Y., Hollister, S.J., 2001. Trabecular surface remodeling simulation for cancellous bone using microstructural voxel finite element models. Journal of Biomechanical Engineering 123, 403–409. Bagge, M., 2000. A model of bone adaptation as an optimization process. Journal of Biomechanics 33, 1349–1357. Bendsøe, M.P., 1989. Optimal shape design as a material distribution problem. Structural Optimization 1, 193–303. Bendsøe, M.P., Kikuchi, N., 1988. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 71, 197–224. Carter, D.R., Harris, W.H., Vasu, R., Caler, W.E., 1981. The mechanical and biological response of cortical bone to in vivo strain histories. In: Cowin, S.C. (Ed.), Mechanical Properties of Bone, American Society of Mechanical Engineers. Cook, R.D., Malkus, D.S., Plesha, M.E., 1989. Concepts and Applications of Finite Element Analysis, third ed. Wiley, New York. Cowin, S.C., Hegedus, D.H., 1976. Bone remodeling I: a theory of adaptive elasticity. Journal of Elasticity 6, 313–326. Fernandes, P.R., Rodrigues, H.C., 1999. A material optimization model for bone remodeling around cementless hip stems. In: Proceedings of the Ninth European Conference on Computational Mechanics, Munich, Germany. Franklin, G.F., Powell, J.D., Emami-Naeini, A., 2002. Feedback Control of Dynamic Systems. Prentice-Hall, New Jersey. Fyhrie, D.P., Carter, D.R., 1986. A unifying principle relating stress to trabecular bone morphology. Journal of Orthopaedic Research 4, 304–317. Gray, H., 1918. Anatomy of the human body. Philadelphia, Lea & Febiger. (Revised and re-edited by Lewis W.H., 2000, twentieth ed., Bartleby.com, New York. Haftka, R.T., Gu¨rdal, Z., 1992. Elements of Structural Optimization, third ed. Kluwer Academic Publishers, Dordrecht. Harrigan, T.P., Hamilton, J.J., 1994. Bone remodeling and structural optimization. Journal of Biomechanics 27, 323–328. Hart, R.T., 2001. Bone modeling and remodeling: theories and computation. In: Cowin, S.C. (Ed.), Bone Mechanics Handbook. CRC Press, Boca Raton, FL. Hollister, S.J., Brennan, J.M., Kikuchi, N., 1994. A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress. Journal of Biomechanics 27, 433–444. Huiskes, R., 2000. If bone is the answer, then what is the question? Journal of Anatomy 197, 145–156. Huiskes, R., Ruimerman, R., van Lenthe, G.H., Janssen, J.D., 2000. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 405, 704–706. Huiskes, R., Weinans, H., Grootenboer, J., Dalstra, M., Fudala, M., Slooff, T.J., 1987. Adaptive bone remodelling 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. Jang, I.G., Kim, I.Y., Kwak, B.M., 2009. Analogy of strain energy density based boneremodeling algorithm and structural topology optimization. Journal of Biomechanical Engineering 131 (1), 011012. Jang, I.G., Kwak, B.M., 2006. Evolutionary topology optimization using design space adjustment based on fixed grid. International Journal for Numerical Methods in Engineering 66 (11), 1817–1840. Kabel, J., Odgaard, A., van Rietberge, B., Huiskes, R., 1999. Connectivity and elastic properties of cancellous bone. Bone 24, 115–120. Kim, I.Y., Kwak, B.M., 2002. Design space optimization using a numerical design continuation method. International Journal for Numerical Methods in Engineering 53, 1979–2002. Kneissel, M., Roscher, P., Steiner, W., Schamall, D., Kalchhauser, G., Boyde, A., TeschlerNicola, M., 1997. Cancellous bone structure in the growing and aging lumbar spine in a historic Nubian population. Calcified Tissue International 61, 95–100. Korstjens, C.M., Geraets, W.G.M., van Ginkel, F.C., Prahl-Andersen, B., van Der Stelt, P.F., Burger, E.H., 1995. Longitudinal analysis of radiographic trabecular pattern by image processing. Bone 17, 527–532.
ARTICLE IN PRESS 580
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Nafei, A., Kabel, J., Odgaard, A., Linde, F., Hvid, I., 2000. Properties of growing trabecular ovine bone. Part II: architectural and mechanical properties. Journal of Bone and Joint Surgery (Br) 82-B, 921–927. Nowak, M., 2006. Structural optimization system based on trabecular bone surface adaptation. Structural Multidisciplinary Optimization 32, 241–249. Parfitt, A.M., 1994. Osteonal and hemi-osteonal remodeling: the spatial and temporal framework for signal traffic in adult human bone. Journal of Cellular Biochemistry 55, 273–286. Pauwels, F., 1965. Gesammelte Abhandlungen zur Funktionellen Anatomic des Bewegungsapparates. Springer, Berlin. Petersson, J., 1999. Some convergence results in perimeter-controlled topology optimization. Computer Methods in Applied Mechanics and Engineering 171, 123–140. Ruimerman, R., Hilbers, P., van Rietbergen, B., Huiskes, R., 2005. A theoretical framework for strain-related trabecular bone maintenance and adaptation. Journal of Biomechanics 38, 931–941. Silva, M.J., Gibson, L.J., 1997. Modeling the mechanical behavior of vertebral trabecular bone: effects of age-related changes in microstructure. Bone 21 (2), 191–199. Skedros, J.G., Baucom, S.L., 2007. Mathematical analysis of trabecular ‘trajectories’ in apparent trajectorial structures: The unfortunate historical emphasis on the human proximal femur. Journal of Theoretical Biology 244, 15–45.
Subbarayan, G., Bartel, D.L., 2000. A reconciliation of local and global models for bone remodeling through optimization theory. Journal of Biomechanical Engineering 122, 72–76. Tanck, E., Homminga, J., van Lenthe, G.H., Huiskes, R., 2001. Increase in bone volume fraction precedes architectural adaptation in growing bone. Bone 28 (6), 650–654. Tovar, A., 2004. Bone remodeling as a hybrid cellular automaton optimization process. Ph.D. Dissertation, Department of Aerospace and Mechanical Engineering, University of Notre Dame, USA. Tsubota, K., Adachi, T., Tomita, Y., 2003. Effects of a fixation screw on trabecular structural changes in a vertebral body predicted by remodeling simulation. Annals of Biomedical Engineering 31, 733–740. Van Rietbergen, B., Huiskes, R., Weinanas, H., Odgaard, A., Kabel, J., 1998. Relationships between bone morphology and bone elastic properties can be accurately quantified using high-resolution computer reconstructions. Journal of Orthopaedic Research 16, 23–28. Ward, F.O., 1838. Outlines of Human Osteology. Henry Renshaw, London. Weinans, H., Huiskes, R., Grootenboer, H.R., 1992. The behavior of adaptive boneremodeling simulation models. Journal of Biomechanics 25, 1425–1441. Whitehouse, W.J., Dyson, E.D., 1974. Scanning electron microscope studies of trabecular bone in the proximal end of the human femur. Journal of Anatomy 118 (3), 417–444.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Computational simulation of trabecular adaptation progress in human proximal femur during growth In Gwun Jang a, Il Yong Kim b, a b
Department of Mechanical and Materials Engineering, Jackson Hall 213, 35 Fifth Field Company Lane, Queen’s University, Kingston, Ontario, Canada K7L 3N6 Department of Mechanical and Materials Engineering, McLaughlin Hall 305, 130 Stuart Street, Queen’s University, Kingston, Ontario, Canada K7L 3N6
a r t i c l e in fo
abstract
Article history: Accepted 18 December 2008
There are a large number of clinical and experimental studies that analyzed trabecular architecture as a result of bone adaptation. However, only a limited amount of quantitative data is currently available on the progress of trabecular adaptation during growth. In this paper, we proposed a two-step numerical simulation method that predicts trabecular adaptation progress during growth using a recently developed topology optimization algorithm, design space optimization (DSO), under the hypothesis that the mechanisms of DSO are functionally equivalent to those of bone adaptation. We applied the proposed scheme to trabecular adaptation simulation in human proximal femur. For the simulation, the full trabecular architecture in human proximal femur was represented by a two-dimensional mFE model with 50 mm resolution. In Step 1, we determined a reference value that regulates trabecular adaptation in human proximal femur. In Step 2, we simulated trabecular adaptation in human proximal femur during growth with the reference value derived in Step 1. We analyzed the architectural and mechanical properties of trabecular patterns through iterations. From the comparison with experimental data in the literature, we showed that in the early growth stage trabecular adaptation was achieved mainly by increasing bone volume fraction (or trabecular thickness), while in the later stage of the development the trabecular architecture gained higher structural efficiency by increasing structural anisotropy with a relatively low level of bone volume fraction (or trabecular thickness). We demonstrated that the proposed numerical framework predicted the growing progress of trabecular bone that has a close correlation with experimental data. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Trabecular adaptation Bone growth Human proximal femur Topology optimization Design space optimization
1. Introduction Although there are a number of clinical and experimental studies on trabecular adaptation, only few quantitative data are currently available on the architectural changes in trabecular adaptation during growth. Korstjens et al. (1995) found that fine trabecular patterns in young children develop into coarser ones as they age. Nafei et al. (2000) studied the change of morphological properties in growing trabecular ovine bone. Based on the experiment with five different age groups of female pigs, Tanck et al. (2001) showed that in the early phase of growth the bone density changes significantly, whereas in the later stage of the development mainly the trabecular architecture undergoes adaptation with a small change in overall density. In an attempt to numerically simulate trabecular adaptation progress, phenomenological update algorithms were used with the aid of the finite element method; these studies sought to
Corresponding author. Tel.: +1 613 533 3077; fax: +1 613 533 6489.
E-mail addresses: [email protected] (I.G. Jang), [email protected] (I.Y. Kim). 0021-9290/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2008.12.009
quantitatively describe the mechanical stimulus (cause) and the bone response (effect). In the phenomenological models, a variety of mechanical stimuli has been utilized: stress (Adachi et al., 2001; Pauwels, 1965), strain (Carter et al., 1981; Cowin and Hegedus, 1976), strain-energy density (Fyhrie and Carter, 1986; Huiskes et al., 1987; Weinans et al., 1992), and strain-energy density rate (Huiskes et al., 2000; Ruimerman et al., 2005). The aforementioned mechanical stimuli used in the different models have been shown to be consistent with one another in boneremodeling simulation. Recently developed phenomenological models with so-called mFE models enable us to simulate bone remodeling or adaptation at the tissue level even in threedimensions (Ruimerman et al., 2005; Tsubota et al., 2003). In these studies, mFE models were built to represent trabecular architecture with the resolution of tens of mm, and in situ mechanical stimuli were used for the regulation of the trabecular adaptation. As a different numerical approach, topology optimization methods (Bendsøe, 1989; Bendsøe and Kikuchi, 1988) have been successfully applied and produced a variety of meaningful results in bone-remodeling simulation (Bagge, 2000; Fernandes and Rodrigues, 1999; Hollister et al., 1994). Recently, the authors developed a new topology optimization scheme that implements
ARTICLE IN PRESS 574
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
an extra perimeter constraint for obtaining trabecular architecture in human proximal femur (Jang and Kim, 2008). The numerical results showed a close correlation with clinical data in the literature. It is interesting to note that phenomenological methods share the common underlying assumption with topology optimizationbased methods, in that the aim of bone remodeling is to use the bone material in an ‘‘optimal’’ way. Their similarities can also be revealed from the view point of fully stressed design formulation, where each member of the structure is fully stressed under the multiple loading conditions for the optimum design (Haftka and Gu¨rdal, 1992). Jang and Kim (2008) recently showed that the trabecular architecture in human proximal femur is constructed as a fully stressed structure to support multi-directional and multi-axial load cases in daily activities. Huiskes (2000) summarized this observation: ‘‘It has been shown that computer algorithms for material design optimization, when given the opportunity, show the tendency to produce trabecular patterns, similar to bone.’’ In addition to this qualitative comparison, the phenomenological and topology optimization-based approaches have been explicitly compared (Harrigan and Hamilton, 1994; Subbarayan and Bartel, 2000; Tovar, 2004), and it has been shown that their structural analogy is remarkable (Nowak, 2006). Especially Jang et al. (2009) derived the mathematical relations among the parameters in topology optimization and phenomenological approach. However, as Hart (2001) pointed out, a critical deficiency of topology optimization is that the numerical approach fails to obtain physiological ‘‘intermediate’’ adaptation process and is, therefore, unable to provide an understanding of the ‘‘entire’’ bone-adaptation process. The aforementioned limitation of topology optimization stems from the fact that the design domain does not change during the optimization process; actual trabecular adaptation involves the sequential change of the trabecular surface via osteoclast resorption and osteoblast deposition (Parfitt, 1994). Due to this limitation of conventional topology optimization, we could obtain only ‘‘final’’ trabecular architecture without any information on ‘‘intermediate’’ trabecular adaptation in our previous paper (Jang and Kim, 2008). Compared to the fixed design domain in conventional topology optimization, a recently developed design space optimization (DSO) incorporates an evolutionary feature in order to adaptively change the design domain (Jang and Kwak, 2006; Kim and Kwak, 2002). The objectives of this paper were to numerically predict the entire trabecular adaptation progress in human proximal femur during growth by using DSO for the first time in the field, and to compare the numerical results to the experimental data in the literature. We developed a two-step simulation procedure: determining a reference value for trabecular adaptation (Step 1) and then simulating trabecular adaptation (Step 2). We showed that in the early growth stage where the trabecular architecture was less optimal, trabecular adaptation was achieved mainly by increasing bone volume fraction (or trabecular thickness), whereas in later stage of the development the trabecular architecture achieved higher structural efficiency by increasing structural anisotropy with a relatively low level of bone volume fraction (or trabecular thickness).
2. Methods 2.1. Design space optimization and its application to trabecular adaptation Design space optimization consists of double loops: the inner loop which is the same as the conventional topology optimization and the outer loop which controls design domain change, as shown in Fig. 1. When an optimization process in the inner loop with the current design domain reaches convergence status, the outer loop seeks out a better design domain by expanding the design space (at the start of the inner loop) and by shrinking the design space (at the end of the inner loop).
Fig. 1. Flow chart of design space optimization (DSO).
Fig. 2. Schematic diagram for trabecular adaptation and design space optimization. This evolutionary nature of DSO allowed us to hypothesize that the mechanisms of DSO are functionally equivalent to those of the trabecular adaptation, as illustrated in Fig. 2. Osteoblast bone formation in trabecular adaptation may be considered to be equivalent to design space expansion in DSO, in which process one layer of new elements with very low density is added along the boundary of the current whole design domain. Osteoclast bone resorption is conceptually similar with design space reduction in the sense that inefficiently used elements are deleted from the current design domain. In this paper, starting from initial trabecular bone (Fig. 3(a)), we added new elements with the lowest density (0.05 in this paper) on all over the boundary of the bone structure (Fig. 3(b)). This means we did not use any specific criterion—such as stress, strain, or strain-energy density requirement—to add elements. After expanding the trabecular bone domain this way, we determined new, balanced morphology using ‘‘optimization’’. Then, inefficiently used elements which have bone densities below 0.1 were removed at the end of inner loop (Fig. 3(c)). The steps shown in Fig. 3(a)–(c) constitute a single iteration in our numerical method. As this numerical routine progressed, trabecular morphology gradually evolved into a more optimized one. 2.2. Two-step simulation for trabecular adaptation 2.2.1. Step 1 The main objective of Step 1 was to determine the reference value k used in the bone-remodeling update equation (Huiskes et al., 2000; Ruimerman et al., 2005;
ARTICLE IN PRESS I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
575
Fig. 3. Application of design space optimization to bone adaptation. (a) Initial trabecular morphology; (b) expanded trabecular morphology right after design space expansion; and (c) resultant trabecular morphology right after design space reduction. element, h1 and h2 are lengths of a finite element (both 50 mm in this paper), n1 and n2 are divisions of the domain along x1 and x2 directions, and ri,j represents the bone density located in the coordinate of (i, j) in the mesh (Petersson, 1999). The mass and perimeter constants in Eq. (2), M0 and P0, were calculated from the initial trabecular structure. Sensitivity equations for Eq. (2) are presented in our previous work (Jang and Kim, 2008). Because the reference value k for bone remodeling in human proximal femur is not known, we approximated the value using the simulation result by Eq. (2) as PN k ffi kapprox ¼
rÞ
i¼1 ðU a;i =
N
,
(3)
which is the average value of Ua,i/r over the trabecular architecture determined at the end of Step 1. It was shown in our previous work that Eq. (2) produced physiological results (Jang and Kim, 2008), and therefore kapprox obtained by Eq. (3) may be considered as an effective estimate for subsequent trabecular adaptation simulation, Step 2. 2.2.2. Step 2 We reformulated Eq. (2) to simulate trabecular adaptation during growth as 3 N X X 1 T cj r i vi uj Kuj þ a 2 j¼1 i¼1
Minimize f ðqÞ ¼ Fig. 4. FE model for an initial trabecular architecture in human proximal femur.
Subject to gðqÞ ¼ P ¼
n2 X n1 X
h2 jri;j ri1;j j
j¼0 i¼1
Weinans et al., 1992)
þ
U a;i dri ¼B k dt ri
Subject to
f ðqÞ ¼
h1 jri;j ri;j1 jXP 0
j¼1 i¼0
(1)
where ri represents the bone density of the ith element, B the constant, and Ua,i the apparent strain-energy density of the ith element. First, the trabecular domain was discretized into finite elements with the size of 50 mm 50 mm. In order to represent initial trabecular architecture, we randomly generated hollow circular patterns with an external diameter of 1400 mm and an internal diameter of 1200 mm (Fig. 4). The circular pattern was chosen to construct initial trabecular architecture with thin and near-isotropic trabecular patterns based on the observation by Nafei et al. (2000). The total number of finite elements needed for the definition of the initial trabecular architecture was 0.75 million. Check more detailed information on FE modeling and analysis in Appendix A. Then, total strain energy was to be minimized under three load cases subject to mass (M) and perimeter (P) constraints. The design variables were the bone densities (r) in the FE elements with the total number N. The bone density vector, r, is composed of each bone density ri of the ith element as q ¼ ½ r1 r2 . . . ri . . . rN T . The perimeter of a structure is the sum of the lengths of all inner and outer boundaries. In this study, a minimum allowable perimeter constraint was applied to effectively represent porosity of the trabecular bone and produce physiological numerical results (Jang and Kim, 2008). Using the information, we formulated the optimization problem as Minimize
n2 X n1 X
3 X 1 T uj Kuj cj 2 j¼1
g 2 ðqÞ ¼ P ¼
N X
i vi pM 0 i¼1 n n 2 1 XX
g 1 ðqÞ ¼ M ¼
r
n2 X n1 X
h1 jri;j ri;j1 jXP 0
j¼1 i¼0
0:05pri p1:0
(4)
a ¼ gk ffi gkapprox .
(5)
The reference value k with Eq. (5) represents the balance between two P3 T competing measures in Eq. (4): total strain energy j¼1 c j ðð1=2Þuj Kuj Þ and total PN bone mass i¼1 ri vi . As k increases, more emphasis is placed on the minimization of the second term (total bone mass), and the optimum solution would have a relatively lower total bone mass and higher total strain energy. It can also be seen from Eq. (1) that a large k will slow down the increase in bone density (accordingly bone mass) over the iterations. Sensitivity equations for Eq. (4) were derived as 3 X @f ðqÞ 1 T @K uj ¼ cj uj þ ave @re 2 @re j¼1
h2 jri;j ri1;j j
j¼0 i¼1
þ
0:05pri p1:0.
where a is to be determined on the basis of kapprox obtained in Step 1. Harrigan and Hamilton (1994) mathematically verified that a global function in the form of Eq. (4) is exactly equivalent to a strain-energy density-based boneremodeling equation, Eq. (1). In terms of bone-mass change during simulation, when using Eq. (2), we have ‘‘constant’’ total bone mass during simulation, which is not physiological in trabecular adaption during growth. With Eq. (4), however, we can vary the total bone mass during simulation while balancing the bone mass and total strain energy so as to simulate physiological trabecular adaptation progress. In our previous study (Jang et al., 2009), we mathematically derived the relationships among the reference value k in Eq. (1), the penalty parameter a in Eq. (4), and the exponent g in the relationship between bone density r and Young’s Modulus E of bone material (E ¼ E0rg) as (check more detailed mathematical derivation for Eq. (5) in Appendix B)
(2)
where cj is the normalized weighting factor for three load cases (c1 ¼ 0.6 for onelegged stance, c2 ¼ 0.2 for abduction, and c3 ¼ 0.2 for adduction), uj is the nodal displacement vector of the jth load case, K stiffness matrix, vi the volume of the ith
n2 X n1 X @ri;j @gðqÞ ¼ h2 sgnðri;j ri1;j Þ @re @re j¼0 i¼1
þ
n2 X n1 X j¼1 i¼0
h1
@ri;j @re
sgnðri;j ri;j1 Þ.
(6)
The overall procedure for this two-step simulation is briefly illustrated in Fig. 5.
ARTICLE IN PRESS 576
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Fig. 5. Schematic diagram for two-step simulation for trabecular adaptation during growth.
Trabecular adaptation in Step 2 resulted in the changes of morphological properties of trabecular patterns: bone volume fraction ratio (BV/TV), trabecular thickness (Tb.Th), structural anisotropy, and strain-energy density. In order to observe the histories of the four indices during growth, four volume of interest (VOIs) with the size of 10 mm 10 mm were selected from the epiphysis of proximal femur, the metaphysis of proximal femur, Ward’s triangle, and intertrochanteric trabecular group. The structural anisotropy was determined from the mean intercept length (MILmax/MILmin) (Whitehouse and Dyson, 1974).
3. Results In Step 1, we obtained kapprox ¼ 1.54 102 J g1, and this value was used to calculate a in Eq. (4) in Step 2. If kapprox obtained from the final trabecular architecture in Step 2 is close to the original kapprox in Step 1, we can verify that the original kapprox was an effective estimate for trabecular adaptation during growth in human proximal femur. Indeed, the two ks should be ideally the same. However, they may have some differences due to numerical issues including round-off errors and error accumulation during simulation. When we approximate the reference value k from the result in Step 1 and put it into Eq. (4) with the form of a in Step 2, it implicitly has the meaning that in Step 2 we are expected to obtain bone architecture with similar bone mass and the reference value k calculated in Step 1, because the derived value a with Eq. (5) represents the regulating measure for the total strain energy and bone mass in trabecular architecture in Eq. (4). We obtained kapprox of 1.56 102 J g1 in Step 2; the difference was only 1.3%. In Step 2, the initial trabecular structure in a human proximal femur adapted to its surrounding mechanical environment through 20 iterations as shown in Fig. 6; the final result presented in Fig. 6(d) is very similar to the actual patterns observed in the
literature as shown in Fig. 7 (Gray, 1918; Skedros and Baucom, 2007). Fig. 8 shows the complete histories of the four characteristics during iterations. The bone volume fraction increased rapidly in the early growth phase, reached its maximum at around iteration 5, decreased slightly, and stabilized from iteration 10 for all VOIs. Trabecular thickness in Fig. 8(b) showed a similar trend, which indicates a strong correlation between bone volume fraction and trabecular thickness. On average the early structural anisotropy was relatively low, and there was only a small change until iteration 5. After this point, anisotropy increased gradually until iteration 15, and then it flattened out (Fig. 8(c)). In Fig. 8(d), we investigated the history of strain-energy density, which is a measure of the efficiency of the structural architecture. Strain-energy density in all VOIs reached their minimum values at iteration 5 and then stayed at the same level.
4. Discussion The trabecular adaptation process may be explained from the view point of bone stiffness. Bone stiffness is determined by three parameters: tissue modulus, bone volume fraction, and bone architecture (Kabel et al., 1999; van Rietbergen et al., 1998). In this study, the tissue modulus was assumed to be constant, and therefore the bone stiffness developed only through the combined effect of bone volume fraction and anisotropy. This suggests that, in the early growth stage (iterations 1–5) when trabecular adaptation was achieved mainly by increasing bone volume fraction (or trabecular thickness) (Fig. 8(a) and (b)), the trabecular architecture was less optimal; on the other hand the trabecular architecture gained higher structural efficiency in the later development stage (iterations 5–15) by increasing structural
ARTICLE IN PRESS I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Fig. 6. Trabecular architecture in human proximal femur during growth.
577
anisotropy (Fig. 8(c)) even with a relatively low level of bone volume fraction (or thickness). General trends of the architectural and mechanical properties during growth in our simulation (Fig. 8) can be compared to the experimental data in the literature (Kneissel et al., 1997; Nafei et al., 2000; Tanck et al., 2001). In Fig. 9, bone volume fraction and anisotropy in the vertebra and proximal tibia of female pigs during 6–230 weeks of age are presented (Tanck et al., 2001). Examining the growth plate in the proximal tibia, Tanck et al. (2001) estimated that 104 weeks in pigs is roughly equivalent to late adolescence (about 15 years) for humans. Their experimental data in Fig. 9(a) shows a close correlation with our numerical data in Fig. 8(a), and therefore we estimated 1 iteration in this study to be 12.5 weeks for pigs or 1.8 years for human. In this study, bone volume fraction and trabecular thickness reached their maximum values at around iteration 5 (approximately 9 years for human), and a similar trend was observed in medieval Nubian skeletons by Kneissel et al. (1997): the highest bone volume fraction and trabecular thickness occurred during adolescence. Nafei et al. (2000) also found that architectural properties of trabecular bone in sheep change significantly with skeletal maturity. The bone volume fraction of the skeletally mature group was 22% higher than that of the skeletally immature group, and the structural anisotropy was 25% higher in the skeletally mature than in the skeletally immature sheep. As Huiskes et al. (2000) stated, homeostasis in a living organism is the property of a control system which regulates its internal parameters so as to maintain a stable, constant condition over possibly wide environmental variations. According to the classical control theory (Franklin et al., 2002), a general timeresponse of the system to any input or disturbance consists of rise time, peak overshoot, and others such as settling time, quarterdecay, etc. It is interesting to observe that trabecular adaptation history obtained in this paper (Fig. 8) shows common features of a control system, and it suggests that the proposed method plays the role of a regulator for bone adaptation. There are very few numerical studies that determined the trabecular adaptation progress of full trabecular architecture, because of computational difficulties. Specifically, with Euler forward integration in phenomenological approaches, the numerical stability theory typically requires that information does not propagate across more than one element per time step (Cook et al., 1989). Hence, the time step for a mFE model must be very small, and usually hundreds of iterations are needed whose computing time is very long (Ruimerman et al., 2005). On the other hand, topology optimization-based methods show better performance in convergence and computing time despite of their non-physiological intermediate progresses (Jang et al., 2009). In this paper, we successfully predicted physiological intermediate trabecular adaptation progress by virtue of evolutionary nature of DSO while improving the overall computational efficiency. We should address the limitation and challenge of this research. Due to a very limited amount of clinical data, it was very difficult to directly compare our numerical results to the literature data in terms of species (human in Kneissel et al.’s paper; sheep in Nafei et al.’s paper; pig in Tanck et al.’s paper), age span (0–60 years in Kneissel et al.’s paper; 3–80 months in Nafei et al.’s paper; 6–230 weeks in Tanck et al.’s paper), and the anatomical site of sampling (vertebra in Kneissel et al.’s paper; proximal tibia in Nafei et al.’s paper; proximal tibia and vertebra in Tanck et al.’s paper). Nevertheless, the common feature on these different samples is that they are major weight-bearing bones in the skeletal systems. The sites of sampling were also chosen from parts which are loaded predominantly in axial compression.
ARTICLE IN PRESS 578
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Fig. 7. Photograph and radiograph of human proximal femur (by courtesy of John G. Skedros; Skedros and Baucom, 2007) (Reprinted with permission from Elsevier).
Fig. 8. Architectural and mechanical properties in human proximal femur during growth. (a) Bone volume fraction (BV/TV); (b) Trabecular thickness (Tb. Th, mm); (c) Structural anisotropy; and (d) Strain energy density (J cm3).
Another limitation is two-dimensional aspects of geometry and loads in the FE model. Our previous paper (Jang and Kim, 2008) showed that final trabecular architecture from topology optimization agreed well with the real trabecular architecture even with two-dimensional models. The good correlation in terms of Ward’s classification (Gray, 1918; Ward, 1838; Whitehouse and Dyson, 1974) and intersection angles of trabecular groups is mainly because the hip joint contact force in the analysis is the greatest load on the mediolateral plane of a proximal femur, and moreover its magnitude is at least an order of magnitude higher
than other loads. Even though we obtained trabecular architecture which is structurally similar to photograph and radiograph images of human proximal femur, no quantitative comparisons could be made between the numerically predicted values of bone volume fraction, trabecular thickness, or structural anisotropy and those measured on bone specimens experimentally. Another consequence of using a two-dimensional model is that the mechanical properties would be more sensitive to reductions in the thickness and number of trabeculae than would be expected for a threedimensional model of trabecular bone, because two-dimensional
ARTICLE IN PRESS I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
579
Acknowledgements The authors would like to thank Dr. Krister Svanberg for providing the MMA code and Dr. John G. Skedros for providing the photograph and radiograph images of human proximal femur. The authors would also like to thank Ryan Willing at Queen’s University (Kingston, Canada) for helpful discussion.
Appendix. Supporting Information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.jbiomech.2008.12.009.
References
Fig. 9. History of trabecular architecture of female pig during growth (Tanck et al., 2001) (Reprinted with permission from Elsevier). (a) Bone volume fraction of trabecular architecture in a female pig through age and (b) anisotropy of trabecular architecture in a female pig through age.
networks have fewer connections, on average, at each node than three-dimensional networks (Silva and Gibson, 1997). The proposed method in this paper is mainly based on topology optimization and is enhanced by adding the evolutionary feature of design space optimization to mimic the bone formation and resorption activities in trabecular adaptation. We performed two-step simulation to obtain the intermediate trabecular architecture throughout trabecular adaptation as well as its final architecture in human proximal femur during growth. From the comparison with clinical and experimental data, it is shown that the proposed numerical framework predicted physiological trabecular adaptation progress during growth. It means that, as other researchers already found, topology optimizationbased methods are compatible with the phenomenological methods in bone-remodeling simulation. Moreover, the proposed method can be used to predict the progress of trabecular adaptation unlike other conventional topology optimizationbased approaches. The trends of the architectural and mechanical properties obtained by the simulation would improve the insight into trabecular adaptation during growth. However, with the twodimensional aspects of the model in this paper, we could not quantitatively compare morphometrical indices from our simulation with those from experiments. For thorough and exact comparison, it is required to perform trabecular adaptation simulation with the corresponding three-dimensional model, which would be our future work.
Conflict of interest It is stated that the authors have no conflict of interest including any financial and personal relationships with other people or organizations that could inappropriately influence their work.
Adachi, T., Tsubota, K., Tomita, Y., Hollister, S.J., 2001. Trabecular surface remodeling simulation for cancellous bone using microstructural voxel finite element models. Journal of Biomechanical Engineering 123, 403–409. Bagge, M., 2000. A model of bone adaptation as an optimization process. Journal of Biomechanics 33, 1349–1357. Bendsøe, M.P., 1989. Optimal shape design as a material distribution problem. Structural Optimization 1, 193–303. Bendsøe, M.P., Kikuchi, N., 1988. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 71, 197–224. Carter, D.R., Harris, W.H., Vasu, R., Caler, W.E., 1981. The mechanical and biological response of cortical bone to in vivo strain histories. In: Cowin, S.C. (Ed.), Mechanical Properties of Bone, American Society of Mechanical Engineers. Cook, R.D., Malkus, D.S., Plesha, M.E., 1989. Concepts and Applications of Finite Element Analysis, third ed. Wiley, New York. Cowin, S.C., Hegedus, D.H., 1976. Bone remodeling I: a theory of adaptive elasticity. Journal of Elasticity 6, 313–326. Fernandes, P.R., Rodrigues, H.C., 1999. A material optimization model for bone remodeling around cementless hip stems. In: Proceedings of the Ninth European Conference on Computational Mechanics, Munich, Germany. Franklin, G.F., Powell, J.D., Emami-Naeini, A., 2002. Feedback Control of Dynamic Systems. Prentice-Hall, New Jersey. Fyhrie, D.P., Carter, D.R., 1986. A unifying principle relating stress to trabecular bone morphology. Journal of Orthopaedic Research 4, 304–317. Gray, H., 1918. Anatomy of the human body. Philadelphia, Lea & Febiger. (Revised and re-edited by Lewis W.H., 2000, twentieth ed., Bartleby.com, New York. Haftka, R.T., Gu¨rdal, Z., 1992. Elements of Structural Optimization, third ed. Kluwer Academic Publishers, Dordrecht. Harrigan, T.P., Hamilton, J.J., 1994. Bone remodeling and structural optimization. Journal of Biomechanics 27, 323–328. Hart, R.T., 2001. Bone modeling and remodeling: theories and computation. In: Cowin, S.C. (Ed.), Bone Mechanics Handbook. CRC Press, Boca Raton, FL. Hollister, S.J., Brennan, J.M., Kikuchi, N., 1994. A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress. Journal of Biomechanics 27, 433–444. Huiskes, R., 2000. If bone is the answer, then what is the question? Journal of Anatomy 197, 145–156. Huiskes, R., Ruimerman, R., van Lenthe, G.H., Janssen, J.D., 2000. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 405, 704–706. Huiskes, R., Weinans, H., Grootenboer, J., Dalstra, M., Fudala, M., Slooff, T.J., 1987. Adaptive bone remodelling 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. Jang, I.G., Kim, I.Y., Kwak, B.M., 2009. Analogy of strain energy density based boneremodeling algorithm and structural topology optimization. Journal of Biomechanical Engineering 131 (1), 011012. Jang, I.G., Kwak, B.M., 2006. Evolutionary topology optimization using design space adjustment based on fixed grid. International Journal for Numerical Methods in Engineering 66 (11), 1817–1840. Kabel, J., Odgaard, A., van Rietberge, B., Huiskes, R., 1999. Connectivity and elastic properties of cancellous bone. Bone 24, 115–120. Kim, I.Y., Kwak, B.M., 2002. Design space optimization using a numerical design continuation method. International Journal for Numerical Methods in Engineering 53, 1979–2002. Kneissel, M., Roscher, P., Steiner, W., Schamall, D., Kalchhauser, G., Boyde, A., TeschlerNicola, M., 1997. Cancellous bone structure in the growing and aging lumbar spine in a historic Nubian population. Calcified Tissue International 61, 95–100. Korstjens, C.M., Geraets, W.G.M., van Ginkel, F.C., Prahl-Andersen, B., van Der Stelt, P.F., Burger, E.H., 1995. Longitudinal analysis of radiographic trabecular pattern by image processing. Bone 17, 527–532.
ARTICLE IN PRESS 580
I.G. Jang, I.Y. Kim / Journal of Biomechanics 42 (2009) 573–580
Nafei, A., Kabel, J., Odgaard, A., Linde, F., Hvid, I., 2000. Properties of growing trabecular ovine bone. Part II: architectural and mechanical properties. Journal of Bone and Joint Surgery (Br) 82-B, 921–927. Nowak, M., 2006. Structural optimization system based on trabecular bone surface adaptation. Structural Multidisciplinary Optimization 32, 241–249. Parfitt, A.M., 1994. Osteonal and hemi-osteonal remodeling: the spatial and temporal framework for signal traffic in adult human bone. Journal of Cellular Biochemistry 55, 273–286. Pauwels, F., 1965. Gesammelte Abhandlungen zur Funktionellen Anatomic des Bewegungsapparates. Springer, Berlin. Petersson, J., 1999. Some convergence results in perimeter-controlled topology optimization. Computer Methods in Applied Mechanics and Engineering 171, 123–140. Ruimerman, R., Hilbers, P., van Rietbergen, B., Huiskes, R., 2005. A theoretical framework for strain-related trabecular bone maintenance and adaptation. Journal of Biomechanics 38, 931–941. Silva, M.J., Gibson, L.J., 1997. Modeling the mechanical behavior of vertebral trabecular bone: effects of age-related changes in microstructure. Bone 21 (2), 191–199. Skedros, J.G., Baucom, S.L., 2007. Mathematical analysis of trabecular ‘trajectories’ in apparent trajectorial structures: The unfortunate historical emphasis on the human proximal femur. Journal of Theoretical Biology 244, 15–45.
Subbarayan, G., Bartel, D.L., 2000. A reconciliation of local and global models for bone remodeling through optimization theory. Journal of Biomechanical Engineering 122, 72–76. Tanck, E., Homminga, J., van Lenthe, G.H., Huiskes, R., 2001. Increase in bone volume fraction precedes architectural adaptation in growing bone. Bone 28 (6), 650–654. Tovar, A., 2004. Bone remodeling as a hybrid cellular automaton optimization process. Ph.D. Dissertation, Department of Aerospace and Mechanical Engineering, University of Notre Dame, USA. Tsubota, K., Adachi, T., Tomita, Y., 2003. Effects of a fixation screw on trabecular structural changes in a vertebral body predicted by remodeling simulation. Annals of Biomedical Engineering 31, 733–740. Van Rietbergen, B., Huiskes, R., Weinanas, H., Odgaard, A., Kabel, J., 1998. Relationships between bone morphology and bone elastic properties can be accurately quantified using high-resolution computer reconstructions. Journal of Orthopaedic Research 16, 23–28. Ward, F.O., 1838. Outlines of Human Osteology. Henry Renshaw, London. Weinans, H., Huiskes, R., Grootenboer, H.R., 1992. The behavior of adaptive boneremodeling simulation models. Journal of Biomechanics 25, 1425–1441. Whitehouse, W.J., Dyson, E.D., 1974. Scanning electron microscope studies of trabecular bone in the proximal end of the human femur. Journal of Anatomy 118 (3), 417–444.