# Mechanically-driven bone remodeling simulation: application to LIPUS treated rat calvarial defects

Ilaria Scala, Camille Spingarn, Yves Rémond, Angela Madeo, Daniel Georges

► To cite this version:

Ilaria Scala, Camille Spingarn, Yves Rémond, Angela Madeo, Daniel Georges. Mechanically-driven bone remodeling simulation: application to LIPUS treated rat calvarial defects. Mathematics and Mechanics of Solids, 2016, 10.1177/1081286516651473. hal-02077445

HAL Id: hal-02077445

https://hal.science/hal-02077445

Submitted on 19 Mar 2020

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mechanically-driven bone remodeling simulation: Application to LIPUS treated rat calvarial defects

**Ilaria Scala**

MSME, Université Paris-Est, Paris, France

LGCIE, Université de Lyon-INSA, Lyon, France

**Camille Spingarn**

ICube, Université de Strasbourg-CNRS, Strasbourg, France

**Yves Rémond**

ICube, Université de Strasbourg-CNRS, Strasbourg, France

**Angela Madeo**

LGCIE, Université de Lyon-INSA, Lyon, France

**Daniel George**

ICube, Université de Strasbourg-CNRS, Strasbourg, France

**Abstract**

In this paper we numerically simulate the phenomenon of bone growth in bone defects as driven by external mechanical excitation. Bone growth is accounted for through a continuum model that allows simulation of the filling of a defect. The influence of the model boundary conditions is also discussed. Two and three dimensional simulations are presented, explicitly showing the bone regeneration process inside the cavity on a weekly basis. Numerical results are qualitatively compared with literature experimental data from a rat calvarial defect exposed to low-intensity pulsed ultrasound. The obtained results show trend correlations with the targeted phenomenological observations and allow us to perform a first evaluation of the proposed model parameters to be optimized for clinically relevant situations, even if a systematic experimental campaign is still needed to precisely identify the bio-mechanical parameters involved.

**Keywords**

Mechanical-biological coupling, growth resorption process, continuum mixture model, LIPUS driven bone growth

<footnote>**Corresponding author:**
Angela Madeo, LGCIE, Université de Lyon-INSA, 20 Av. Albert Einstein, 69100 Villeurbanne Cedex, France.
Email: angela.madeo@insa-lyon.fr</footnote>
# I. Introduction

Bone remodeling has been extensively studied since the publication of the well-known Wolff's Law [1], which enabled better understanding of bone architecture organization under applied mechanical load and was later redefined by Frost [2, 3]. Following Wolff and Frost, Carter et al. [4, 5] studied the effects of mechanical load history on the developed mechanical energy for bone architecture reorganization. For various types of loading scenarios, they showed similar mathematical formulations of bone density evolutions as a function of effective applied stress, fatigue damage accumulation and strain energy density, and these could be applied for functional adaptation and bone growth. Then, Carter and Beaupré [6] highlighted the influence of the mechanical factors in bone growth and showed the specific effects of mechanical loads on the biological reaction leading to better understanding of bone mechanobiology. In this paper, following the methods presented in [7, 8], we propose to use a continuum mixture model to describe the physical phenomenon of bone growth inside a cavity as driven by mechanical excitation. In particular, we study the influence of low-intensity pulsed ultrasound (LIPUS) on the ability to globally reproduce the gradual filling of a rat calvarial defect.

LIPUS is employed more and more today to promote the bone repair process in fractured bones. The usual medical practice consists of a LIPUS treatment of the fractured bone which lasts some minutes per day. Typical applications of LIPUS entail a 5–20 min/daily treatment with frequencies ranging from almost 1–10 MHz and with intensities of the order of 10–100 mW/cm². LIPUS is currently applied transcutaneously, although recent experimental studies have proven the efficacy of a trans-osseous application. In fact, Malizos et al. [9] reported on the first trans-osseous application of ultrasound on a sheep tibial osteotomy model and demonstrated that the treated bones became significantly stiffer and stronger with higher bone mineral density compared to the untreated tibiae.

Over the last few years different studies (see [10–13]) have been performed to unveil the efficacy of LIPUS to promote bone growth, such as those devoted to establishing the actual positive effect of LIPUS treatment on the enhanced filling of defects. Duarte [14] observed the evolution of small holes with diameters of 1.5 mm in the cortex of the femur in rabbits subjected daily to a 15 min ultrasound treatment for 2 weeks. He found that the LIPUS treatment stimulated the callus formation inside the holes compared with the contralateral non-treated holes. Lavandier et al. [15] performed similar experiments on calvarial bone defects, but they focused their attention on the study of the effect of the different LIPUS characteristics on bone regeneration. They came to the conclusion that the main parameter influencing bone growth is the *intensity* of the applied load, rather than its frequency. This seems reasonable if one thinks that the bone growth process has characteristic time-scales of the order of weeks, while the characteristic time associated with the load application is of the order of seconds. This follows work from Carter et al. [4], who studied the influence of daily loading histories onto the evolution of trabecular bone density. Although cyclic stresses are influencing bone mass, as long as there is no observable bone density variation, it can be assumed that as long as bone is experiencing neither a net loss or gain of bone the apparent density will be exposed to a constant daily stimulus. Hence, in our model, assuming that frequency is not relevant in our case study [4, 15] and that ultra-sound intensity leads to small strains, we will suppose in our modeling phase that the LIPUS signal can be approximated with a static force representing the averaged value of the applied dynamic load and deforming the considered bone specimen elastically. These hypotheses remain sensible, as the time scale of the applied cyclic mechanical load is many orders of magnitude smaller than the time scale of the biological processes. In fact, if the deformation response to particular applied loads occurs on time scales of the order of seconds or minutes, the characteristic time to observe bone growth due to cellular activities is of the order of weeks [2, 3, 16–18].

