A phase field approach to trabecular bone remodeling

We introduce a continuous modeling approach which combines elastic responds of the trabecular bone structure, the concentration of signaling molecules within the bone and a mechanism how this concentration at the bone surface is used for local bone formation and resorption. In an abstract setting bone can be considered as a shape changing structure. For similar problems in materials science phase field approximations have been established as an efficient computational tool. We adapt such an approach for trabecular bone remodeling. It allows for a smooth representation of the trabecular bone structure and drastically reduces computational costs if compared with traditional micro finite element approaches. We demonstrate the advantage of the approach within a minimal model. We quantitatively compare the results with established micro finite element approaches on simple geometries and consider the bone morphology within a bone segment obtained from $\mu$CT data of a sheep vertebra with realistic parameters.

We introduce a continuous modeling approach which combines elastic responds of the trabecular bone structure, the concentration of signaling mole-cules within the bone and a mechanism how this concentration at the bone surface is used for local bone formation and resorption. In an abstract setting bone can be considered as a shape changing structure. For similar problems in materials science phase field approximations have been established as an efficient computational tool. We adapt such an approach for trabecular bone remodeling. It allows for a smooth representation of the trabecular bone structure and drastically reduces computational costs if compared with traditional micro finite element approaches. We demonstrate the advantage of the approach within a minimal model. We quantitatively compare the results with established micro finite element approaches on simple geometries and consider the bone morphology within a bone segment obtained from µCT data of a sheep vertebra with realistic parameters.

INTRODUCTION
Bone undergoes a continuous renewal process, which helps to maintain its mechanical performance and allows for adaptation to changes in mechanical requirements. Different cells are involved in this remodeling process: osteoclasts, which resorb bone, and osteoblasts, which deposit bone. The process is controlled by mechanosensing osteocytes [1][2][3][4] and provides an example of a homeostatic system where external mechanical loads control the bone mass and structure. Various experimental and theoretical studies analyze the remodeling process on different levels of detail [5][6][7][8][9]. We here introduce a continuous modeling approach with the potential to combine different scales in an efficient way. In an abstract setting we consider bone as a shape changing structure, with concentrations of mechanosensing cells within the bone and resorbing and depositing cells on the bone surface. Depending on the surface concentrations the structure is adapted. In contrast to previous modeling approaches, using micro finite element analysis [10][11][12][13][14], we describe the structure implicitly using a time-dependent phase field function. This not only leads to a more accurate model, as the artificial voxel-roughness of the bone surface can be avoided, but also to a drastic reduction of system size and required computing time. Phase field models have been developed to describe shape changing structures, e.g. in solidification [15] or multiphase fluids [16], or have been used as a general numerical tool to solve problems in complex time-evolving geometries [17] and are today well established in various fields. They also have already been used to compute the mechanical properties of trabecular bone structures [18] and have been validated against micro finite element analysis. The required implicit description of the bone structure can be obtained from imaging tools, such as µCT and µMR by standard algorithms. From the voxel representation of the segmented image a smooth surface triangulation can be constructed, from which a signed distance function and the initial phase field function can be computed. The computation of the remodeling process can then be done on a simple cubic domain, which can easily be decomposed for efficient parallel computations. Fig. 1 shows an example of an implicitly described trabecular structure, the phase field representation and details of the phase field function. The paper is organized as follows: In Section we review the current understanding of the bone remodeling process and motivate approximations for our computational approach. Our goal is to introduce a minimal model to demonstrate the potential of the phase field modeling and to point out possible model extensions. We also compare our assumptions with existing micro finite element analysis. In Section we propose our mathematical model, the numerical approach and required preprocessing steps and provide details on the µCT data. In Section we describe results which relate the bone morphology to an applied force. In Section we discuss the results and give conclusions.

MATERIALS AND METHODS
We briefly review the role of the different cells which are involved in the remodeling process and describe the mechanical properties of bone on the level required for our continuous modeling approach.

Osteocytes
Osteocytes form a network throughout the bone matrix by connecting with each other and the lining cells on the bone surface. They sense mechanical loads and transduce the mechanical signal into a chemical response, which is realized by signaling molecules. On the bone surface these molecules orchestrate the recruitment and activity of osteoblasts and osteoclasts, resulting in the adaptation of bone mass and structure. How the osteocytes sense the mechanical loads and coordinate adaptive alterations in bone mass and architecture is not yet completely understood. For a current review see [19]. For the mechanical sensing several mechanisms have been proposed [20]. Within our modeling approach we consider the volumetric compression and the strain energy density as two options to stimulate the osteocytes. The signaling molecules in the bone are then modeled through a diffusion process, with a decay rate and the mechanical stimulus as a source term. This differs from typical micro finite element analysis approaches, where an exponential decay of the signal molecules is assumed and only contributions within a certain distance are accounted for, e.g. [4,7].

