<figure><img src="image_1.png" /></figure>

HAL
open science

# A thermodynamics framework to describe bone remodeling: a 2D study

Madge Martin, Thibault Lemaire, Guillaume Haiat, Peter Pivonka, Vittorio Sansalone

► To cite this version:

Madge Martin, Thibault Lemaire, Guillaume Haiat, Peter Pivonka, Vittorio Sansalone. A thermodynamics framework to describe bone remodeling: a 2D study. CFM 2017 - 23ème Congrès Français de Mécanique, Aug 2017, Lille, France. hal-03465402

HAL Id: hal-03465402

https://hal.science/hal-03465402

Submitted on 3 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A thermodynamics framework to describe bone remodeling: a 2D study

M. MARTINª, T. LEMAIREª, G. HAIATª, P. PIVONKAᵇ, V. SANSALONEª

a. Laboratoire Modélisation et Simulation Multi Echelle MSME UMR 8208 CNRS
e-mail: madge.martin@u-pec.fr; thibault.lemaire@univ-paris-est.fr; haiat@u-pec.fr;
vittorio.sansalone@u-pec.fr

b. School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology
e-mail: peter.pivonka@qut.edu.au

## Résumé :

Le remodelage osseux est un phénomène complexe traduisant le renouvellement continuel ainsi que l'adaptation du tissu osseux. Cette évolution permanente est d'une importance majeure pour la cationisation et l'ostéo-intégration autour d'un implant. Dans ce travail, le tissu osseux est modélisé comme un matériau orthotrope dont les axes principaux tournent au cours du temps. Le remodelage, exprimé comme une rotation de la microstructure du matériau, est ainsi couplé avec les contraintes macroscopiques. Ce modèle est ici appliqué à un élément macroscopique de tissu osseux afin de suivre l'évolution de l'orientation locale du matériau jusqu'à l'équilibre.

## Abstract :

Bone remodeling is a complex phenomenon that characterizes the lifelong turnover and adaptation of bone tissue. This perpetual evolution is of major importance, as it governs healing as well as osseointegration around implants. In this work, bone tissue is modeled as an orthotropic material whose principal axes rotate over time due to changes of the mechanical environment. The remodeling process is represented in the form of a rotation of the bone microstructure which is linked to the macroscopic loading. This model is applied here to a macroscopic unit of bone tissue in order to follow the evolution of the local orientation of the material until equilibrium.

**Keywords: bone; remodeling; rotation; thermodynamics**

## 1 Introduction

Bone is a living material which is continuously reorganized by bone cells in response to their mechanical and biochemical environment. This process, known as bone remodeling, is of major importance in everyday life, in case of fractures and to allow osseointegration phenomena around implants. Bone remodeling can be described as a stress- and chemistry-driven evolution of the mechanical properties of
bone tissue. The interplay between mechanics and biochemistry, as well as the multiple scales involved in this process, represent serious challenges for the development of realistic bone adaptation models [1-5].

This study describes a novel, thermodynamically sound model of bone remodeling at the tissue scale which, while being mechanistic in nature, is able to account for the above biological processes.

# 2 Materials and Methods

## 2.1 Description of the medium

In the present study, we model bone tissue as a microstructured continuum [5]. We model bone as a body $B$ experiencing remodeling as a rotation of its microstructure. The complete kinematical description of a body point is thus provided by its position in space and a rotation tensor describing the orientation of its microstructure.

When calculating the virtual power of internal and external forces, two virtual velocities are necessary to fully describe the behavior of the body: one relative to the displacement $\hat{v}$ and another one relative to the micro-rotation $\hat{W}$. We assume here the continuum to be a first gradient medium in $\hat{v}$ and zero-gradient in $\hat{W}$. Consequently, the internal and external virtual powers read, respectively:

$$ 
\begin{aligned}
P_i(\hat{v}, \hat{W}) &= \int_B (b_i \cdot \hat{v} - S \cdot \nabla \hat{v} + \hat{T}^i \cdot \hat{W}), \\
P_e(\hat{v}, \hat{W}) &= \int_B (b \cdot \hat{v} + \hat{T}^\circ \cdot \hat{W}) + \int_{\partial B} t_\partial \cdot \hat{v},
\end{aligned} 
\quad (1) $$

