Multiphysics of bone remodeling: A 2D mesoscale activation simulation

Abstract

In this work, we present an evolutive trabecular model for bone remodeling based on a boundary detection algorithm accounting for both biology and applied mechanical forces, known to be an important factor in bone evolution. A finite element (FE) numerical model using the Abaqus/Standard^® software was used with a UMAT subroutine to solve the governing coupled mechanical-biological non-linear differential equations of the bone evolution model. The simulations present cell activation on a simplified trabeculae configuration organization with trabecular thickness of $200 µ m$ . For this activation process, the results confirm that the trabeculae are mainly oriented in the active direction of the principal mechanical stresses and according to the principal applied mechanical load directions. The trabeculae surface activation is clearly identified and can provide understanding of the different bone cell activations in more complex geometries and load conditions.

Keywords

Mechanobiology bone remodeling FEM trabecula multiphysics

1. Introduction

It is well know that bone density distribution evolve as a function of its mechanical environment [1] but the rules governing this evolution are not well understood and are depending mainly on the bone biology. This has been well accepted and developed since the famous Wolff’s Law [2]. Nowadays, detailed understanding of these phenomena has a real interest for the prevision of bone remodeling and bone reconstruction in many orthopedic and surgical applications.

Bone remodeling is a mechanobiological phenomenon occurring at different length scales from the macroscopic to the nanoscopic. Modeling the existing links between these scales is sometimes complex [3] but building precise bone remodeling numerical models requires multiphysics [4] and multiscale [5–7] approaches in coordination with the cells activity since they are the main biological driving force of bone remodeling. In addition, bone is considered as a poroviscoelastic material sometimes grafted with bioresorbable materials when fractured [8–14] with a kinetic evolution depending both on the mechanical properties and on the feeding of biological nutrients [15,16].

At the mesoscopic scale, the trabeculae are important in the remodeling of bone where osteoblasts and osteoclasts are acting on their surfaces. We use here the coordination of the three types of bone cells being the osteoblasts, osteoclasts and osteocytes [17] to propose a finite element model of bone remodeling at the trabecular level solving a set of ordinary differential equations coupling the mechanics and biology [18–25]. According to the mechanical information, the surface of the trabeculae will either be resorbed by osteoclasts or (re)built by osteoblasts [26]. The simulation is used to predict the bone evolving density and microstructure during remodeling.

2. Numerical model development

2.1. Bone remodeling model