Osteoblasts and osteoclasts
A large number of hypotheses have been postulated regarding bone-cell communication and the role played by various receptor-ligand pathways. Various modeling approaches are concerned with the biochemical signaling between active osteoclasts and osteoblasts [21][22][23]. However, they all work on spatial averages and are not yet incorporated in a computational bone remodeling process. We therefore here only use an effective modeling approach by directly considering the concentration of the signaling molecules at the bone surface as an indicator for resorption and deposition. Similar approximations are made in typical micro finite element analysis approaches with various functional forms. The influences of different remodeling rules like linear, step functions or the originally proposed profile by Frost [24] with a lazy zone have been compared in [5]. While the exact form of the remodeling rule is of importance for quantitative predictions, all forms lead to the emergence of trabecular-like patterns.

Trabecular bone structure and mechanical behaviour
Also from a mechanical perspective, trabecular bone is a highly complex material, being anisotropic with different strengths in tension, compression, and shear and with mechanical properties that vary widely across anatomic sites, and with aging and disease [25]. Various of these material properties remain uncertain. However, different experiments have shown a linear behavior [26] and also have related the anisotropy of the mechanical properties to the bone structure [27]. We therefore also for this part of the model consider the simplest possible approach, an isotropic material with prescribed elastic modulus and Poisson ratio but resolve the trabecular bone structure. The same approximations are made in most micro finite element analysis approaches [28][29][30].

Mathematical model
Let Ω bone be the bone structure and Ω a cuboid computational domain of size L = (L x , L y , L z ) such that Ω bone ⊂ Ω. The mechanical material properties of the trabeculae are supposed to be isotropic with Young's modulus E and Poisson's ratio ν. The bone deformation is described by the displacement vector u obeying the partial differential equation with the linear elastic stress tensor where I is the identity matrix and µ, λ are the Lame coefficients Compression is achieved by applying Dirichlet boundary conditions to the normal component of u on ∂Ω, such that u x,y,z = 0 if x, y, z = 0 and u x,y,z =ū x,y,z if x = L x , y = L y , z = L z , while the tangential displacement components are free andū x,y,z are adapted such that the normal forces (F x , F y , F z ) are kept constant. At the remaining boundaries ∂Ω bone we assume no outer forces and set σ · n = 0, where n is the normal to Ω bone . We assume that osteocytes are continuously distributed through Ω bone and that a growth stimulating signal is produced in the osteocytes stimulated by either the volumetric compression or the strain energy density The signal is propagated by a diffusive process with a constant decay rate. Denoting the concentration of the signaling molecule by c, this leads to the equation with boundary condition n · ∇c = 0 on ∂Ω bone , where k is the decay rate and d the diffusion coefficient, both assumed to be constant and S = S vc,sed . Due to the different time scale, if compared with the remodeling process, we consider the stationary solution.
Finally growth is triggered directly according to the concentration c at the trabecular surface. We assume that the growth velocity in normal direction V = V lin,lazy depends linearly on the signal strength, without or with a lazy zone with positive constants α and β and an intermediate range for c where no growth occurs, controlled by the threshold value T . Mathematically the resulting system of equations (1) - (7) is a free boundary problem. According to the proposed normal forces (F x , F y , F z ) on δΩ the displacement u and the concentration c have to be computed in Ω bone . The obtained signal c at δΩ bone is then used to update the time dependent bone structure Ω bone . This non-linear problem has to be iterated until the solution for u, c and Ω bone converge.