where $b_i$ and $b$ represent the inner and outer body forces, respectively, $t_\partial$ is the boundary traction, $S$ is the stress tensor, and $\hat{T}^i$ and $\hat{T}^\circ$ the inner and outer skew-symmetric remodelling couples.

The principle of material frame-indifference to change in observer induces that the power of internal forces is null in the vectorial space of rigid body velocities, which leads classically to $b_i = 0$ and to the symmetry of the tensor $S$. Moreover, the principle of virtual power enunciates that the total virtual power, here the sum of the external and internal forces $P(\hat{v}, \hat{W}) = P_i(\hat{v}, \hat{W}) + P_e(\hat{v}, \hat{W})$ should be null for any admissible virtual velocity $(\hat{v}, \hat{W})$, leading to the following balance laws:

$$ 
\begin{aligned}
\operatorname{div} S + b &= 0 \quad \text{on } \mathcal{B}, & S n_\partial &= t_\partial \quad \text{on } \partial\mathcal{B}, \\
\hat{T}^i + \hat{T}^\circ &= 0 \quad \text{on } \mathcal{B}.
\end{aligned} 
\quad (2) $$

One can note that in the previous set of equations (2), $b$, $t_\partial$ and $\hat{T}^i$ are the fields that depict the impact of the mechanical and biochemical environment on the medium.

A constitutive theory is then developed from the specification of a quadratic strain energy density $\psi$:

$$ \psi = \frac{1}{2} C E : E, \quad (3) $$

where $E$ is the infinitesimal strain tensor and $C$ is the elastic tensor. This latter is assumed to be able to evolve in time, which corresponds to material remodeling. Since we focus on rotary remodeling, the actual elastic tensor $C$ can be obtained through the action of a rotation tensor $R$ (the state variable describing material remodeling) on the initial elastic tensor $C_0$.
A generalized dissipation principle is enforced to obtain thermodynamical restrictions on the constitutive mappings of the inner actions S and $\dot{T}$. No dissipation is assumed to be related to the gross (elastic) deformation, leading to $S = \partial\Psi/\partial E = CE$. (Note that S depends on both E-explicitly-and R-implicitly, via C.) Then, the dissipation principle binds the micro-rotation velocity $W = \dot{R} \cdot R^T$, the displacement velocity $v$ and the material properties via the statement of the positivity of the intrinsic dissipation $D_{int}$:

$$D_{int} = -p_i(v, W) - \dot{\psi} \ge 0, \\ p_i(v, W) = -S \cdot \nabla v + \dot{T} \cdot W = -S \cdot \dot{E} + \dot{T} \cdot (\dot{R} \cdot R^T), \quad (4) \\ \dot{\psi} = \frac{1}{2}(S:E)' = S \cdot \dot{E} - [S, E] : (\dot{R} \cdot R^T).$$

where $p_i(v, W)$ is the power of internal forces per unit volume and the double dot stands for tensor double dot contraction $(A, B) \mapsto A:B = A_{ij}B_{ij}$. Subsequently, we can exhibit from Eq. (4) a structure for the inner remodeling couple $\dot{T}$, involving the dissipation antisymmetric tensor $D$ which interprets as a resistance to remodeling:

$$D_{int} = ([S, E] - \dot{T}) : (\dot{R} \cdot R^T) = D (\dot{R} \cdot R^T) : (\dot{R} \cdot R^T) \ge 0 \quad (5) \\ \dot{T} = [S, E] - D (\dot{R} \cdot R^T),$$

where the superimposed dot stands for time differentiation and the brackets denote the commutator operator $(A, B) \mapsto [A, B] = A.B - B.A$.

## 2.2 Local behaviour

We study here the evolution of the material principal axes as a passive remodeling, which implies that the outer remodeling couple $\dot{T}$ is null. We assume here that dissipation due to remodeling does not depend on the strain E, the microrotation R nor the remodeling velocity $\dot{R}$. Therefore, at a given position, $D$ is constant and equals $D = D_0$. Hence, the remodeling is coupled to local stress and strain by the following equation:

