
==== Front
Theor Biol Med ModelTheoretical Biology & Medical Modelling1742-4682BioMed Central London 1742-4682-2-61573332810.1186/1742-4682-2-6ResearchVolume-based non-continuum modeling of bone functional adaptation Wang Zhengyuan 1wzhengyuan@gmail.comMondry Adrian 1mondry@hotmail.com1 Medical and Clinical Informatics Group, Bioinformatics Institute, #07-01 Matrix, 30 Biopolis Street, 138671 Singapore2005 28 2 2005 2 6 6 28 9 2004 28 2 2005 Copyright © 2005 WANG and MONDRY; licensee BioMed Central Ltd.2005WANG and MONDRY; licensee BioMed Central Ltd.This is an Open Access article distributed under the terms of the Creative Commons Attribution License (), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Background
Bone adapts to mechanical strain by rearranging the trabecular geometry and bone density. The common finite element methods used to simulate this adaptation have inconsistencies regarding material properties at each node and are computationally demanding. Here, a volume-based, non-continuum formulation is proposed as an alternative. Adaptive processes corresponding to various external mechanical loading conditions are simulated for the femur.

Results
Bone adaptations were modeled for one-legged stance, abduction and adduction. One-legged stance generally results in higher bone densities than the other two loading cases. The femoral head and neck are the regions where densities change most drastically under different loading conditions while the distal area always contains the lowest densities regardless of the loading conditions. In the proposed formulation, the inconsistency of material densities or strain energy densities, which is a common problem to finite element based approaches, is eliminated. The computational task is alleviated through introduction of the quasi-binary connectivity matrix and linearization operations in the Jacobian matrix and is therefore computationally less demanding.

Conclusion
The results demonstrated the viability of the proposed formulation to study bone functional adaptation under mechanical loading.
==== Body
Background
Much research effort has been devoted to understanding the functional adaptation of bone under physiological loading ever since the idea of bone functional adaptation was proposed by Wolff more than one hundred years ago [1-14]. Various computational models have been put forward in the past decades and the methods describing the changing rate of bone density corresponding to strain energy density, with finite element implementation, have become the most popular of them [6,15-29].

The common finite element approach is to take the element densities as the state variables and define elements with either constant or varying densities, then update the material densities for the next step of computation according to the computed strain energy density [22,23,26,27]. With more and more powerful desktop computers and commercial finite element analysis software available, this approach is widely used today. Yet some specific problems of this approach are not well addressed so far, although decades have passed, and the numerical results are inevitably affected.

One common problem is the inconsistency of material densities on element boundaries [30]. During the updating of material densities in each step, different elements may take different densities due to strain energy density, thus often leading to conflicting material properties at the boundaries shared by more then one element. Since this conflict affects the integration points, which always come from the element boundaries, the errors are carried forward and cannot be eliminated by smoothing techniques. So it is not surprising that, if the program is allowed to run for a certain time, most of the elements tend to become either saturated or completely resorbed, leading to checker-board patterns especially in the proximal area of the femur [30].

Some effort has been made trying to solve this problem. For example, a node-based variant of the finite element method was tried with focus on the densities of the nodes rather than densities of the whole elements [21,30]. The node densities are then interpolated across the whole elements before the next step of computing begins. The results are improved, but the stress and strain quantities are still conflicting at the nodes.

Other previous work has used Voronoi structures [31] to study the effects of crack growth on trabecular bone, tapered strut models [32] to study the ageing effect through a parametric approach or continuum FEM [33] to compute the strain energy density in order to overcome individual drawbacks of the common method described. Their potential impact on the formula proposed here is discussed below.