Numerical solution
In micro finite element analysis the time dependent bone structure is accounted for by adding and removing voxels to Ω bone . This is not only costly as it requires a high spatial resolution, it also leads to an artificial roughness of the bone surface. Various numerical methods have been proposed to avoid such manipulation of the computational domain. One of the most successful approaches is the phase field or diffuse domain approach [17]. Here Ω bone is described implicitly through a smooth phase field function φ = 0.5(1 − tanh(r/( √ 8 ))) in Ω, with a small parameter determining the width of the diffuse interface and the signed distance function r, with r = 0 at δΩ bone , r > 0 in Ω bone and r < 0 in Ω \ Ω bone . For → 0, φ converges to the characteristic function for Ω bone . Using φ we can now reformulate the problem in the time-independent domain Ω. The diffuse domain approximation of eqs. (1) and (2) reads The boundary condition σ · n = 0 is implicitly included, see [18]. For eq. (5) we obtain which follows along the same arguments and also already includes the boundary condition n · ∇c = 0. Adapting the bone domain with normal velocity V can be realized by solving with a mobility factor γ, see [31]. The terms on the right hand side essentially guarantee the tanh-profile of φ and the left hand side models the transport of φ by V . Using matched asymptotic analysis one can show, see [17], that the system of equations (8) -(10) converge for → 0 to the original problem (1) -(5) with the boundary moved with velocity V and boundary conditions σ · n = 0 and ∇c · n = 0. We can now use standard discretization techniques to solve for u, c and φ in Ω, again iteratively until the system converges. For a computationally efficient treatment adaptive mesh refinement is used. A high resolution within the diffuse interface is required, which needs to be of order and is comparable to the initial voxel scale. However, away from the diffuse interface a coarser mesh can be used, which drastically reduces the computational cost. The efficiency can further be improved by using differently refined meshes for the components [32,33]. We here use a parallel adaptive finite element approach on unstructured meshes which is implemented in the open source toolbox AMDiS [34,35].