Hasuike et al. [19] demonstrated the bone regenerative effect of LIPUS treatment on rat 2.7 mm non-critical calvarial defects, as confirmed with in vivo micro-CT. The rats were divided into a LIPUS rat group (exposed for 20 min daily to a typical LIPUS application) and a control rat group (not exposed to mechanical excitation) and the beneficial effect of LIPUS treatment on bone regeneration was clearly demonstrated.

In the remainder of this paper, we will use the experimental results of [19] in order to show that our model allows us to properly simulate bone growth inside a cavity in accordance with available experimental evidence.

The development of the theoretical model used here to simulate bone growth inside a cavity is based on some mixture models originally conceived to describe the phenomenon of bone remodeling in the presence of an artificial bio-resorbable graft. In such models [7, 8, 20–26] the first constituent of the mixture represents natural bone tissue, while the second represents the bio-resorbable material which is filling the void defect.

In this paper, we use a similar model to simulate the phenomenon of bone remodeling when the defect is not filled with an artificial material after surgery. To reproduce this situation, we will assume in our numerical simulations that the geometrical domain occupied by the defect is filled with a “fictive material” with almost
vanishing density and Young's modulus. In this way, the subregion of the considered geometrical domain, which is supposed to be filled with such “fictive material”, can be fully interpreted as a void region. When performing numerical simulations, we introduce a “fictive material” inside the cavity to avoid the problem related to the modeling of a moving material interface. This reduces the difficulty of conceiving a model with an additive surface growth to a problem of volume growth within the adjacent domain. The latter situation is evidently more easily treatable in a finite element framework. Considering an extremely soft material in the defect region is actually much more realistic than considering a hole in which no matter at all is present. Indeed, when producing a defect, some soft tissues are always present in the cavity from the earliest phases after surgery. Such soft tissues are then mineralized to be transformed in cartilage and finally in bone [27]. In the results section, we present some two-dimensional (2D) numerical simulations showing that the proposed model is able to catch the main features of bone regeneration in a calvarial bone defect as driven by mechanical excitation. In addition, we test different types of loading conditions to try to establish some comparisons with the effect of the application of a LIPUS treatment.

Finally, we will introduce some preliminary three-dimensional (3D) results in order to show that the proposed model can be further optimized in order to simulate fully realistic growth patterns in three dimensional structures.

## 2. Low-intensity pulsed ultrasounds (LIPUS) as a source of mechanical excitation for bone remodeling

Ultrasound is a form of energy that is transmitted through biological tissues in the form of acoustic waves, and is today widely used in medicine as a diagnostic and therapeutic tool. Indeed, low-intensity pulsed ultrasounds can be seen as an effective source of mechanical excitation, which is able to promote the process of bone remodeling. This is, in turn, known to be strongly dependent on the application of external loads [1, 28, 29]. Based on the evidence of mechanical stimulation affecting the evolution of healthy bones, different studies have recently focused their attention on the possibility of improving the process of bone regeneration via the use of LIPUS in the presence of cavities or defects. However, although there is a common agreement that LIPUS pressure waves can cause biochemical events at the cellular levels that promote bone formation, the precise bio-physical mechanisms at the root of such behavior have not yet been fully clarified. In particular, it is not clear whether such LIPUS-induced remodeling improvement is caused by the energetic content of the traveling waves [15] or by their frequency [30, 31].

In this paper, following [15], we suppose that the main mechanical factor that may influence the remodeling process is the intensity rather than the frequency of the traveling wave. We are in some sense assuming that the activation of the cells is mainly due to the intensity of the load that is sensed by the osteocytes rather than to the frequency of the applied signal. In other words, we are assuming that there exists a sort of load threshold which needs to be exceeded in order to activate the cells which are deputed to the formation of new bone tissue. More precisely, we are considering a situation in which when the intensity of the applied load is lower than a given threshold value, the associated deformation of the bone is low; this means that the osteoblasts do not need to create new bone tissue to correctly sustain the externally applied load. Osteoclasts, on the other hand, can resorb what can be interpreted as an excess of bone tissue, which can be more conveniently used where the applied loads are more intense. When the intensity of the applied load lies on the threshold, a sort of bio-mechanical equilibrium takes place in which neither creation nor resorption of bone tissue is needed to effectively sustain the applied external loads. When the intensity of the applied load is higher than the given threshold, then the associated deformation of the bone is high and, as a result, osteoblasts start forming new bone tissue in order to render the bone more compact. On the basis of such remarks it is reasonable to assume that the LIPUS treatment can be thought to be optimized by suitably choosing, at fixed frequency, the intensity of the energy carried by the traveling waves in order to maximize bone regeneration. Following our previous hypothesis of different time scales for the mechanical and biological phenomena, (see Lavandier et al. [15]), we will study the effects of an “average” static load and observe the evolution of the system on a weekly basis in order to follow the evolution of the bone mass density.