The long existing problems and the limitations of assuming a continuum drive this new effort to explore the possibility of a non-continuum formulation of bone functional adaptations through nodal analysis in the hope of eliminating the errors present in the previous approaches. In the proposed non-continuum formulation, neighboring nodes are connected by struts that are defined with invariant material densities with respect to time but strut volumes are defined as state variables indicating different configurations of bone structure. The updating of strut volume will depend on the strain energy density in the strut in the previous step. As a result, there is no conflict either in density or in strain energy density. The shift of state variables from bone densities to bone volumes not only eliminates the errors inherent to the density-based finite element approaches but also transforms the continuum formulation to a non-continuum formulation [34]. The advantages of a volume-based non-continuum formulation may be appreciated by looking at the bone volume ratios in osteoporotic bones. The ever-increasing resolution of modern imaging techniques now allows to take a much closer look at the trabecular structure of the bone. In the trabecular network, trabeculae with different lengths and thicknesses are connected with each other to form a scaffold serving both mechanical and biological functions [35,36]. They are well connected in normal bones but poorly connected in osteoporotic bones in addition to reduced thickness. To characterize the trabecular structure, two terms are often used: bone volume/tissue volume (BV/TV) ratio and bone material orientation [15,25]. Although the cortical bone is densely packed with mineralized material, the trabecular bone dominates the inside space of the bone, highly exposed to bone marrow, highly distributed in volume, and highly involved in bone remodeling. The ratio of trabecular bone volume over tissue volume can be below 30% in osteoporotic bones, which means most space is taken up by void or bone marrow and this questions the appropriateness of a continuum formulation [35]. Besides elimination of the errors mentioned earlier, the small physiological range of bone deformation during normal activities allows linearization operations in the volume-based non-continuum formulation. This saves computation time and alleviates the high demand on hardware resources.

The proposed volume-based non-continuum formulation shows computational advantages in modeling bone functional adaptations and has much potential for clinical applications in this field.

Results
Changes in trabecular structure
Fig. 1, 2, 3 and 4 show the node-based representation of bone adaptation results according to loading cases of one-legged stance, abduction, adduction and the combined loading case respectively.

Figure 1 Adaptation Results: one-legged stance. Results of bone functional adaptation. The color bar shows percentage of actual bone density against maximum bone density, which is 1.74 g/cm3. The density of bone structure is not indicated by the number of sample nodes selected in that region, but by the density (converted from volume) of each node, which is expressed as degree of "red" in this illustration.

Figure 2 Adaptation Results: abduction. Results of bone functional adaptation. The color bar shows percentage of actual bone density against maximum bone density, which is 1.74 g/cm3. The density of bone structure is not indicated by the number of sample nodes selected in that region, but by the density (converted from volume) of each node, which is expressed as degree of "red" in this illustration.

Figure 3 Adaptation Results: adduction. Results of bone functional adaptation. The color bar shows percentage of actual bone density against maximum bone density, which is 1.74 g/cm3. The density of bone structure is not indicated by the number of sample nodes selected in that region, but by the density (converted from volume) of each node, which is expressed as degree of "red" in this illustration.

Figure 4 Adaptation Results: the combined loading case. Results of bone functional adaptation. The color bar shows percentage of actual bone density against maximum bone density, which is 1.74 g/cm3. The density of bone structure is not indicated by the number of sample nodes selected in that region, but by the density (converted from volume) of each node, which is expressed as degree of "red" in this illustration.

The nodal density is a percentage relative to its maximum value (0~100%). It is 'converted' based on volume information (BV/TV) for the purpose of easy visual inspection. As suggested by Zhu [34], the value of 1.74 g/cm3 is used as the maximum value in the present formulation.

In the combined loading case and one-legged stance, the nodal densities are generally high in the proximal area where a considerable number of nodes reach the highest density due to the relatively high load, while in the large distal area, the densities become lower. In the case of abduction, the external load is the smallest out of the three loading cases and the nodal density in this case seldom reaches the maximum. Although the high densities also appear in the proximal area of femoral head, the densities are lower than those of other loading cases. In the large distal area, the node densities generally range at very low levels. In the case of adduction, the highest densities still appear in the proximal area of femoral head, but in the large diaphyseal area, the densities range at very low levels.

In summary, one-legged stance and the combined loading case generally result in higher bone densities than the other two loading cases due to higher mechanical loading. The femoral head and neck are the regions where densities change most drastically under different loading conditions while the distal area always contains the lowest densities regardless of the loading conditions.