$$D_0 (\dot{R} \cdot R^T) = [S, E]. \quad (6)$$

In 2D, the rotation tensor is parametrized by a scalar (the angle $\alpha$) and the vectorial space of fourth-order tensors is $Skw \otimes Skw = \text{Vect}(\ast \otimes \ast)$ where $\ast = e_2 \otimes e_1 - e_1 \otimes e_2$. Hence, $D_0$ is uniquely described by a coefficient $d_0$.

It follows that, with given initial material properties, the rotation rate $\dot{\alpha}$ only depends on the dissipation $d_0$, the strain E and the orientation $\alpha$. Remodeling equilibrium is achieved when material properties no longer evolve, corresponding to a stationary state of the rotation.

It is worth noting that this model predicts the principal axes of the strain and stress tensors to be locally collinear at the remodeling equilibrium [6]. Thus, a physically sound condition for remodeling equilibrium [7] is recovered without any ad-hoc assumption.

## 2.3 Finite element model

The model presented here describes the 2D rotation of bone principal axes. In order to account for rotary remodeling, we modeled cortical bone as an orthotropic material [8]. The values of the material elastic properties are displayed in Table 1. We monitored the evolution of a square piece of bone subjected
<table><caption>Table 1: Two-dimensional material properties</caption><thead><tr><th>C<sub>11</sub></th><th>C<sub>22</sub></th><th>C<sub>12</sub></th><th>C<sub>66</sub></th><th>d<sub>0</sub></th></tr></thead><tbody><tr><td>30 GPa</td><td>20.85 GPa</td><td>11.49 GPa</td><td>12.44 GPa</td><td>1.0 MPa.s</td></tr></tbody></table>

to a boundary load under plane strain with a finite element analysis. In order to both account for the local response to the local stress field and update the material orientation, two overlapping meshes were necessary.

A coarse mesh defined the rotation of the material axes on specified points $i \in \mathbb{I}$. From these points, a linear interpolation enabled the creation of an orientation field throughout the material. A finer mesh specified the elements of the finite element analysis.

An initial orientation was given to the macroscopic continuum structure. At each step $n \in \mathbb{N}$, the orientation field established from the coarse mesh $\{\alpha_{i,n}\}_{i \in \mathbb{I}}$ leads to a finite element analysis built from the new material properties. The calculated strain field $E_n$ obtained from the analysis is then used to update the orientation of the material axes on the coarse mesh, according to the following equation deriving from Eq. (6).

$$ \forall n \in \mathbb{N}, \forall i \in \mathbb{I}, \alpha_{i,n+1} = \alpha_{i,n} + f(\alpha_{i,n}, d_0, E_n) \Delta t. \quad (7) $$

Note that the function $f$ in Eq. (7) qualifies the approximation method used. In this study, fourth-order Runge-Kutta approximation was preferred to Euler-Cauchy method in order to optimize the computational cost.

Hence, the orientation of the material on the coarse mesh interpolation points is determined by the previous orientation, the strain field $E$ and a fixed time step $\Delta t$. Therefore, the time step represents the response time of the material to a new mechanical environment. On the other hand, the finite element study was chosen to be stationary. Hence, the change in material properties instantaneously impacts the strain field.

We considered the following criterion to characterize equilibrium:

$$ \sup_{i \in \mathbb{I}} |\dot{\alpha}_{i,n_{eq}}| \le \dot{\alpha}_{crit}. \quad (8) $$

In other words, the material has reached equilibrium if and only if none of the coarse mesh interpolation points $i \in \mathbb{I}$ are evolving. This criterion is consistent with the definition of remodeling equilibrium established in Section 2.2.

# 3 Results and Discussion

## 3.1 Homogeneous boundary conditions