## 3. Theoretical model

In this section we recall the main peculiarities of the continuum mixture model proposed in [7, 8] and we show how this model can be suitably adapted to the particular case of bone growth inside a defect.
On the basis of the experimental results obtained by [15], we suppose in this paper that the effect of the frequency of the applied load on the process of remodeling can be, in a first instance, neglected, and that the main parameter which intervenes in such remodeling processes is indeed the intensity of the traveling wave, i.e. the energetic content of the wave itself. We will hence introduce a static model neglecting inertial effects.

### 3.1. Continuum mixture-model describing the mechanical behavior of a natural-bone/fictive-material mixture

We describe the deformation of a two-solids continuum mixture by introducing a Lagrangian (or reference) configuration $B_L \subset \mathbb{R}^3$ and a suitably regular kinematical field $\chi(\mathbf{X}, t)$ that associates to any material point $\mathbf{X} \in B_L$ its current position $\mathbf{x}$ at time $t$. The kinematics of the system is completed by introducing two Lagrangian mass densities $\rho_b(\mathbf{X}, t)$ and $\rho_m(\mathbf{X}, t)$, which represent the Lagrangian mass density of the natural bone tissue and of the artificial material respectively. Such Lagrangian mass densities are supposed to evolve in time to allow for the possibility of describing the effect of the underlying cellular activity. The image of the function $\chi$ gives, at any instant $t$, the current shape of the body, also called the Eulerian configuration of the mixture. Since we will use it in the following, we also introduce the displacement field $\mathbf{u}(\mathbf{X}, t) := \chi(\mathbf{X}, t) - \mathbf{X}$, the tensor $\mathbf{F} := \nabla\chi$ and the Green–Lagrange deformation tensor $\boldsymbol{\epsilon} := (\mathbf{F}^T \cdot \mathbf{F} - \mathbf{I})/2$. We consider a quasi-static case, i.e. a mechanical problem in which no kinetic energy is considered. In this case, we can write the power of internal forces as the first variation of a suitable action functional $\mathcal{A}$ as follows:

$$ \mathcal{P}^{\text{int}} = \delta\mathcal{A} = -\delta \int_{B_L} U(\boldsymbol{\epsilon}) d\mathbf{X}, \qquad (1) $$

where $U(\boldsymbol{\epsilon})$ is the bulk strain energy density. We stress the fact that expression (1) for the power of internal forces is valid only in the quasi-stationary case, when the kinetic energy does not intervene in the description of the mechanical problem. When the inertia terms are relevant, then the kinetic energy density must be added to the integrand in equation (1) and an integration in time must also be performed. The expression of the action functional presented here is of the same type as the one considered in [8], in which a simple one dimensional application was targeted.

In this paper, assuming material isotropy due to the small thickness of rat calvaria mainly composed of cortical bone, we consider the following constitutive form for $U$ as studied in linear elasticity (the Einstein notation of the sum over repeated indices is used if not otherwise specified):

$$ U(\boldsymbol{\epsilon}, \rho_b, \rho_m) = \frac{1}{2} \frac{E(\rho_b, \rho_m) \nu}{(1+\nu)(1-2\nu)} \epsilon_{ii}^2 + \frac{1}{2} \frac{E(\rho_b, \rho_m)}{(1+\nu)} \epsilon_{ij} \epsilon_{ij} $$

where $E$ and $\nu$ are the Young's modulus and the Poisson coefficient of the bone/bio-material mixture, respectively. The Young's modulus is assumed to be dependent on the densities of both bone and bio-material according to the following power law:

$$ E(\rho_b, \rho_m) = \left( E_b \left( \frac{\rho_b}{\rho_{\text{max}}} \right)^{\beta} + E_m \left( \frac{\rho_m}{\rho_{\text{max}}} \right)^{\beta} \right), \qquad (2) $$

where $E_b$ and $E_m$ are characteristic constant moduli, $\beta$ is a suitable constant exponent and $\rho_{\text{max}}$ is the maximum possible density that the bone-bio-material mixture can reach during the remodeling process. A reasonable value for this maximum density can be fixed to be approximately $\rho_{\text{max}} = 2000 \text{ kg/m}^3$, which is the density of compact bone. As will be explained, since the “fictive material” must indeed represent a void cavity, its mass density $\rho_m$ and Young's modulus $E_m$ will be considered to be almost vanishing.

We explicitly remark that by treating the void region as a “fictive material”, we are drastically improving the simplicity of the proposed model. In fact, a model which accounts for the possibility of introducing a real void region should also be able to describe a surface growth inside the cavity. This would result in the need to set up a model which accounts for the description of the motion of a deformable body with evolving material interface; this is a real theoretical challenge. Moreover, simulating the existence of a void region by setting its density to be very small allows for a really easy and effective numerical implementation.
As for the expression of the power of external forces, we assume that forces per unit area can act on the external boundary of the body (or on some portion of it), while we assume that no external volume forces are present. More precisely, we consider the following expression for the power of external forces:

$$ \mathcal{P}^{\text{ext}} = \int_{\partial B_L} f_i^{\text{ext}} \delta u_i \quad d\mathbf{X}. \qquad (3) $$

Also for the power of external forces, we recall that expression (3) is valid in the quasi-static case.