Preprocessing
To use the phase field approach for real bone structures requires various preprocessing steps. From the segmented data on the voxel scale, first a surface triangulation is constructed, from which in a second step the signed distance and phase field function can be computed. Fig. 2 shows the three steps. Various approaches are available for achieving these steps. We here use ParaView (http://www.paraview.org) for the generation and smoothing of the surface mesh. Mesh generation has become a standard tool. However, the quality of these meshes in terms of regularity is often very poor. Automatic construction of high quality surface meshes is still not possible for arbitrary surfaces and is an ongoing research topic. The available algorithms in ParaView provide a reasonable compromise between usability and mesh quality. To compute the signed distance function we embed the structure in a cube and adaptively refine the mesh until we obtain the proposed accuracy to resolve the surface. For each node we compute the shortest distance to the surface using raytracing. The approach is implemented in MeshConv (https://gitlab.math.tudresden.de/wir/meshconv) and produces the initial mesh and signed distance function for the computation in AMDiS (https://gitlab.math.tu-dresden.de/wir/amdis). All these software tools are available under an open source license.

Validation
Before we simulate bone remodeling on real bone structures, we consider two simple examples, a cylinder and a cross, for which we compare the results with a micro finite element analysis approach [4,7], and in the cylinder case   also with a semi-analytical solution. We consider both types of stimuli, the volumetric compression and the strain energy density and consider the linear growth law eq. (6).
For a cylinder with radius R and a compression force F z , we obtain for the stimuli S vc = (1 − 2ν) Fz E 1 πR 2 and S sed = 1 2 F 2 z E 1 π 2 R 4 . Both are constant in Ω bone , which allows to solve eq. (5). We obtain c = S k and thus from eq. (6) the growth law V lin = α k S − β from which the evolution of the cylinder can be obtained. The micro finite element analysis is based on [4,7]. The domain Ω is divided into a regular voxel mesh with mesh size h. Each voxel has a continuous value g ∈ [0, 1], where bone is associated with g = 1 and bone marrow corresponds to g = 0. Transient states correspond to voxels which are partly bone. The elastic modulus for a voxel and the stimuli are defined as gE and gS, respectively. The Poisson ratio remains independent of g. Instead of directly solving for the normal velocity V we first compute the change rate of a voxel's g value at ∂Ω bone and the neighboring marrow voxels by solving ∂ t g = αc − β and restricting g ∈ [0, 1]. The normal velocity then follows by multiplying the change rate with the voxel size V = h ∂ t g and applying corrections for the lattice anisotropy. For the numerical solution of eqs. (1) and (5), each voxel with g = 1 is converted to a hex-8 brick element and the resulting systems are solved by standard finite element analysis.
For the comparison we only consider non-dimensional values. The computational domain is the unit cube with L x = L y = L z = 1. The voxel size, as well as the minimum grid spacing in the phase field approach is set to h = 1/128. Time steps were chosen adaptively. Other parameters are d = 1, k = 1/2, F z = 1 for the cylinder and F x = F y = F z = 1 for the cross, as well as β = 1 and α depending on the considered stimulus α sed = 0.36 and α vc = 0.9 for the growth experiment and α sed = 0.04 and α vc = 0.3 for the shrinkage simulations. These parameters lead to a stationary state with R = 0.6/π and R = 0.2/π for growth and shrinkage, respectively. As initial condition we thus use the opposed value, i.e. R = 0.2/π (growth), R = 0.6/π (shrinkage). For the simulation of   Table I according to the guidelines for assessment of bone microstructure [36]. A 0.5 mm aluminium filter was used for beam filtration. Specimens were stored in phosphate buffered saline (PBS) during the scan. Reconstruction of cross sections were computed using the NRecon-Software (Bruker microCT, Kontich, Belgium). Therefore, a gaussian filter (smoothing kernel = 2) was employed. Segmentation, i.e. binarization of the data set, was performed using a locally adaptive thresholding technique (CTAn, Bruker microCT, Kontich, Belgium).

RESULTS
We now apply our model to a segment of a trabecular bone, which is obtained from tomography data of a sheep vertebra. Instead of the linear growth law used for validation we consider eq. (7) with a lazy zone determined by T = 0.2. Volumetric compression was chosen as stimulus and parameters were set to = 0.003, γ = 10, k = 1, d = 0.01, α = 690, β = 1. The elastic properties are chosen as E = 6.829GP a and ν = 0.33, see [7,37,38]. The cubic region Ω is defined by L x = L y = L z = 2.84mm. To account for the dominant loading in z-direction in the animal, the applied force is set twice as high as in the other directions F x = F y = 8N and F z = 16N . The computations are done in parallel on 24 cores with approximately 7.5 Mio degrees of freedom in each time step. The computational time for each setting was approximately 1 hour. Fig. 6 shows the evolution obtained with these parameters. The evolution is shown up to a state for which the main morphological changes have been completed and the concentration c has been mainly equilibrated. A change in morphology is hardly visible, however the computed average microstrain in Ω bone , < | zz | > = 2004, with = 0.5(∇u + (∇u) T ) agrees well with physiological strains measured in bone [39].
To highlight the morphological adaption to anisotropic forces we consider compression forces which are increased by a factor 4 in one direction. Fig. 7 shows the results. The stronger force leads to larger values for c and the adaption of the morphology strongly depending on the direction of the increased force. The change in morphology is clearly visible with regions which are formed and regions which are resorbed. The dependency of the morphology on the direction of the increased force is highlighted in Fig. 8, which shows slices of the bone morphology along the xy-, xz-and yzplane, through the center of the domain. This visualization clearly shows that structures grow in the direction of the enhanced (4-fold) force and thus provide a proof of consistency of the modeling approach. Further validation would require in vivo µCT data of the trabecular bone segment and a calibration of the applied forces. Approaches in this direction can be found in [13,14], which could reproduce changes in bone volume fraction and other global parameters of bone structure but failed to reproduce local bone formation and resorption. One reason for this discrepancy, which is reported in [13,14] are local areas of strong thickening and bone resorption in the experimental images in contrast to more homogeneous layers in their simulations. Fig. 8 clearly shows non-homogeneous morphology changes, see e.g. the level curves for an increased force in x-direction (red curves) in the first and third figure. However, due to lack of in vivo data for the considered bone segment, further validation has to be postponed and we here can only conclude the qualitative consistence of our approach. DISCUSSION We have introduced a continuous modeling approach for bone remodeling which has the potential to combine different scales in an efficient way. In the current form it combines elastic responds of the trabecular bone structure to an applied force, the concentration of signaling molecules within the bone and a mechanism how this concentration at the bone surface is used for local bone formation and resorption. In contrast to previous modeling approaches, using micro finite element analysis, the bone morphology is implicitly described by a time-dependent phase field function. This not only leads to a more accurate model, as the artificial voxel-roughness of the bone surface can be avoided, but also to a drastic reduction of system size and required computing time. The goal of this paper is to provide a minimal model for bone remodeling which demonstrates the advantages of the phase field description. We therefore not only introduce the model, but also quantitatively compare the results with established micro finite element approaches on simple geometries and consider the bone morphology within a segment of 2.84mm 3 obtained from µCT data of a sheep vertebra with realistic parameters. Systematic studies with an enhanced force in one direction clearly demonstrate that the structures grow in the direction of the enhanced force and lead to strong local variations in thickness. These results clearly demonstrate the applicability of the phase field approach. However, quantitative validation would require in vivo µCT data of the trabecular bone segment and a calibration of the applied forces, which are currently not available. But already without such a validation the approach can be used to provide a deeper understanding of the mechanisms underlying bone remodeling. A possible extension considers the incorporation of osteoblast and osteoclast concentrations on the bone surface and their biochemical signaling, which will allow to compute the influence of various signaling molecules on the bone morphology. The ability to predict changes in bone morphology might eventually lead to a better prediction of individual fracture risk in osteoporotic patients or to improved implant materials. The introduced phase field description is ideally suited for this task, as it drastically reduces the computational cost, allows for extensions of additional effects and does only require standard numerical methods, which can be parallelized to deal with significantly larger systems. computing resources provided by JSC within project HDR06 as well as resources provided by ZIH at TU Dresden.