Program convergence and bone fracture probability
Starting from a uniform material distribution, the program computes the strain energy density, adjusts the strut configurations, and continues on to the next iteration. Fig. 5 shows the convergence of the model, which demonstrates the adapting process of the model evaluated with the normalized error residue against the final solution. The solution vector starts with a state of uniform material distribution, it then moves toward the final state with uniform strain energy density in the solution space. As the residue of solution drops, the residue of bone volume/tissue volume ratio also drops toward trivial while the program progresses.

Figure 5 Program Convergence. Program convergence evaluated by residues between successive solutions.

A prediction of bone fracture risk has been proposed using stress range levels [37] from which analysis of the material fatigue strength under given loading conditions can be derived. The material is subjected to a range of stress levels due to external load, then the stress leads to fatigue damage and finally leads to collapse of the material. This prediction method was used to estimate the fracture probability of the simulation described here. The BV/TV ratios, stress levels and fracture probability are shown in Fig. 6, 7 and 8 respectively. While the program progresses toward the final solution, the volume of bone material used, indicated by the BV/TV ratio, increases by a few percents then slowly decreases again; meanwhile, the stress level moves down to a low level in the final phase and accordingly, the estimated fracture probability drops from around 90% to around 2% in the final step. It is safe to say that the program finally simulates a reasonable configuration of bone internal structure as a result of the physiological adaptation process.

Figure 6 BV/TV Evolution During Adaptation Process. BV/TV evolution during adaptation process.

Figure 7 Evolution of Apparent Normal Stress Level Within Bone Tissue. Apparent stress level evolution during adaptation process.

Figure 8 Fracture Probability. Bone facture probability during adaptation process.

Discussion
A volume-based non-continuum formulation has been developed that describes the adaptation of bone to various mechanical loading situations. In the finite element approach to bone adaptation simulation, the integration of the entries in the coefficient matrix can be a heavy computing task [22,23,26,27]. In the model proposed here, this is alleviated through the introduction of the simpler connectivity matrix and Jacobian matrix [38]. In daily physiological activities, bone deformations are small and the Jacobian matrix can be linearized.

As mentioned in the introduction, one common problem is the inconsistency of material densities or strain energy densities on element boundaries [30]. Since all the integration points come from the boundaries, the resulting errors will essentially affect the computation. In the improved node-based implementation of the finite element method, the stress and strain quantities are still conflicting at the nodes. In the volume-based non-continuum formulation proposed here, this conflict is eliminated. The material density and strain energy density are all consistent in each individual strut, and the computing becomes less demanding.

As described above, previous work [31-33] addresses some of the weaknesses of the common FEM model. Makiyama [31] employed the "Voronoi structure" to study the effects of crack growth on trabecular bone. The method for generating a Voronoi structure could be quite useful when it calibrates the artificially constructed structure against the physical trabecular structure scanned from a patient. This might then serve as the starting state of bone configuration before adaptation begins. Moore [39] proposed a model to replace the partially damaged trabecula with another trabecula reduced in thickness. If this concept is combined with that of the strut structure, one may also derive the model proposed here, that is, a strut model with either varied modulus due to bone mineralization or adaptive cross-section/volume or even tapered struts as proposed by Kim [32].

Hip fracture is one typical manifestation of osteoporosis, and the results obtained by the simulation indicate that considerable changes of bone structure take place in the regions of femoral head and neck, where the stress level is normally higher than that of distal regions. The variations in stress level as shown in Fig. 7 reflect the adaptive process of the bone internal structure and different structural configurations will yield different stress levels in spite of little change in bone volume / tissue volume ratio.

In the current literature, the time scale for adaptive processes is not very well defined. This general lack of knowledge poses a problem for any experimental proof of concept – while the numbers of strain repetition can be predefined, they must be done within biologically suitable time frames. If a given strain comes too sudden, the bone may break instead of remodel; if the strain is applied over too long a period, it may not be a sufficient stimulus to activate adaptive processes. The lack of well-defined temporal constraints, however, is common. Kim's approach [32] is very interesting in as much as it may allow to integrate the effect of time in the model proposed here. At present, however, there is not sufficient data available to allow to integrate time effects of the clinically interesting mid-range scale, i.e. weeks to months. Kim's model looks at the process of ageing of 35 years and more.