We imposed displacements to achieve uniform stress/strain conditions to the microstructured continuum. The system underwent rotary remodeling until achieving equilibrium. Remodeling equilibrium states were found to correspond to strain energy minima. In Figure 1, we subjected the material to a constant longitudinal strain $E = e_1 \otimes e_1$ while varying the initial orientation of the medium. When subjecting the material to a longitudinal strain, the medium tends to minimize the strain energy to reach equilibrium. Two equilibrium states are observed in the present loading configuration: one corresponds to an alignment of the elastic principal axes with the loading direction ($\alpha_{0,\infty} = 0$), and the other corresponds
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the strain energy of a 2D body under constant longitudinal strain, with varying initial states {$α_{0,0} = k \frac{\pi}{12}, k \in \{0, 1, 2, 3, 4, 5, 6\}$}.</figcaption></figure>

to a rotation towards the normal to that direction ($α_{0,∞} = \frac{\pi}{2}$).

Furthermore, the minimization of strain energy in these imposed strain configurations reads as a minimization of potential energy. Therefore, we have here a sensible condition to the achievement of equilibrium.

## 3.2 Finite element study

<figure><img src="image_3.png" /><figcaption>Figure 2: 2D body under boundary traction. White lines represent material orientation before (bottom) and after (top) remodeling. Left: strain energy density (color map, arbitrary units). Right: Von Mises stress (color map, arbitrary units).</figcaption></figure>

Finite element analyses depicted the elastic principal axes of the material in response to prevailing stresses until a steady state was reached. Figure 2 displays the adaptation of the material when subjected to a boundary load. When changing the orientation, the change of material properties leads to an
evolution of the strain energy as well as von Mises stress.

The prediction of the rotation of the principal axes of the material in simple loading configurations is
consistent with the superimposed boundary conditions. These results confirm the ability of the model
to simulate the material response to non-uniform stress configurations.

4 Conclusion

A novel, thermodynamically sound model of bone remodeling was proposed. A toy problem was used to investigate the rotary remodeling of bone at the tissue scale. Elastic principal axes of bone tissue were observed to rotate according to the superimposed mechanical forces, which is interpreted as a reorien- tation of the material microstructure. Preliminary results obtained in 2D are promising and show the potential of this approach. Model predictions for in vivo biomechanical loading configurations need to be further tested. Suitable experimental data will be identified. The particular case of tissue surrounding an implant will be also studied. Our model can also integrate the mechanobiological phenomena regu- lating bone remodeling. However, this would require a reliable description of the biochemical stimuli of bone remodeling. This matter is out of the scope of this paper and will be addressed in future works.

5 Acknowledgement

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 682001, project ERC Consolidator Grant 2015 BoneImplant).
References

[1] Y. K. Kim, Y. Kameo, S. Tanaka, T. Adachi, Capturing microscopic features of bone remodeling into a macroscopic model based on biological rationales of bone adaptation. *Biomech Model Mechanobiol*, in press, 2017.

[2] I. Goda, J.F. Ganghoffer, G. Maurice, Combined bone internal and external remodeling based on Eshelby stress. *Int J Sol Struct*, Vol. **94**, 138–157, 2016.

[3] M. Colloca, R. Blanchard, C. Hellmich, K. Ito, B. van Rietbergen, A multiscale analytical approach for bone remodeling simulations: linking scales from collagen to trabeculae. *Bone*, Vol. **64**, 303–313, 2014.

[4] I. Scala, C. Spingarn, Y. Rémond, A. Madeo, D. George. Mechanically-driven bone remodeling simulation: Application to LIPUS treated rat calvarial defects. Mathematics and Mechanics of Solids, 1081286516651473. The elastic anisotropy of bone. *Mth Mech Sol*, 1081286516651473, 2016.

[5] A. Di Carlo, S. Naili, S. Quiligotti. Sur le remodelage des tissus osseux anisotropes. *C. R. Mecanique*, Vol. **334**, 651-661, 2006.

[6] V. Sansalone, S. Naili, A. Di Carlo. On the rotary remodeling equilibrium of bone. *Comput Methods Biomech Biomed Engin*, Vol. **14(s1)**, 203–204, 2011.

[7] S. Cowin, Wolff's law of trabecular architecture at remodeling equilibrium. *J Biomech Eng*, Vol. **108**, 83–88, 1986.

[8] J.L. Katz, A. Meunier. The elastic anisotropy of bone. *J Biomechanics*, Vol. **20**, 1063–1070, 1987.