The mechanical governing equations in weak form can be directly expressed by imposing the validity of the principle of virtual powers $\mathcal{P}^{\text{int}} = \mathcal{P}^{\text{ext}}$ with suitable kinematical or natural boundary conditions.

## 3.2. Bio-mechanical coupling: Stimulus and growth-resorption laws

The bio-mechanical coupling is assured by the biological stimulus $S$ defined on the Lagrangian configuration $B_L$ as

$$ S(\mathbf{X}, t) = \left( \int_{B_L} U(\mathbf{X}_0, t) d(\mathbf{X}_0, t) \exp(-\|\chi(\mathbf{X}) - \chi(\mathbf{X}_0)\|^2 / D^2) d\mathbf{X}_0 \right) - S_0(\mathbf{X}, t), \qquad (4) $$

where $S_0$ is a threshold of the stimulus discriminating between resorption and synthesis, according to the hypothesis of the existence of a load threshold made in Section 2, and $d$ is a function which allows us to specify the concentration of osteocytes at any point of the domain. The stimulus $S$ can be thought of as a bio-mechanical signal associated to the osteocyte activity, whose intensity is directly proportional to measured deformation and to the density of osteocytes.

Since osteocytes can be present only if natural bone is present (they are formed from osteoblasts when they are completely surrounded by natural bone tissues), we suppose that the concentration of osteocytes is proportional to the mass density of natural bone tissue at any point of the considered domain. More particularly, we consider the following expression for the concentration of osteocytes:

$$ d = \eta \frac{\rho_b}{\rho_{\max}}, \quad \eta \in [0, 1], $$

where $\eta$ is a suitable constant which allows us to account for a higher or lower concentration of osteocytes at a given point of the mixture. From now on we will also make the hypothesis of small deformations, so that we will consider that in equation (4) $\|\chi(\mathbf{X}) - \chi(\mathbf{X}_0)\| \approx \|\mathbf{X} - \mathbf{X}_0\|$.

The Lagrangian mass densities $\rho_b$ and $\rho_m$ of bone and bio-material are assumed to vary in time according to the following ordinary differential equations:

$$ \frac{\partial (\rho_b / \rho_{\max})}{\partial t} = A_b(S)H(\varphi), \qquad \frac{\partial (\rho_m / \rho_{\max})}{\partial t} = A_m(S)H(\varphi), \qquad (5) $$

where

$$ H(\varphi) = k \varphi(1 - \varphi), \quad \varphi = 1 - \theta \frac{(\rho_b + \rho_m)}{\rho_{\max}}, \quad \theta \in [0, 1], $$

$$ A_b(S) = \begin{cases} s_b S & \text{for } S > 0 \\ r_b S & \text{for } S < 0 \end{cases}, \quad A_m(S) = \begin{cases} 0 & \text{for } S > 0 \\ r_m S & \text{for } S < 0 \end{cases}. $$

In these formulas, we introduce some suitable bio-mechanical parameter that is used in our model in order to simulate the activity of the underlying cells. In particular, the value of the parameter $k$ weights the influence that porosity has on the overall remodeling process; such a parameter thus indirectly accounts for the porosity patterns of the considered bone. On the other hand, the parameters $r_b$, $r_m$ and $s_b$ account respectively for the resorption rate for bone and bio-material and for the synthesis rate of natural bone tissue.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Regenerating process in a rat model without LIPUS treatment measured every week over 1 month from initial configuration (from Hasuike et al. 2011) [19].</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Regenerating process in a rat model with LIPUS treatment measured every week over 1 month from initial configuration (from Hasuike et al. 2011) [19].</figcaption></figure>

# 4. A 2D case simulating growth of bone tissue in the presence of a bone cavity.

Different authors (see [32–35]) have been interested in the experimental study of bone growth in the presence of bone defects. A common protocol is to produce a circular-shaped bone defect in the calvarium of a rat (or of another control species), thus observing in vivo, by an imaging procedure, the regeneration of bone in the weeks following surgery. In this section, we show how our model can be used to describe such a regenerative phenomenon.

In Figure 1, we show the evolution of bone growth in the defect for a rat which has not been treated by LIPUS. A slow bone regeneration process is observed, which can be related to a basic cellular activity triggered by mechanical excitations that are compatible with everyday activity such as chewing (or even with the simple presence of weight forces). Figures 1 and 2 show the evolution of bone growth around the defect as observed during a period of approximately 1 month for the control group (not treated with ultrasound) and for the group treated with LIPUS. From the results observed by Hasuike (2011) [19, 33], it turns out that the LIPUS treatment undoubtedly has a beneficial effect on the process of bone healing. In the following we will show how our theoretical model can be applied to reproduce this behavior, after suitable calibration of the constitutive parameters.

## 4.1. Numerical simulations of regeneration in the presence of bone defect

In order to describe the physical situation in which a hole is produced inside bone tissue, the continuum mixture model proposed in Section 3 can be used. We consider the geometry reproduced in Figure 3: our specimen is composed of a rectangular-shaped bone sample 1 cm × 1 cm, with a circular-shaped hole in the center with a diameter of 2 mm. Since no bio-material is present inside the cavity, we will assume that no resorption of considered “fictive material” occurs by setting $r_m = 0$ (see Table 1). Furthermore, it was assumed that the thickness *d* of the specimen was equal to 1 mm.