The simulation model presented here may, beyond theoretical calculations, be applied to look at two clinical questions. Firstly, the simulation can be adjusted so that a realistic density distribution is the starting point, and outcomes following certain loading conditions, such as a predefined number of load cycles can then be predicted. Secondly, the program can integrate the measured bone density of a given patient to estimate the fracture risk based on stress level calculations.

Conclusion
By eliminating the common inconsistencies at each node, the formulation presented here shows good numerical performance and successfully predicts reasonable bone structure changes under different loading conditions. It is viable to serve as an alternative method apart from the traditional finite element based approached to study bone adaptations. In conclusion, the volume based non-continuum formulation is a new approach to bone adaptation study and has its own advantages.

Methods
Volume-based representation of the trabecular bone structure
In the volume-based non-continuum formulation used here, the trabecular structure is represented by a connected strut system and each strut can take different sizes according to the mechanical loading requirements, that is, strain energy density. The strut representation is shown in Fig. 9, which resembles a small volume of the trabecular structure. In this setting, the BV/TV ratio can be directly obtained from the ratio of the strut volumes over the unit volume, and material orientation  can be obtained though the resultant of the vectorial material components of the struts, as described by equation (1) and (2).

Figure 9 Bone Structure Decomposition. Representation of bone material by struts. Each strut can assume different geometric dimensions and material properties. Apparent mechanical property of bone material is based on the strut configuration.





where vi is the volume of the i-th strut in the j-th basic unit, Vj is the volume of the j-th basic unit, Rj is the material orientation of the j-th basic unit,  is the orientation of the i-th strut in the j-th basic unit and N is the number of struts in the j-th basic unit.

For the formulation proposed here, the bone structure is decomposed into and represented by a connected network of struts. These struts are a mathematical abstraction of the physical bone structure, from which the BV/TV ratio can be derived. The thickness of a strut is adapted during the bone adaptation process in mimicry of the physiological processes.

Volume adaptation under mechanical loading
The bone mass will vary under mechanical loading. In engineering, the general relationship between varying mass, density and volume is described as:



Since the density here is taken as constant regarding time, the second term on the right hand side, , simply vanishes and the mass variation is realized through volume variation under mechanical loading. Based on the density-based adaptation proposed by Zhu X. et al (39), which is stated as:



the volume-based adaptation can thus be stated as:



where βi = Ui / ρik, which is a comparative coefficient describing the comparison of a given mechanical stimulus in each sensor cell with the reference value k, and Ui represents the strain energy density for the I-th sensor unit; N is the number of sensor cells and fi(x) is the spatial influence function; B(t) is a remodeling coefficient; α indicates the remodeling power of strain energy density [34].

Non-Continuum formulation
With the whole bone represented by a volume-based strut system, the non-continuum formulation can be noted as follows:



where A is a connectivity matrix describing the connecting relationship between the struts, α is the linearized Jacobian matrix,  is the nodal displacement vector to be solved for and  is the loading vector derived from the external mechanical load. The generalized conjugate residue method is used to solve this formulation [38,40].

Connectivity matrix A is the matrix to show the relationship between connected struts with the entries of 1, -1 or 0. A strut starts from the node with the index corresponding to the entry 1 and ends at the node with the index corresponding to -1. It is further illustrated in Fig. 10 and equation (7).

Figure 10 Connectivity of Struts. Physical connectivity relationship between struts is indicated by the connectivity matrix.



Finally, the different loading conditions to be applied are shown in Fig. 11.

Figure 11 Loading cases. Quantitative information of different loading cases. Four loading cases are considered: one-legged stance, abduction, adduction and the combined loading case (weighted based on their respective daily occurrence cycles).

Competing interests
The author(s) declare that they have no competing interests.

