# Volume-based non-continuum modeling of bone functional adaptation

- Zhengyuan Wang
^{1}and - Adrian Mondry
^{1}Email author

**2**:6

https://doi.org/10.1186/1742-4682-2-6

© WANG and MONDRY; licensee BioMed Central Ltd. 2005

**Received: **28 September 2004

**Accepted: **28 February 2005

**Published: **28 February 2005

## Abstract

### 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.

## 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

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/cm^{3} 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

## 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

where *v*_{
i
}is the volume of the *i*-th strut in the *j*-th basic unit, *V*_{
j
}is the volume of the *j*-th basic unit, *R*_{
j
}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:

*et al*(39), which is stated as:

the volume-based adaptation can thus be stated as:

where *β*_{
i
}= *U*_{
i
}/ *ρ*_{
i
}*k*, which is a comparative coefficient describing the comparison of a given mechanical stimulus in each sensor cell with the reference value *k*, and *U*_{
i
}represents the strain energy density for the *I*-th sensor unit; *N* is the number of sensor cells and *f*_{
i
}(*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].

## Declarations

### 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.

## Authors’ Affiliations

## References

- 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-Google Scholar
- DB Burr, RB Martin: Errors in bone remodeling: toward a unified theory of metabolic bone disease. Am J Anat. 1989, 186: 186-216. 10.1002/aja.1001860208.View ArticleGoogle Scholar
- DR Carter: Mechanical loading history and skeletal biology. J Biomech. 1987, 20: 1095-1109. 10.1016/0021-9290(87)90027-3.View ArticleGoogle Scholar
- DR Carter, TE Orr, DP Fyshrie: Relationships between loading history and femoral cancellation bone architecture. J Biomech. 1989, 22: 231-244. 10.1016/0021-9290(89)90091-2.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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. 10.1016/S0021-9290(97)84505-8.View ArticleGoogle Scholar
- RL Duncan, CH Turner: Mechanotransduction and the functional response of bone to mechanical strain. Cal Tissue Int. 1995, 5: 344-358. 10.1007/BF00302070.View ArticleGoogle Scholar
- HM Frost: Skeletal structural adaptations to mechanical usage (SATMU): Redefining Wolff's law: THe remodeling problem. Anat Rec. 1989, 6: 414-422.Google Scholar
- 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. 10.1016/S0021-9290(96)00106-6.View ArticleGoogle Scholar
- MG Mullender, R HUiskes: A proposal for regulatory mechanism of Wolff's law. J Orthop Res. 1995, 13: 503-512. 10.1002/jor.1100130405.View ArticleGoogle Scholar
- TE Orr: The role of mechanical stress in bone remodeling. 1990, , Stanford UniversityGoogle Scholar
- H Weinans, PJ Prendergast: Tissue adaptation as a dynamical process far from equilibrium. Bone. 1996, 19: 143-149. 10.1016/8756-3282(96)00143-3.View ArticleGoogle Scholar
- G Subbarayan, DL Bartel: Fully stressed and minimum mass structures in bone remodeling. 1989, 173-176.Google Scholar
- J Wolff: Das Gesetz der Transformation der Knochen. 1892, Hirschwald Berlin, HirschwaldGoogle Scholar
- M Baggie: A model of bone adaptation as an optimization process. J Biomech. 2000, 33: 1349-1357. 10.1016/S0021-9290(00)00124-X.View ArticleGoogle Scholar
- 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. 10.1002/jor.1100080507.View ArticleGoogle Scholar
- 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. 10.1002/jor.1100080506.View ArticleGoogle Scholar
- 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.Google Scholar
- R Huiskes, H Weinans, HJ Grootenboer: The behavior ofadaptive bone remodeling simulation models. J Biomech. 1992, 25: 1425-1441. 10.1016/0021-9290(92)90056-7.View ArticleGoogle Scholar
- 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. 10.1016/0021-9290(87)90030-3.View ArticleGoogle Scholar
- CR Jacobs, ME Levenson, GS Beaupre, JC Simo: A new implementation of finite element-based bone remodeling. Computer Methods in Biomechanics & Biomedical Engineering. Edited by: J Middleton, GN Pande and KR Williams. 1992, Clydach, Books & Journals INternationalGoogle Scholar
- 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.View ArticleGoogle Scholar
- JH Keyak, HB Skinner: Three-dimensional finite element modeling of bone: Effects of element size. J Biomech Eng. 1992, 14: 483-489.View ArticleGoogle Scholar
- 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. 10.1016/0021-9290(94)90049-3.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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. 10.1016/0021-9290(95)80008-5.View ArticleGoogle Scholar
- 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.Google Scholar
- H Weinans, R Huiskes, HJ Grootenboer: Convergence and uniqueness of adaptive bone remodeling. Trans Orthop Res Soc. 1989, 14: 310-Google Scholar
- 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.View ArticleGoogle Scholar
- 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. 10.1115/1.1503792.View ArticleGoogle Scholar
- HS Kim, STS Al-Hassani: A morphological model of vertebral trabecular bone. J Biomech. 2002, 35: 1101-1114. 10.1016/S0021-9290(02)00053-2.View ArticleGoogle Scholar
- 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. 10.1002/jor.1100150416.View ArticleGoogle Scholar
- 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. 10.1016/S0021-9290(02)00028-3.View ArticleGoogle Scholar
- JA Kanis: Maintenance of Bone Mass. Textbook of osteoporosis. Edited by: Kanis JA. 1996, Cambridge, Blackwell ScienceGoogle Scholar
- MC de Vernejoule, JM Pierre: Chapter 1: New aspects of bone biology. The spectrum of renal osteodystrophy. Edited by: BD Tilman and S Isidro. 2001, Oxford, Oxford Univ PressGoogle Scholar
- D Taylor, JH Kuiper: The prediction of stress fractures using a 'stressed volume' concept. J Orthop Res. 2001, 19: 919-926. 10.1016/S0736-0266(01)00009-2.View ArticleGoogle Scholar
- J White: SAM5211: Introduction to simulation. High Performance Computation for Engineered Systems. 2001, , Singapore-MIT AllianceGoogle Scholar
- TL Arthur Moor, LJ Gibson: Modeling Modulus Reduction in Bovine Trabecular Bone Damaged in Compression. J Biomech Eng. 2001, 123: 613-622. 10.1115/1.1407828.View ArticleGoogle Scholar
- H William, SAT Saul, TV William, PF Brian: Numerical Recipes in C++. The Art of Scientific Computing. 2002, Cambridge, Cambridge University Press, 2ndGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.