These hypotheses allow us to adopt the theoretical model introduced in the previous section to describe bone growth inside a hole. The model was implemented within the COMSOL Multiphysics(R) software (see Figure
<table><caption>**Table 1.** Values of the parameters used for the numerical simulation reproducing bone growth inside a hole.</caption><thead><tr><th>E<sub>b</sub> (GPa)</th><th>E<sub>m</sub> (GPa)</th><th>β</th><th>k</th><th>θ</th><th>S<sub>0</sub>(Nm)</th><th>D (mm)</th><th>η</th></tr></thead><tbody><tr><td>20</td><td>0</td><td>1.6</td><td>4</td><td>1</td><td>10<sup>-8</sup></td><td>1</td><td>1</td></tr><tr><td>ρ<sub>b</sub> / ρ<sub>max</sub> at t = 0</td><td>ρ<sub>m</sub> / ρ<sub>max</sub> at t = 0</td><td></td><td></td><td>s<sub>b</sub> (s/m<sup>5</sup>)</td><td>r<sub>b</sub> (s/m<sup>5</sup>)</td><td></td><td>r<sub>m</sub> (s/m<sup>5</sup>)</td></tr><tr><td>0.8</td><td></td><td>10<sup>-5</sup></td><td></td><td>10<sup>5</sup></td><td>2 × 10<sup>4</sup></td><td></td><td>0</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Geometry of the considered bone defect and different types of applied boundary displacements.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Initial distribution of bone mass density: 80% of maximum bone density $\rho_{max}$ outside the hole; (b) adopted mesh used with triangular-shaped elements: curvature factor of 0.4, range of elements size $2 \times 10^{-5}/9 \times 10^{-4}$, maximal growth rate of the elements 1.4.</figcaption></figure>

4). Considering that mesh density directly influences the results, we paid a great deal of attention to this choice. The characteristic length of the biological parameters (particularly for the *D* parameter which is associated with the range of action of osteocytes) covers several mesh elements for a better numerical interpolation. Table 1 summarizes the values of the parameters which have been used to perform the numerical simulations for this case.

In Figures 5 to 10 we show the time evolution of bone density and of the associated bio-mechanical stimulus that follows the application of three types of displacement boundary conditions to the considered specimen, as shown in Figure 3. The order of magnitude of the applied displacements was chosen in order to be compatible with those displacements which are reasonably caused by a LIPUS treatment. In these simulations, the comparable volume force with the studied situation was $F_v = -10^8$ N/m$^2$.

It can be seen that, after the application of the external displacement, bone growth can be observed on a weekly basis. If the observation period is increased, then complete replacement can be observed as shown in Figure 5. Experimental studies (see e.g. [15, 19]) suggest that the effect of external mechanical excitation on regeneration of bone defects can be crucial, above all for the case of “critical size defects” i.e. of defects that would not heal without external mechanical excitation. The model proposed in the present paper allows us to directly account for this effect of mechanical loading on bone growth and, for the boundary conditions shown in Figure 3(a), an almost complete replacement is predicted after 10 weeks (see Figure 5). Similar observations
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Time evolution of tissue growth under the application of the boundary conditions shown in Figure 3(a).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Time evolution of biological stimulus under the boundary condition shown in Figure 3(a).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time evolution of tissue growth under the application of the boundary condition shown in Figure 3(b).</figcaption></figure>

can be given in the other cases of displacements. In particular, the model represented in Figure 3(b) leads to the simulation shown in Figure 7, where an almost complete replacement is predicted after 9 weeks, with a prescribed displacement even smaller than that of case (a) (i.e. 0.00005 mm). In addition, the application of the asymmetric prescribed displacement (0.0001 · y mm) illustrated in Figure 3(c) implies, as shown in Figure 9, an almost complete replacement is predicted after 40 weeks. If we focus on the stimulus (see Figures 6, 8 and 10), we can appreciate the propagation of the osteocyte signal that moves toward the critical size defect, its intensity increasing with the increasing of the observation period. In particular, a more intense red zone corresponds to a higher value of the stimulus and vice versa for the blue region. It is worth remarking that the imposed displacement is very small and corresponds to strains of the order of $10^{-5}$–$10^{-6}$, which are definitely compatible with characteristic displacements due to LIPUS application. When comparing the simulations shown in Figure 5 to the experimental observations in Figure 1, it can be seen that the qualitative behavior of remodeling is completely caught by the proposed model. Some slightly asymmetric growth patterns are present in the experimental case, which could be due to different imperfections, such as the non-perfect circular shape of the initial hole, asymmetry of the applied load, etc. For this reason the qualitative regenerative pattern associated with an asymmetric load shown in Figure 9 appears to be closer to the real situation. The simulated regeneration is nevertheless slower due to the intensity of the externally applied load. The process could be accelerated by increasing the intensity of the applied asymmetric load.

With particular reference to the evolution of the biological stimulus, we can infer that bone growth is always triggered in those regions where the stimulus is much higher. This is sensible if we think that a high stimulus is associated with a high deformation, which induces osteoblasts to synthesize new bone tissue more easily.