We develop a continuous theoretical numerical model where the osteocytes perceive a mechanical information, and translate it into biological signals [18–24] so that for an increasing stimulus, osteoblasts produce bone and on the contrary, for a low stimulus, osteoclasts degrade bone. The developed model is based on a strain energy density function where the elastic deformation energy influences directly the biology of the bone. The biological stimulus S received at any point x of the geometry and time t is given by $S (x, t) = \int_{Ω} U (x_{0}, t) * d (x_{0}, t) * exp (- ‖ x - x_{0} ‖ D^{- 1}) d x_{0} - S_{0} (x, t)$ where U is the strain energy, d is the mechanosensibility of the osteocyte as a function of bone density and D is a normalized distance identified for osteocyte activity. In addition, we define bone density $ρ_{bone}$ evolution through kinematic parameters according to the following rule $\begin{matrix} \frac{\partial ρ_{bone}}{\partial t} = H (ϕ) * \{\begin{matrix} {synth}_{bone} * S & for S > 0 \\ {resorb}_{bone} * S & for S ⩽ 0 \end{matrix} \end{matrix}$ where ${synth}_{bone}$ and ${resorb}_{bone}$ are the coefficients parameterizing the kinetics of resorption and synthesis of bone, and $H (φ)$ illustrates the surface available for reconstruction and resorption. $S_{0} (x, t)$ is the osteocyte activation threshold determining if resorption or formation of bone will occur.

2.2. Trabecular surface search algorithm

Spongy bone is constituted by trabeculae, small beams interconnected and filled with bone marrow. The bone remodeling process begins at the trabeculae surface by the action of osteoclasts receiving the signal transmitted by the osteocytes forming a sealing zone and dissolving bone [25,26]. Then osteoblasts are recruited by the increase of the osteocytes signals sent through the osteocytic canalicular network to the surface where reconstruction occurs. This phenomenon is directly integrated in the bone model remodeling equations acting on the trabecular surface of the bone.

We present here the algorithm, to detect the trabeculae boundary elements (see Fig. 1(a)). This algorithm (developed in the Abaqus UMAT subroutine) enables the localization of the bone cell activation on the trabeculae surface (and therefore will better control the macroscopic bone density evolution). The first step of the algorithm is to create an array that lists in our model the adjacent element of each element (who share two nodes at least). Then we set the existing physical bone material with a given numerical value (1) and empty spaces (void) within the trabeculae with another numerical value (0) as presented on Fig. 1(b) according to: $\begin{matrix} {Flag}_{i} = \{\begin{matrix} 1 & if ρ_{i} > 0.90 \\ 0 & if ρ_{i} ⩽ 0.90 \end{matrix} \end{matrix}$ where ${Flag}_{i}$ represents the identifier (bone or void) and $ρ_{i}$ is the bone density. The geometrical model is then completely identified through a given array of numbers. Each element is compared with its neighbors, and if one of them has a different numerical value, the element is considered included on the surface of the trabecula and consequently becomes an active element, in which remodeling can occur. These interfacial elements are detected according to the rule $\begin{matrix} {IsSurface}_{i} = \{\begin{matrix} 0 & if \sum_{j = neighbors} ρ_{j} = 0 \\ 0 & if \sum_{j = neighbors} ρ_{j} = number of neighbors \\ 1 & otherwise . \end{matrix} \end{matrix}$ Here, ${IsSurface}_{i}$ is the identifier of whether the element number is considered as bone (1) or void (0) and $ρ_{j}$ is the bone density. Thus, trabeculae surfaces are identified at the interfaces between values 0 and 1. During analysis, the marked active elements (construction or resorption at the trabecula surface) are specifically defined (see Fig. 1(c)) to be activated (or not) depending on the calculated stimulus $S (x, t)$ .

Fig. 1.

Theoretical definition of the boundary element search algorithm: (a) General configuration. Red illustrates bone element and blue void element. (b) Configuration after flagging. (c) Active element (marked X).

The constitutive bone remodeling equations are then calculated only for the active area elements so remodeling will occur only on these elements. The numerical model calculates the new bone density for each element at each time increment and the algorithm of boundary element search is used again in order to determine the new active area for the next iteration.

3. Results

A simplified bone microstructure is presented here as an arbitrary non oriented grid. This configuration is loaded on each of the four sides at 45° degrees (Fig. 2, top left corner). Only bone remodeling is allowed within the geometry without any external displacement.

Fig. 2.

Trabeculae kinetics. First and second lines: evolution of the bone density (%) and of the stimulus (J) according to time with a 45° load applied on the 4 sides. Third and fourth lines: zoom on bone density (in %), and drawing of the algorithm evolution.

First and second line of Fig. 2 show the bone density and the stimulus evolution as a function of time (from left to right). Using the defined detection algorithm enables to observe good correlation of the bone density evolution regarding the trabeculae orientations in the applied mechanical load directions. Localization of the cell activation is clearly identified through detection of the active zone when the stimulus is specifically localized in this area.

Cell activity is highlighted on lines three and four of Fig. 2 through the algorithm activation where evolutions are shown as a function of time (from left to right). Bone reconstruction (3rd line) occurs specifically on the trabeculae surface. On the fourth line, the algorithm detects the activation zone and FE model elements are identified through an interpolation of the bone density evolution. This defines the new trabeculae boundaries for the new increment of density calculation.

Using this approach shows the principal difference between a macroscopic evolution model (accounting for an average bone density evolution through the entire geometry) and a more localized evolution (accounting for localized biological effects). In turn, more patient specific problems can be determined accounting for their own biology and materials kinetics (prosthesis, illness, …).

4. Conclusion

Phenomenological models available in the literature provide mathematically optimized bone microstructures as a function of phenomenological parameters not being able to account for patient specific microstructure evolutions. In this paper, an algorithm implemented in an Abaqus/Standard^® UMAT subroutine is presented for known trabecular bone orientation. Results show a reorientation of the trabeculae network in the main principal applied mechanical stress. Both the local mechanical information and cells activation are identified and can provide accurate bone reconstruction specific to each patient issue.

Conflict of interest

The authors have no conflict of interest to report.

References

Becquart

Cruel

Hoc

Sudre

Pemelle

Bizios

Logeart-Avramoglou

Petite

and Bensidoum

, Human mesenchymal stem cell responses to hydrostatic pressure and shear stress, European Cells and Materials31 (2016), 160–173. doi:10.22203/eCM.v031a11.

Wolff

, Das Gesetz der Transformation der Knochen, 1st edn, Hirschwald, Berlin, 1892. Reprint: Pro Business.

Rémond

Ahzi

Baniassadi

and Garmestani

, Applied RVE Reconstruction and Homogenization of Heterogeneous Materials, Wiley–ISTE, 2016. ISBN 978-1-84821-901-4.

Bala

Lefevre

Roux

J.P.

Baron

Lesague

Pithioux

Kaftandjian

and Follet

, Pore network microarchitecture influences human cortical bone elasticity during growth and aging, Journal of the Mechanical Behavior of Biomedical Materials63 (2016), 164–173. doi:10.1016/j.jmbbm.2016.05.018.

Goda

Ganghoffer

J.F.

Czarnecki

Wawruch

and Lewinski

, Optimal internal architecture of femoral bone based on relaxation by homogenization and isotropic material design, Mechanics Research Communications76 (2016), 64–71. doi:10.1016/j.mechrescom.2016.06.007.

Sansalone

Gagliardi

Descelier

Haiat

and Naili

, On the uncertainty propagation in multiscale modeling of cortical bone elasticity, Computer Methods in Biomechanics and Biomedical Engineering18 (2015), 2054–2055. doi:10.1080/10255842.2015.1069619.

Janc

Kaminski

Tarasiuk

Bonnet

A.S.

and Lipinski

, Homogenization of trabecular bone microstructure based on finite element method and micro computed tomography, in: 11th World Congress on Computational Mechanics, Barcelona, Spain, 2014, pp. 1012–1017.

Giorgio

Andreaus

and Madeo

, The influence of different loads on the remodeling process of a bone and bioresorbable material mixture with voids, Continuum Mechanics and Thermodynamics28 (2016), 21–40. doi:10.1007/s00161-014-0397-y.

Giorgio

Andreaus

Dell’Isola

and Lekszycki

, Viscous second gradient porous materials for bones reconstructed with bio-resorbable grafts, 2016. hal-01306433.

10.

Giorgio

Andreaus

Scerrato

and Dell’Isola

, A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials, Biomechanics and modeling in mechanobiology15 (2016), 1325–1343. doi:10.1007/s10237-016-0765-6.

11.

Jaramillo-Isaza

Mazeran

P.E.

El-Kirat

and Ho Ba Tho

M.C.

, Effects of bone density in the time-dependent mechanical properties of human cortical bone by nanoindentation, Computer Methods in Biomechanics and Biomedical Engineering17 (2014), 34–35. doi:10.1080/10255842.2014.931090.

12.

and Lekszycki

, Modeling of an initial stage of bone fracture healing, Continuum Mechanics and Thermodynamics27 (2015), 851–859. doi:10.1007/s00161-014-0380-7.

13.

Andreaus

Giorgio

and Lekszycki

, A 2-D continuum model of a mixture of bone tissue and bio-resorbable material for simulating mass density redistribution under load slowly variable in time, ZAMM Zeitschrift fur Angewandte Mathematik und Mechanik94 (2014), 978–1000. doi:10.1002/zamm.201200182.

14.

Andreaus

Giorgio

and Madeo

, Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids, Zeitschrift fur Angewandte Mathematik und Physik66 (2014), 209–237. doi:10.1007/s00033-014-0403-z.

15.

and Lekszycki

, A novel coupled system of non-local integro-differential equations modelling Young’s modulus evolution, nutrients’ supply and consumption during bone fracture healing, Zeitschrift fur Angewandte Mathematik und Physik67 (2016), 111. doi:10.1007/s00033-016-0708-1.

16.

Bednarczyk

and Lekszycki

, A novel mathematical model for growth of capillaries and nutrient supply with application to prediction of osteophyte onset, Zeitschrift fur Angewandte Mathematik und Physik67 (2016), Art. no. 94.

17.

Bonewald

L.F.

, The amazing osteocyte, Journal of Bone Mineral Research26 (2011), 229–238. doi:10.1002/jbmr.320.

18.

Scala

Spingarn

Rémond

Madeo

and George

, Mechanically-driven bone remodeling simulation: Application to LIPUS treated rat calvarial defects, Mathematics and Mechanics of Solids (2016), 1–13. In press. doi:10.1177/1081286516651473.

19.

Madeo

George

and Rémond

, Second-gradient models accounting for some effects of microstructure on remodelling of bones reconstructed with bioresorbable materials, Computer Methods in Biomechanics and Biomedical Engineering16 (2013), 260–261. doi:10.1080/10255842.2013.815856.

20.

Madeo

George

Lekszycki

Nierenberger

and Rémond

, A second gradient continuum model accounting for some effects of micro-structure on reconstructed bone remodelling, Comptes Rendus Mécanique340 (2012), 575–589. doi:10.1016/j.crme.2012.05.003.

21.

Lekszycki

and Dell’Isola

, A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials, ZAMM Zeitschrift fur Angewandte Mathematik und Mechanik92 (2012), 426–444. doi:10.1002/zamm.201100082.

22.

Madeo

Lekszycki

and dell’Isola

, A continuum model for the bio-mechanical interactions between living tissue and bio-resorbable graft after bone reconstructive surgery, Comptes Rendus – Mecanique339 (2011), 625–640. doi:10.1016/j.crme.2011.07.004.

23.

Lekszycki

, Modelling of bone adaptation based on an optimal response hypothesis, Meccanica37 (2002), 343–354. doi:10.1023/A:1020831519496.

24.

George

Spingarn

Madeo

and Rémond

, Effects of mechanical loading conditions on 3D bone reconstruction: A theoretical numerical study for application to maxillo-facial surgery, in: 9th European Solid Mechanics Conference, Madrid, Spain, 2015.

25.

Wang

and Liu

, A microscale bone remodeling simulation method considering the influence of medicine and the impact of strain on osteoblast cells, Finite Elements in Analysis and Design104 (2015), 16–25. doi:10.1016/j.finel.2015.04.007.

26.

Ruimerman

Hilbers

van Rietbergen

et al., A theoretical framework for strain-related trabecular bone maintenance and adaptation, Journal of Biomechanics38 (2005), 931–941. doi:10.1016/j.jbiomech.2004.03.037.