Authors' contributions
Zhengyuan Wang developed the formulation and partly prepared the manuscript, Adrian Mondry participated in the adaptation controls and partly prepared the manuscript.

Acknowledgements
This project is supported by the BioMedical Research Council of Agency for Science, Technology and Research, Singapore. Thanks also go to SMA5211 lecturers from Singapore-MIT Alliance, for helpful advice on nodal formulation.
==== Refs
YP Arramon SC Cowin G Luo AM Sadegh D Zhang  Strain rate indicated as a remodeling stimulus by animal experiments Trans Orthop Res Soc 1994 19 280 
DB Burr RB Martin  Errors in bone remodeling: toward a unified theory of metabolic bone disease Am J Anat 1989 186 186 216 2683714 10.1002/aja.1001860208 
DR Carter  Mechanical loading history and skeletal biology J Biomech 1987 20 1095 1109 3323201 10.1016/0021-9290(87)90027-3 
DR Carter TE Orr DP Fyshrie  Relationships between loading history and femoral cancellation bone architecture J Biomech 1989 22 231 244 2722894 10.1016/0021-9290(89)90091-2 
TJ Chambers TN Gardner A Turner-Smith JWM Chow  Induction of bone formation in rat tail vertebrae by mechanical loading Bone and Mineral 1993 20 167 178 8453332 
HT Charles A Vital MVP Ramana  A uniform strain criterion for trabecular bone adaptation: Do continuum-level strain gradient drive adaptation? J Biomech 1997 30 555 563 9165388 10.1016/S0021-9290(97)84505-8 
RL Duncan CH Turner  Mechanotransduction and the functional response of bone to mechanical strain Cal Tissue Int 1995 5 344 358 10.1007/BF00302070 
HM Frost  Skeletal structural adaptations to mechanical usage (SATMU): Redefining Wolff's law: THe remodeling problem Anat Rec 1989 6 414 422 
RE Guldberg M Richards NJ Caldwell CL Kuelske SA Goldstein  Trabecular bone adaptation to variations in porous-coated implant topology J Biomech 1997 30 147 153 9001935 10.1016/S0021-9290(96)00106-6 
MG Mullender R HUiskes  A proposal for regulatory mechanism of Wolff's law J Orthop Res 1995 13 503 512 7674066 10.1002/jor.1100130405 
TE Orr  The role of mechanical stress in bone remodeling 1990 , Stanford University 
H Weinans PJ Prendergast  Tissue adaptation as a dynamical process far from   equilibrium Bone 1996 19 143 149 8853858 10.1016/8756-3282(96)00143-3 
G Subbarayan DL Bartel  Fully stressed and minimum mass structures in bone remodeling 1989 173 176 
J Wolff  Das Gesetz der Transformation der Knochen 1892 Hirschwald Berlin, Hirschwald 
M Baggie  A model of bone adaptation as an optimization process J Biomech 2000 33 1349 1357 10940393 10.1016/S0021-9290(00)00124-X 
GS Beaupre TE Orr DR Carter  AN approach for time-dependent boen remodeling and remodeling application: a preliminary remodeling simulation J Orthop Res 1990 8 662 670 2388106 10.1002/jor.1100080507 
GS Beaupre TE Orr DR Carter  An approach for time-dependent bone remodeling and remodeling  - theoretical development: a preliminary remodeling simulation J Orthop Res 1990 8 651 661 2388105 10.1002/jor.1100080506 
OT Chavez TM Keaveny  High-resolution finite element models of trabecular bone: the dependence of tissue strains and apparent modulus on image resolution: ; Stanford, CA. 1995 23 24 
R Huiskes H Weinans HJ Grootenboer  The behavior ofadaptive bone remodeling simulation models J Biomech 1992 25 1425 1441 1491020 10.1016/0021-9290(92)90056-7 
R Huiskes H Weinans HJ Grootenboer M Dalstra B Fudala TJ Sloof  Adaptive bone remodeling theory applied to prosthetic-design analysis J Biomech 1987 20 1135 1150 3429459 10.1016/0021-9290(87)90030-3 
CR Jacobs ME Levenson GS Beaupre JC Simo  J Middleton, GN Pande and KR Williams  A new implementation of finite element-based bone remodeling Computer Methods in Biomechanics & Biomedical Engineering 1992 Clydach, Books  & Journals INternational 
JH Keyak JM Meagher HB Skinner CD Mote  Automated three-dimensional finite element modeling of bone: A newmethod J Biomech Eng 1990 12 389 397 
JH Keyak HB Skinner  Three-dimensional finite element modeling of bone: Effects of element size J Biomech Eng 1992 14 483 489 
MG Mullender R Huiskes H Weinans  A physiological approach to the simulation of bone remodeling as a self-organizational control process J Biomech 1994 27 1389 1394 7798290 10.1016/0021-9290(94)90049-3 
K Tsubota A Taiji T Yoshihiro  Functional adaptation of cancellous bone in human proximal femur predicted by trabecular surface remodeling simulation toward uniform stress state J Biomech 2002 25 1541 1551 10.1016/S0021-9290(02)00173-2 
B van Rietbergen S Majumdar D Newitt B MacDonald  High resolution MRI and micro-FE for the evaluation of changes in bone mechanical properties during longitudinal clinical trials: application to calcaneal boen in postmenopausal woman after one year of idoxifene treatment Clinical Biomech 2002 17 81 88 10.1016/S0268-0033(01)00110-3 
B van Rietbergen H Weinans R Huiskes A Odgaard  A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element methods J Biomech 1995 28 69 81 7852443 10.1016/0021-9290(95)80008-5 
B van Rietbergen H Weinans BJ Polman R Huiskes  A fast solving method for large-scale FE-models generated from computer images, based on a row-by-row matrix-vector multiplication scheme ASME CED 1994 6 47 52 
H Weinans R Huiskes HJ Grootenboer  Convergence and uniqueness of adaptive bone remodeling Trans Orthop Res Soc 1989 14 310 
CR Jacobs EL Marc SB Gary CS Juan RC Dennis  Numerical instibilities in bone remodeling simulations: The advances of a node-based finite element approach J BIomech 1994 28 449 459 10.1016/0021-9290(94)00087-K 
AM Makiyama S Vajjhala LJ Gibson  Analysis of Crack Growth in a 3D Voronoi Structure: A Model for Fatigue in Low Density Trabecular Bone J Biomech Eng 2002 124 512 520 12405593 10.1115/1.1503792 
HS Kim STS Al-Hassani  A morphological model of vertebral trabecular bone J Biomech 2002 35 1101 1114 12126669 10.1016/S0021-9290(02)00053-2 
TS Smith RB Martin M Hubbard BK Bay  Surface Remodeling of Trabecular Bone Using a Tissue Level Model Journal of Orthopaedic Research 1997 15 593 600 9379270 10.1002/jor.1100150416 
XH Zhu H Gong D Zhu BZ Gao  A study of the effect of non-linearities in the euqation of bone remodeling J Biomech 2002 35 951 960 12052397 10.1016/S0021-9290(02)00028-3 
JA Kanis  Kanis JA  Maintenance of Bone Mass Textbook of osteoporosis 1996 Cambridge, Blackwell Science 
MC de Vernejoule JM Pierre  BD Tilman and S Isidro  Chapter 1: New aspects of bone biology The spectrum of renal osteodystrophy 2001 Oxford, Oxford Univ Press 
D Taylor JH Kuiper  The prediction of stress fractures using a 'stressed volume' concept J Orthop Res 2001 19 919 926 11562142 10.1016/S0736-0266(01)00009-2 
J White  SAM5211: Introduction to simulation High Performance Computation for Engineered Systems 2001 , Singapore-MIT Alliance 
TL Arthur Moor LJ Gibson  Modeling Modulus Reduction in Bovine Trabecular Bone Damaged in Compression J Biomech Eng 2001 123 613 622 11783733 10.1115/1.1407828 
H William SAT Saul TV William PF Brian  Numerical Recipes in C++ The Art of Scientific Computing 2002 2nd Cambridge, Cambridge University Press