Considering the last loading case depicted in Figure 3(c), we note that the corresponding bone growth represented in Figure 9 arrives on larger time-spans. The process would of course be accelerated by increasing the intensity of the applied displacement.
<figure><img src="image_9.png" /><figcaption>Figure 8. Time evolution of biological stimulus under the boundary condition shown in Figure 3(b).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Time evolution of tissue growth under the application of the boundary condition shown in Figure 3(c).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Time evolution of biological stimulus under the boundary condition shown in Figure 3(c).</figcaption></figure>

## 5. Three-dimensional case

In order to show that the introduced continuum model is actually able to capture the main features of bone growth in real living systems, we present in this section some preliminary results concerning bone growth in a three dimensional structure. The geometry and adopted loading conditions are presented in Figure 11. Such a situation mimics the creation of a circular hole in the rat calvaria. The boundary conditions are then set to be completely fixed on the external radius R (corresponding to the radius up to which stresses are applied) and the small hole of radius r (corresponding to the drilled hole). These boundary conditions define the circular plate with fixed external radius and flexure on the top surface (applied pressure).

<figure><img src="image_12.png" /><figcaption>Figure 11. Geometry and load distribution for the 3D case: R = 0.5 cm, r = 1.5 mm, h = 1 mm, F = 45 kPa.</figcaption></figure>

With the applied force F = 45 kPa, we explicitly remark that, with reference to the geometry of the considered specimen given in Figure 11, a measure of the energy per unit area that is furnished to the system can be given by

$$|F| R = 22.5 \frac{\text{mJ}}{\text{cm}^2}.$$

Assuming that the mechanical phenomena occur instantaneously when compared to the biological ones, such a quantity is directly comparable to an “averaged power” of the LIPUS (usually called LIPUS intensity) which, as we mentioned before, is of the order of 10–100 mW/cm².
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Three dimensional time evolution of tissue growth inside a cavity for $D = 1$ mm.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Three dimensional time evolution of the biological stimulus for $D = 1$ mm.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Three dimensional time evolution of tissue growth inside a cavity for $D = 0.4$ mm.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Three dimensional time evolution of the biological stimulus for $D = 0.4$ mm.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Three dimensional time evolution of tissue growth inside a cavity for $D = 0.2$ mm.</figcaption></figure>

All the results presented in this section are relative to a 3D specimen whose geometry and loading conditions are given in Figure 11. The values of the constitutive parameters are the same as those chosen for the 2D simulations and given in Table 1, except for the parameter $D$ for which different values were tested in order to show the sensitivity of the model to such a crucial parameter. Indeed, $D$ represents a characteristic length corresponding to a measure of the range of action of osteocytes. For all of the 2D simulations, as well as for the first 3D simulation depicted in Figures 12 and 13, we chose $D = 1$ mm (equivalent to the thickness of the model). The reconstruction occurs all around the hole with symmetrical distribution through the thickness. This is due to the fact that the parameter $D$ and the thickness of the specimen are comparable.

The bone growth and corresponding biological stimulus are shown on Figures 14 and 15 for a value of the parameter $D = 0.4$ mm and similarly on Figure 16 and 17 for a value of the parameter $D = 0.2$ mm. It can be inferred that when the range of action of osteocytes is smaller, the process of bone growth becomes slower and the observed growth patterns appear to be different to the first case. Since we imposed a general flexure of the plate on the top surface, this is where the bone reconstruction occurs first, progressing further down the thickness of the specimen towards the bottom surface. Changing the value of $D$ from 0.4 mm to 0.2 mm creates a thinner bone reconstruction area, which depends directly on the mechanical stimulus developed.
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Three dimensional time evolution of the biological stimulus for $D = 0.2$ mm.</figcaption></figure>

Results obtained using the 3D model show a clear influence of the $D$ parameter corresponding to the “osteocyte influence radius”. In the case of $D = 1$ mm (corresponding to the model/specimen thickness), we can observe a homogeneous reconstruction with a slightly higher influence in the center of the plate. This is to be expected since the osteocyte radius of influence will distribute the bone kinetic reconstruction homogeneously through the thickness. For the two other cases, we have a different reconstruction through the thickness. In the case of $D = 0.4$ mm (corresponding to about half of the plate thickness), the bone reconstruction occurs mainly in the top surface initially, since our loading condition is a bending state with load applied on the top surface. Reconstruction goes down through the thickness as the analysis progresses, and full reconstruction is observed after about 120 weeks. In the final case for which $D = 0.2$ mm, we clearly have a smaller influence on the cell radius and bone reconstruction, which initially develops on the top surface near the hole edge and progresses towards the center. Even after 120 weeks (similar to $D = 0.4$ mm), full reconstruction is not complete. This is mainly due to the correlation existing between cell radius activity, bone kinetics and time. Even if time plays an important role in the different studied cases linked with other material parameters, we show that the geometry of bone reconstruction kinetics is different for each case. Further experimental studies targeted at establishing a reliable characteristic length for the range of action of osteocytes are needed in order to calibrate the parameter $D$ of our model.

# 6. Conclusions

Nowadays both regenerative medicine and tissue engineering are focused on the development of new therapeutic approaches which may improve the regeneration of living organs and tissues, like bones, affected by trauma or degenerative diseases. The development of theoretical and numerical tools which are able to simulate the process of bone healing as driven by mechanical excitation would be of great help for the choice of the optimal LIPUS treatment to be delivered to the patient, in order to reach optimal healing in the most efficient way.

In this paper, we use a continuum model that is able to reproduce the natural behavior of bone tissue regeneration and growth inside a defect, as influenced by LIPUS treatment. We assume that the real dynamical force exerted by a LIPUS treatment is replaced with an averaged static force. This is rather sensible if we think that the time-scale of the intervening mechanical phenomena associated with the LIPUS application is much smaller than the time-scale of bone growth. The introduced model is then used to simulate the phenomenon of bone growth inside a 2D cavity for different types of externally applied loads. Finally some results concerning a first 3D case are also reported in order to show the generality and potential of the introduced model.

Even if extended experimental campaigns accompanied by a systematic optimization of the parameters of the introduced model are still needed in order to render the proposed tools useful for medical uses, the obtained results are encouraging and in agreement with experimental observation. The obtained patterns of bone growth are seen to be close to the experimental ones for applied load values, which are compatible with the intensities of typical LIPUS treatments.

From the point of view of the mechanical modeling, a possible extension of this work could be to include the effects due to the microstructure of the bone. In this context, strain gradient generalized continuum models are proven to be a possible modeling choice [7, 36, 37] for introducing characteristic lengths and size and anisotropy effects related to the microstructural characteristics of the material.

Moreover, higher grade continuum models also allow for the use of enhanced boundary conditions, useful to represent connections between materials with microstructure [38–40].

It could also be interesting to investigate the effect of osteoporosis on the process of remodeling. To this end, classical [41, 42] or strain gradient [43–45] poroelastic models could be a suitable modeling choice.
Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

[1] Wolff, J. *Das Gesetz der Transformation der Knochen*. Berlin: Verlag von August Hirschwald, 1892.

[2] Frost, HM. Skeletal structural adaptations to mechanical usage (SATMU): 1. Redefining Wolff's law: The bone modeling problem. *Anat Rec* 1990; 226: 403–413.

[3] Frost, HM. Wolff's law and bone's structural adaptations to mechanical usage: An overview for clinicians. *Angle Orthod* 1994; 64(3): 175–188.

[4] Carter, DR, Fyhrie, DP, and Whalen, RT. Trabecular bone density and loading history: Regulation of connective tissue biology by mechanical energy. *J Biomech* 1987; 20(8): 785–794.

[5] Carter, DR, Orr, TE, and Fyhrie, DP. Relationship between loading history and femoral cancellous bone architecture. *J Biomech* 1989; 22(3): 231–244.

[6] Carter, DR, Van der Meulen, MCH and Beaupré, GS. Mechanical factors in bone growth and development. *Bone* 1996; 18(1) 5S–10S.

[7] Madeo, A, George, D, Lekszycki, T, et al. A second gradient continuum model accounting for some effects of micro-structure on reconstructed bone remodelling. *CR Mecanique* 2012; 339(10): 625–640.

[8] Madeo, A, Lekszycki, T, and dell'Isola, F. A continuum model for the bio-mechanical interactions between living tissue and bio-resorbable graft after bone reconstructive surgery. *CR Mecanique* 2011; 339(10): 625–640.

[9] Malizos, KN, Hantes, ME, Protopappas, V, et al. Low-intensity pulsed ultrasound for bone healing: An overview. *Injury* 2006; 37S: S56–S62.

[10] Azuma, Y, Ito, M, Harada, Y, et al. Low-intensity pulsed ultrasound accelerates rat femoral fracture healing by acting on the various cellular reactions in the fracture callus. *J Bone Miner Res* 2001; 16(4): 671–680.

[11] Schortinghuis, J, Ruben, JL, Raghoebar, GM, et al. Ultrasound to stimulate mandibular bone defect healing: A placebo-controlled single-blind study in rats. *J Oral Maxil Surg* 2004; 62(2): 194–201.

[12] Yang, KH and Park, SJ. Stimulation of fracture healing in a canine ulna full-defect model by low intensity pulsed ultrasound. *Yonsei Med J* 2001; 42(5): 503–508.

[13] Yoshida, A, Sasaki, H, Furuya, Y, et al. Effect of low-intensity pulsed ultrasound on bone-healing process in murine low-turnover osteoporosis model. *J Hard Tissue Biol* 2013; 22(3): 301–310.

[14] Duarte, LR. The stimulation of bone growth by ultrasound. *Arch Orthop Traum Su* 1983; 101: 153–159.

[15] Lavandier, B, Gleizal, A, and Béra, JC. Experimental assessment of calvarial bone defect re-ossification stimulation using low-intensity pulsed ultrasound. *Ultrasound Med Biol* 2009; 35(4): 585–594.

[16] Beaupré, GS, Orr, TE, and Carter, DR. An approach for time-dependent bone modeling and remodeling - Application: A preliminary remodeling simulation. *J Orthop Res* 1990; 8: 662–670.

[17] Beaupré, GS, Stevens, SS, and Carter, DR. Mechanobiology in the development, maintenance, and degeneration of articular cartilage. *J Rehabil Res Dev* 2000; 37(2): 142–152.

[18] Carter, DR, and Wong, M. Mechanical stresses in joint morphogenesis and maintenance. In: *Biomechanics of diarthrodial joints*. New York: Springer-Verlag, 1990, pp.157–174.

[19] Hasuike, A, Sato, S, Udagawa, A, et al. In vivo bone regenerative effect of low-intensity pulsed ultrasound in rat calvarial defects. *Oral Surg Oral Med O* 2011; 111: e12–e20.

[20] Andreaus, U, Giorgio, I, and Madeo, A. Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids. *Z Angew Math Phys* 2014; 66(1): 209–237.

[21] Andreaus, U, Giorgio, I, and Lekszycki, T. A 2-D continuum model of a mixture bone tissue and bioresorbable material for simulating mass density redistribution under load slowly variable in time. *ZAMM-Z Angew Math Me* 2014; 94(12): 978–1000

[22] Giorgio, I, Andreaus, U, Lekszycki, T, et al. The influence of different geometries of matrix/scaffold on the remodeling process of a bone and bioresorbable material mixture with voids. *Math Mech Solids* 2015. DOI: 10.1177/1081286515616052.

[23] Giorgio, I, Andreaus, U, Scerrato, D, et al. A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials. *Biomech Model Mechan* 2016. DOI: 10.1007/s10237-016-0765-6.

[24] Lekszycki, T. Functional adaptation of bone as an optimal control problem. *J Theor Appl Mech* 2005; 43(3): 555–574.

[25] Lekszycki, T, and dell'Isola, F. A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials. *ZAMM-Z Angew Math Me* 2012; 92(6): 426–444.

[26] Lu, Y, and Lekszycki, T. Modeling of an initial stage of bone fracture healing. *Continuum Mech Therm* 2015; 27: 851–859.

[27] Carter, DR, and Beaupré, GS. *Skeletal function and form: Mechanobiology of skeletal development, aging and regeneration*. Cambridge: Cambridge University Press, 2007.

[28] Mullender, MG, and Huiskes, R. Proposal for the regulatory mechanism of Wolff's Law. *J Orthop Res* 1995; 13: 503–512.

[29] Wolff, J. *The law of bone remodelling*. Berlin: Springer-Verlag, 1986.

[30] Gomez-Benito, MJ, Gonzales-Torres, LA, Reina-Romo, E, et al. Influence of high-frequency cyclical stimulation on the bone fracture-healing process: Mathematical and experimental models. *Philos T R Soc A* 2011; 369(1954): 4278–4294.
[31] Gonzales-Torres, LA, Gomez-Benito, MJ, Doblaré, M, et al. Influence of the frequency of the external mechanical stimulus on bone healing: A computational study. *Med Eng Phys* 2010; 32(4): 363-371.

[32] Fındık, Y, and Baykul, T. Effect of low-intensity pulsed ultrasound on autogenous bone graft healing. *Oral Surg Oral Med O* 2014; 117(3): e225-e260.

[33] Hasuike, A, Sato, S, Makino, N, et al. The potency of low-intensity pulsed ultrasound in a rat calvarial guided bone-regeneration model. *J Hard Tissue Biol* 2011; 20(3): 217-224.

[34] Logeart-Avramoglou, D, Anagnostou, F, Bizios, R, et al. Engineering bone: Challenges and obstacles. *J Cell Mol Med* 2005; 9(1): 72-84.

[35] Petite, H, Viateau, V, Bensaïd, W, et al. Tissue-engineered bone regeneration. *Nat Biotechnol* 2000; 18(9): 959-963.

[36] Auffray, N, dell'Isola, F, Eremeyev, VA, et al. Analytical continuum mechanics à la Hamilton-Piola least action principle for second gradient continua and capillary fluids. *Math Mech Solids* 2013; 20(4): 375-417.

[37] Auffray, N, Dirrenberger, J and Rosi, G. A complete description of bi-dimensional anisotropic strain-gradient elasticity. *Int J Solids Struct* 2015; 69: 195-206.

[38] dell'Isola, F, Madeo, A, and Placidi, L. Linear plane wave propagation and normal transmission and reflection at discontinuity surfaces in second gradient 3D continua. *ZAMM-Z Angew Math Me* 2011; 92(1): 52-71.

[39] Placidi, L, Rosi, G, Giorgio, I, et al. Reflection and transmission of plane waves at surfaces carrying material properties and embedded in second-gradient materials. *Math Mech Solids* 2013; 19(5): 555-578.

[40] Neff, P, Ghiba, ID, Madeo A, et al. A unifying perspective: The relaxed linear micromorphic continuum. *Continuum Mech Therm* 2014; 26(5): 639-681.

[41] Biot, M. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range. *J Acoust Soc Am* 1956; 28(2): 168-178.

[42] dell'Isola, F, Madeo, A, and Seppecher, P. Boundary conditions at fluid-permeable interfaces in porous media: A variational approach. *Int J Solids Struct* 2009; 46(17): 3150-3164.

[43] Rosi, G, Nguyen, V-H, and Naili, S. Numerical investigations of ultrasound wave propagating in long bones using a poroelastic model. *Math Mech Solids* 2015; 21(1): 119-133.

[44] Rosi, G, Madeo, A, and Guyader, J. Switch between fast and slow Biot compression waves induced by 'second gradient microstructure' at material discontinuity surfaces in porous media. *Int J Solids Struct* 2013; 50(10): 1721-1746.

[45] Sciarra, G, dell'Isola, F, and Coussy, O. Second gradient poromechanics. *Int J Solids Struct* 2007; 44(20): 6607-6629.