Interaction of ultrasound waves with bone remodelling: a multiscale computational study Cécile Baron, Vu-Hieu Nguyen, Salah Naïli, Carine Guivier-Curien To cite this version: Cécile Baron, Vu-Hieu Nguyen, Salah Naïli, Carine Guivier-Curien. Interaction of ultrasound waves with bone remodelling: a multiscale computational study. Biomechanics and Modeling in Mechanobi- ology, Springer Verlag, 2020, �10.1007/s10237-020-01306-7�. �hal-02911706� HAL Id: hal-02911706 https://hal.archives-ouvertes.fr/hal-02911706 Submitted on 14 Dec 2020 HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non, lished or not. The documents may come from émanant des établissements d’enseignement et de teaching and research institutions in France or recherche français ou étrangers, des laboratoires abroad, or from public or private research centers. publics ou privés. Noname manuscript No. (will be inserted by the editor) Interaction of ultrasound waves with bone remodelling: a multiscale computational study Cécile Baron · Vu-Hieu Nguyen · Salah Naili · Carine Guivier-Curien Received: date / Accepted: date Abstract Ultrasound stimulation is thought to influ- present work investigates the effect of LIPUS on the cor- ence bone remodeling process. But recently, the effi- tical bone from the tissue to the osteocytes, considering ciency of ultrasound therapy for bone healing has been that the impact of the ultrasound stimulation applied at questioned. Despite an extensive literature describing the tissue scale is related to the mechanical stress exper- the positive effect of ultrasound on bone regeneration - imented by the bone cells. To do that, simulations based cell cultures, animal models, clinical studies - there are on the finite element method are carried out in the com- more and more reviews denouncing the inefficiency of mercial software Comsol Multiphysis to assess the wall clinical devices based on low intensity pulsed ultrasound shear stress levels induced by the LIPUS on the osteo- stimulation (LIPUS) of the bone healing. One of the cytes. Two formulations of the wall shear stress were in- reasons to cause controversy comes from the persistent vestigated based on two IF flow models inside the LCN misunderstanding of the underlying physical and bio- and associated with two different values of the LCN logical mechanisms of ultrasound stimulation of bone permeability. The wall shear stress estimate is very dif- repair. As ultrasonic waves are mechanical waves, the ferent depending on the assumption considered. One process to be studied is the one of the mechanotrans- of these two models provides wall shear stress values duction. Previous studies on the bone mechanotrans- in accordance with previous works published on bone duction have demonstrated the key-role of the osteo- mechanotransduction. This study presents the prelimi- cytes in bone mechanosensing. Osteocytes are bone cells nary results of a computational model that could pro- ubiquitous inside the bone matrix, they are immersed in vide keys to understanding the underlying mechanisms the interstitial fluid (IF) inside the lacuno-canalicular of the LIPUS. network (LCN). They are considered as particularly Keywords Bone remodelling · ultrasound stimula- sensitive to a particular type of mechanical stress: wall tion · osteocyte · poroelasticity · wall shear stress shear stress on osteocytes due to the IF flow in the LCN. Inspired from these findings and observations, the 1 Introduction Cécile Baron Aix-Marseille Univ, CNRS, ISM, Marseille, France Bone is a living tissue that is constantly being remod- Aix-Marseille Univ, APHM, CNRS, ISM, Sainte-Marguerite Hospital, Institute for Locomotion, Department of Or- elled, adapting to its mechanical environment and ca- thopaedics and Traumatology, Marseille, France pable of repair. Bone is one of few tissues that can heal E-mail: cecile.baron@univ-amu.fr without forming a fibrous scar. Trauma is the most ex- Vu-Hieu Nguyen pensive medical condition after heart conditions, cost- Univ Paris-Est, CNRS, MSME, Créteil, France ing the United States 56 billion dollars every year. Of Salah Naili that, 21 billion is used for the treatment of fractures. Univ Paris-Est, CNRS, MSME, Créteil, France For these reasons, the efficacious and expedient treat- Carine Guivier-Curien ment of fractures is of paramount importance to the Aix-Marseille Univ, CNRS, Centrale Marseille, IRPHE, Mar- patient, physician, and healthcare system as a whole seille, France (Buza and Einhorn, 2016). Unfortunately, sometimes 2 Cécile Baron et al. the healing process fails and non-unions or delayed as the dominant factor, in bone mechanotransduction unions occur. Delayed fracture healing and nonunion more than streaming for example (Klein-Nulend et al., are observed in up to 5–10% of all fractures, and can 2005, 2013). Therefore, in this preliminary study we present a challenging clinical scenario for the treating chose to focus on the WSS. The main assumption is physician. The beneficial effect of ultrasound on bone that the bone remodelling is partly triggered by the remodelling was discovered in the 1950s (Corradi and level of fluid-induced wall shear stress especially on the Cozzolino, 1953; Yasuda, 1977). Since the 1980s, sev- process of the osteocytes (Bakker et al., 2001; Anderson eral authors have studied this phenomenon in differ- et al., 2005). ent processes of adaptation of the bone to its envi- ronment: growth (Duarte, 1983), targeted remodeling A lot of numerical studies showed the impact of (Chan et al., 2006) and healing (Schortinghuis et al., the physiological loading on the bone to interpret the 2003). mechanobiology of bone healing (Claes and Heigele, 1999; Lacroix and Prendergast, 2002; Isaksson et al., However, within the recent years, several systematic 2006; González-Torres et al., 2010; Gómez-Benito et al., reviews and meta-analyses (Busse et al., 2014; Schan- 2011; Nguyen et al., 2009, 2010a, 2011). These numeri- delmaier et al., 2017; Griffin, 2016) provide moderate cal models usually considered a cyclic compressive load- to high-certainty evidences against the contribution of ing applied to bone at frequency corresponding to daily USS in bone healing. The controversy on the subject activity (walk, run, jump) and related the mechanical continues and many questions remain unanswered (As- stimulus induced at the cell scale to a biological re- penberg, 2017; Mortazavi et al., 2016). Actually, the sponse in terms of proliferation and differentiation of real issue is that the underlying physico-biological mech- bone cells following the model of Prendergast and col- anisms of the of ultrasound on bone remodeling re- leagues (Prendergast et al., 1997). Ultrasound stimula- main unclear, as reported in several reviews (Claes and tion is a type of mechanical loading and as such may Willie, 2007; Martinez de Albornoz et al., 2011; Padilla trigger reactions similar to physiological loading. How- et al., 2014). One of the objectives of this work is to ever the characteristics of ultrasonic loading are very gain insight into the mechanical impact of ultrasound different from those of physiological loading and this stimulation on bone from tissue to cell and explore its raises a number of questions such as: do the mechan- capacity to trigger bone remodeling. otransduction processes identified for physiological load- ing also act on osteocyte at US frequencies? Are there On one hand, clinical studies, animal models and any others phenomena which can be implied in ultra- cell cultures put on evidence an effect of LIPUS on bone sound stimulation? To understand how ultrasound could repair, but the underlying physico-biological mechanisms affect bone remodelling, a computational model is de- remain unclear and feed an important literature as re- veloped in order to estimate the level of mechanical ported in several paper reviews (Claes and Willie, 2007; stimulus applied on the osteocytes by the mesoscopic Martinez de Albornoz et al., 2011; Padilla et al., 2014). ultrasound stimulation. To our knowledge, it is the first On the other hand, the influence of a mechanical load- numerical study on the multiscale interaction of ultra- ing on bone remodelling has been widely explored by sound with the bone remodelling process. Cowin et al. (1991); Weinbaum et al. (1994); Claes and Heigele (1999); Mullender et al. (2004); Klein-Nulend This paper does the link between biological observa- et al. (2013), for instance. Some of these works have re- tions and physical interpretation integrating both bone vealed the osteocytes as the cornerstone of the mechano- levels of porosity: the vascular network and the LCN sensing and the pilots of the bone mechanotransduc- network. However, the aim is not to infer the biological tion. The osteocytes are ubiquituous bone cells sited mechanisms acting inside the cell, but to concentrate on in the lacuno-canalicular network (LCN), immerged in the characterization of the mechanical forces applied to the interstitial fluid (IF) which fullfills the voids be- the cell via the study of the WSS. To appropriately tween osteocytes and the extra-cellular matrix (ECM). tackle the question, this mechanical model integrates Authors consider that when mechanical loading is ap- two levels of porosity: the vascular network (Havers and plied to bone at macroscale, it induces pressure gradi- Volkman canals - 100 µm) reconstructed from CT scans ents inside the bone, resulting in IF flow around the and the lacuno-canalicular network containing the os- osteocytes which are in turn submitted to several me- teocytes (1 µm) taken into account through the poroe- chanical loadings such as radiation forces, hydrostatic lasticity theory. On the basis of previous works, this pa- pressure or wall shear stress (WSS). The latter one has per proposes to estimate the osteocyte wall shear stress been introduced by Weinbaum et al. (1994) and exper- identified as the relevant stimulus in bone mechanotrans- imentally demonstrated by Klein-Nulend et al. (1995) duction. To do this properly, one of the key parameters Ultrasound and bone remodeling : a computational study 3 is the permeability of the LCN, the value of which re- one is the LCN where the osteocytes are sited sur- mains controversial (Cardoso et al., 2013). rounded by the interstitial fluid. This level of porosity As a consequence, the goal of this paper is focused is constituted of an ubiquitous pore network of lacunae on the mechanical aspect of the bone mechanotrans- (pores ∅ ∼ 15 µm) and canaliculi (pores ∅ ∼ 0.5 µm). duction under ultrasound stimulation. The relevancy The size of these voids are too small to be extracted of the proposed mechanical model is tested looking at from the above-mentioned µCT images, consequently the influence of the LCN permeability value on the wall the bone matrix around the vascular pores is modeled shear stress induced on the osteocytes by the ultrasound as a poroelastic medium using Biot’s theory to describe stimulation. The wall shear stress levels assessed for two the wave propagation in an equivalent medium (Biot, permeability values are compared and discussed in re- 1955, 1957). To be close to the in-vivo configuration, lation to the results of the literature. the model takes into account the soft tissues around After this introduction, the paper is organized as the bone considered as water in first approximation. follows. The materials and methods are presented in section 2. We give in subsection 2.1 the in vivo configu- ration, the geometry description is presented in subsec- tion 2.2, the governing equations are given in subsec- 2.2 Geometry description tion 2.3, the weak formulation is given in subsection 2.4 and the computational model is given in subsection 2.5. For this problem, we consider a two dimensional (2D) We show and discuss the main results on the estima- plane strain model as shown in Fig. 1. In this model, the tion of wall shear stress applied to the osteocyte by poroelastic bone domain is denoted Ω b . The soft tissue, ultrasonic stimulation in section 3. Section 4 finishes the marrow and the fluid-filled lacuna are all consid- the manuscript with the conclusion. ered as fluid and denoted by Ω0f , Ωjf (j = 1, .., N − 1) f and ΩN , respectively. The interfaces between the bone Ω and the fluid domains Ωjf are denoted by Γj (j = b 2 Materials and Methods 0, .., N ). The external surface of the soft tissue domain Ω0f consists of 2 parts denoted by Γ0F (free surface) and 2.1 In vivo configuration Γ0P (imposed by a pressure), respectively. Ultrasound stimulation is mainly dedicated to fracture healing in clinical applications, however, the goal of this preliminary study is to investigate the effect of ultra- sound as a trigger for bone remodelling process and not its direct influence on the healing callus. Consequently, we choose to represent a slice of intact cortical bone which can be above or below the fracture site and still in interaction with ultrasonic waves since the US trans- ducer has a diameter of 1 cm, greater than the fracture gap size (few mm). Considering that the remodelling would be managed by the osteocytes, the bone heal- ing should be initiated from the intact parts of bone surrounding the fracture site. Moreover, the ultrasound stimulation of bone remodelling could be used for skele- tal events other than fracture (localized bone loss, bone Fig. 1: Geometry description metastases, osteolytic tumors) and could also be repre- sented by the model developed in this study. This model is therefore a 2D model considering a bone cross-section perpendicular to the bone axis. Cortical bone is con- In what follows, we denote the coordinates of a point sidered as a biphasic medium (fluid-saturated poroe- x by (x1 , x2 ) and the time by t. A superposed dot stands lastic solid) taking into account two levels of poros- for time derivative, ∇ and ∇· stand for the gradient ity. The first one is the vascular porosity level (pores and divergence operators in 2D space, respectively. The ∅ ∼ 100 µm). The vascular pores geometry is extracted symbol “·” denotes the scalar product and the symbol from µCT images (pixel size = 22.5 µm). The vascu- “:” between tensors of any order denotes their double lar pores are supposed to be full of fluid. The second contraction. 4 Cécile Baron et al. 2.3 Governing equations The constitutive equations for the anisotropic linear poroelastic material are given by: 2.3.1 Equations in the fluid domains Ωjf σ = C : ǫ− αp, (5) The fluid occupying the domain Ωjf is an acoustic fluid 1 − p = ∇· w + α : ǫ, (6) whose mass density and bulk modulus at equilibrium M are denoted by ρfj and Kjf , respectively. By using the where C is the elasticity fourth-order tensor of drained linear acoustic theory and by neglecting the body forces porous material; α, which is a symmetric second-order (other than inertia), the wave equation for the fluid in tensor, is the Biot effective tensor; M is the Biot scalar the domain Ω j may be expressed only in terms of the modulus; ǫ(x, t) is the infinitesimal strain tensor which pressure field pfj (x, t) as follows: is defined as the symmetric part of ∇us . 1 1 p̈fj − ∇2 pfj = 0, ∀x ∈ Ωjf (j = 0, ..., N ). (1) Kjf ρfj 2.3.3 Boundary and interface conditions Note that the index j is used for indicating the number of the fluid domain, then there is no summation over j The surface of the soft tissue consists of 2 parts: Γ0 = in the above equation as well as in the following. Γ0P ∪Γ0F which correspond to the zone excited by pulsed The velocity vector v j (x, t) is relied to the gradient pressure and the free surface, respectively: of pressure field pfj (x, t) following the Euler equation: pf0 (x, t) = P0 g(t), ∀x ∈ Γ0P , (7) ρfj v̇ j + ∇pfj = 0, ∀ x ∈ Ωjf , (j = 0, ..., N ). (2) pf0 ( x, t) = 0, ∀x ∈ Γ0F , (8) 2.3.2 Equations in the anisotropic poroelastic bone (Ω b ) where P0 is the amplitude and g(t) is pulsed time function with a central frequency fc : g(t) = sin(2πfc t). The cortical bone material is assumed to consist of a At interfaces Γj ((j = 0, ..., N )), the continuity of pres- solid skeleton (with mass density ρs ) and a connected sure and stress fields between the poroealstic medium pore network saturated by fluid (with mass density ρf ). and the fluid domains imposes: The Biot’s anisotropic poroelastic model (Biot, 1957; pf = pfj , ) Coussy, 2005; Thompson and Willis, 1991) could be ∀x ∈ Γj (j = 0, ..., N ), (9) used to describe the wave propagation problem in bone σnj = −pfj nj , (Cowin, 2003; Nguyen and Naili, 2012). where nj is the normal unit vector to Γj pointing from Neglecting the body forces (other than inertia), the the porous domain Ω b towards the fluid domain Ωjf (see equations describing the wave propagation within the Fig. 1). bone domain (Ωb ) read: The periosteal interface (Γ0 ) is assumed to be impervi- ous (Goulet et al., 2008), requiring: ∇ · σ = ρ üs + ρf ẅ , (3) −∇p = ρf üs + κ −1 ẇ + b ẅ , (4) w = 0, ∀x ∈ Γ0 . (10) where ρ = φ ρf + (1 − φ) ρs is the mixture density, φ is At the other interfaces Γj (j = 1, ..., N ), open-pore the porosity, σ(x, t) is the total stress tensor and p(x, t) condition is assumed, requiring: is the interstitial fluid pressure in the pores; the vectors of displacement of the solid skeleton and of the fluid are denoted by us (x, t) and uf (x, t), respectively; the ẇ · nj = (v j − u̇s ) · nj , ∀x ∈ Γj (j = 1, ..., N ). (11) vector of relative displacement between the fluid and the solid frame weighted by the porosity is denoted by w = φ(uf − us ); κ is the permeability tensor defined In view of the Euler equation (2), the interface con- by κ = k/µ where k is the intrinsic permeability ten- dition (11) may be rewritten as: sor and µ is the fluid’s dynamic viscosity; the tensor b ! is defined as b = (ρf /φ) a where a is the tortuosity 1 ∇pfj + ẅ + üs · nj = 0, ∀x ∈ Γj (j = 1, ..., N ). tensor. Note that both tensors k and b are symmetric ρfj second-order tensors. (12) Ultrasound and bone remodeling : a computational study 5 2.4 Weak formulation Green-Gauss theorem (Nguyen et al., 2010b): Z Z Z 2.4.1 Weak formulation in the fluid domains Ωjf ũs · ρüs + ũs · ρf ẅ + ∇ũs : σ Ωb Ωb Ωb X Z   + ũs · pfj nj = 0, The pressure field in the fluid occupying the domain Ωjf Γj j is described by Eq. (1). The weak form of this equation may be obtained by multiplying it by a scalar valued ∀ũs ∈ Cad , (16) test function p̃fj (with support in Ωjf ) and then integrat- Z Z Z
w̃ · ρf üs + w̃ · (bẅ) + w̃ · κ−1 ẇ ing over the domain Ωjf . By applying the Green-Gauss Ωb Z Ωb Ωb theorem and taking into account the boundary condi- XZ   + ∇w̃ : (p I) + w̃ · pfj nj = 0, tions (12), the weak formulation of the acoustic wave Ωb j Γj problem in Ωjf reads: ∀w̃ ∈ Cad , (17) where I is the second-order identity tensor. 1 1 Z Z p̃fj p̈fj + ∇p̃fj · ∇pfj 2.5 Computational model Ωjf Kjf Ωjf ρf Z + p̃fj (ẅ + üs ) · nj = 0, ∀p̃fj ∈ Cad , Ultrasound signal. The ultrasound signal emitted from Γj Γ0P (see Fig. 1) is a pulsed signal with a central fre- (13) quency fc =1 MHz, a duty cycle of 20% and a repetition frequency fr =1 kHz. The spatial average-time average acoustic intensity delivered, ISATA , is of 30 mW/cm2 corresponding to an acoustic pressure amplitude where Cad is the admissible function space of the prob- P0 =67 kPa (assuming pressure amplitude is constant lem constituted by the sufficiently differentiable real- over the transducer face). The transducer diameter is valued functions. equal to 2 cm. Note that the nj is pointing toward the interior of the fluid domain Ωjf . Material properties. – Fluid phase: soft tissue, marrow and fluid inside the vascular pores are supposed to be perfect fluids mechanically equivalent to water: speed of sound 2.4.2 Weak formulation in the bone domain Ω b cf =1500 m/s; mass density ρf =1000 kg/m3 ; bulk modulus K=2.25 GPa. In order to resolve the time-dependent problem which is a high frequency problem, the finite element analy- – Poroelastic phase: bone matrix around the vascular sis were based on the (us − w) formulation. The bal- pores is considered as a transverse isotropic poroe- ance equations of linear momentum (3)-(4) are already lastic medium. The LCN porosity is taken from the written in terms of us and w. Accordingly, constitutive literature φ=0.05 (Smit et al., 2002; Lemaire et al., equations (5)-(6) may be restated as: 2012), the fluid inside the pores is considered as wa- ter with a dynamic viscosity µ=10−3 Pa.s−1 . Around the LCN, the extra-cellular matrix (ECM) σ = Cu ǫ + M α ∇ · w, (14) has a mass density ρECM equal to 2 g/cm3 . The p = −M (α : ǫ + ∇ · w) , (15) model is constructed in two dimensional space in a plane perpendicular to the bone axis where the where Cu = C + M α ⊗ α is the undrained elasticity elastic properties of the ECM are assumed to be tensor. (The symbol ⊗ designates the tensorial product isotropic. The elasticity tensor C in this plane is ex- between two tensors.) pressed from the Voigt notation as (Scheiner et al., 2016): Let ũs and w̃ be two vector valued test functions with support in Ω b . The weak form of balance equations   C1111 C1122 0.0 (3)-(4) may be obtained by multiplying them by ũs and C =  C1122 C1111 0.0  , (18) w̃, respectively, integrating over Ω b and applying the 0.0 0.0 C1212 6 Cécile Baron et al. with C1212 = 12 (C1111 − C1122 ). The values used in the velocity of the interstitial fluid relative to the ECM the simulations are C1111 = 22.88 GPa and C1122 = obtained from the Biot’s theory (ẇ). We do this using 10.14 GPa which lead to the Young modulus E = two different models of fluid flow associated with two 16.56 GPa and Poisson ratio ν = 0.308. The Biot’s different permeability values of the LCN. parameters are calculated from this elasticity tensor. The component values of the Biot effective tensor and the Biot scalar modulus are given by α11 = KC-model. In line with the first permeability value stud- α22 = 0.15 and M = 35.6 GPa (for more details see ied (kKC ), the wall shear stress τKC can be firstly as- Rosi et al. (2016)). sessed using Kozeny-Carman theory, through equation Concerning the permeability of this porous medium, (19): two values are investigated. Inspired from the per- Ri µ|ẇ| meability estimation of the vascular porous network τKC = , (19) 2kKC in cortical bone (Zhang et al., 1998; Swan et al., 2003), a first value of LCN permeability is calcu- with Ri is the typical radius of canaliculi in the LCN lated from the Kozeny-Carman relation, assuming (50×10−9 m), µ is the IF dynamic viscosity (10−3 Pa.s−1 ), a Poiseuille flow in a network of cylindrical cap- kKC is the LCN permeability (1.56×10−17 m2 ) and ẇ is illaries aligned in the direction ei of radius Ri = the velocity of the interstitial fluid relative to the ECM 50 × 10−9 m and representing a porosity φ=5%: calculated from the Biot’s theory (see §2.3.2). Equa- kKC = φRi2 /8 = 1.56 × 10−17 m2 . A second one tion (19) assumes that the LCN is modelled as parallel is taken from literature (Smit et al., 2002), estima- straight tubes full of fluid embedded in a solid matrix. tion based on the best fit between finite element predictions and experimental data on canine bone: kS = 2.2 × 10−22 m2 . WT-model. Wall shear stress levels can also be defined from Wang and Tarbell (1995) as: Computation parameters. The convergence of the nu- µ|ẇ| merical results was tested to choose a mesh size (∆x) τWT = √ , (20) kS smaller than the tenth of the wavelength λ in each material. To ensure a relevant sampling of the time- where kS is the LCN permeability equals to 2.2×10−22 m2 dependent phenomena, the solver time stepping ∆t was (Smit et al., 2002). This equation given in Goulet et al. fixed to T/40 where T=1/fc (fc is the central frequency (2008) considers the fluid flow in an annular space through of pulsed time function and ∆t = 2.5 × 10−8 s). Note a net of transverse fibers taking into consideration the that this time step value verifies the calculation stabil- presence of the osteocyte process inside the canalicu- ∆x lus and the glycosaminoglycans (GAG) fibers which are ity condition ∆t < √2V , where Vmax is a maximum max velocity estimated. transverse fibrils which anchor and center the cell pro- Numerical simulations were run with the commercial cess in its canaliculus (see Fig. 2). finite-element software Comsol Multiphysics. 3.1 KC-model vs WT-model 3 Results and Discussion We calculate the moving average of the WSSIF over a As it has been widely reported in litterature, the osteo- signal period (1 µs). Data are smoothed with the “loess” cytes are the cornerstone of bone remodelling. One of method (LOcally Estimated Scatterplot Smoothing), their main characteristics is that they are mechanosen- which is a method using linear regression of smooth sitive cells (see for instance Bonewald (2011)). It has data. The method is altered to assign a zero weight been demonstrated that the dendritic processes of the to data outside six mean absolute standard deviations. osteocytes are more sensitive to mechanical loading than We compare the time evolution of WSSIF (τKC or τWT ) the cell body (Han et al., 2004; Burra et al., 2010) and in over 3 cycles for both LCN permeability values kKC and particular to fluid shear stress (Bakker et al., 2001; An- kS (see Fig. 3). In what follows, the quantities τKC and derson et al., 2005). Inspired from these assumptions, τWT are evaluated in points shown on the bone geom- this study focuses on the IF wall shear stress (WSSIF ) etry. induced by the ultrasound stimulation into the canali- Firstly, the computational model developed in this culi. However, the present poroelastic approach does study provides results on WSSIF induced by ultrasound not allow to directly assess WSSIF acting on individual stimulation at cell scale which can be considered as rea- osteocytes, but we can still calculate the WSSIF from sonable. Comparing both models (KC and WT) used to Ultrasound and bone remodeling : a computational study 7 3.2 Focus on WT-model WSSIF distribution over 10 cycles. The stimulation du- ration is extended to 10 cycles (10 ms) for the WT- model. On Fig. 4, the WSSIF distribution is represented for three different times: 1 ms, 5 ms and 10 ms. The most stimulated areas are located around the endos- teum and the vascular pores interfaces. Globally the WSSIF levels increase with time. To further explore evolution in time of the WSSIF potentially experienced by the osteocytes, τWT is evaluated on 12 points dis- tributed on the bone surface during 10 cycles (see Fig. 5). The curves are obtained in the same procedure as de- scribed in § 3.2. Some points, located near the endos- teum interface and nearby the vascular pores, reach higher τWT values than the six points investigated in Fig. 2: Schematic methodology to calculate WSSIF from the Fig. 3: 0 < τWT < 1.5 Pa, still in the range defined by output parameter ẇ given by the poroelastic theory imple- mented in the numerical model. Two assumptions are inves- Weinbaum et al. (1994). This WSSIF distribution may tigated: the KC-model of the fluid flow inside the canaliculus, raise questions about the possible role of lining cells ignoring the existence of the osteocyte process, into the bone remodelling stimulated by ultrasound. surrounded by the ECM following the Kozeny-Carman rela- Lining cells are bone cells covering over 90% of the tion (left) and the WT-model including the process and the GAG fibers immersed in interstitial fluid to match up to rela- surfaces of adult bone, and which are also recognized tion given by Wang and Tarbell (1995) (right). The osteocytes as mecanosensory cells (Mullender and Huiskes, 1997; is in light gray, the interstitial fluid in dark gray and the GAG Robling and Turner, 2009). They are liable to play a fibers in black; the ECM is hatched. role in bone remodelling triggering, especially since ul- trasonic stimulation transmits compression waves per- pendicular to the bone axis and thus a form of mechan- represent the canalicular flow, a very significant differ- ical stress different from physiological loading, prefer- ence can be observed: the WSSIF levels τKC , range from entially oriented parallel to the bone axis. 0 to 600 Pa i.e. 3 orders of magnitude higher than those of τWT between 0 and 0.6 Pa. It is therefore clear that The obtained values (see Fig. 5) and the spatial het- the way to calculate WSSIF in the porous medium has erogeneity of the WSSIF distribution (see Fig. 4) are in a great influence on the results. In comparison with lit- agreement with the results reported in Fan et al. (2016) erature, only the WT-model provides WSSIF values in on mouse bone with a similar canalicular flow model as agreement with a potential mechanical signal to acti- the WT-model, even if the loading conditions are not vate osteocyte: 0.8<τ <3 Pa (Weinbaum et al., 1994), the same, in particular, the frequency loading (around even if they are close to the lower limit of the inter- 0.5 Hz) which is lower than the pulsed ultrasound fre- val. These results confirm that the KC-model, bundle quencies. of aligned cylindrical pores full of fluid, neglecting the Another factor that can influence the distribution presence of the process, is inappropriate to correctly of the mechanical stimulus is absorption that could represent the interstitial fluid behavior in the LCN, as lead to wave attenuation. The Biot’s poroelastic model it has already been proven by Weinbaum et al. (1994). used in this work to represent the interaction of the They showed that one has to consider the fibers sur- LCN with ultrasonic waves assumes that the solid phase rounding the cell process. The WT-model, firstly devel- is elastic and the fluid is viscous (µ = 10−3 Pa/s). oped to model the flow through muscle cells, has been This absorption is taken into account via the fluid’s already used for modeling the fluid shear stress inside dynamic viscosity in the permeability tensor (α). On the LCN of the cortical bone (Goulet et al., 2008) using the other hand, absorption related to the solid phase of the analogy of the system osteocyte process and GAG the porous medium is neglected in this first approach. fibers with muscle cell and muscle fibers but also in scaffolds designed for bone regeneration (Ouyang et al., Frequency loading dependency. Indeed, the frequency of 2019; Ali and Sen, 2018) including under LIPUS (Wu cyclic loading is regarded as a decisive factor for bone et al., 2015). formation (Weinbaum et al., 1994; Bacabac et al., 2004; According to these observations, only results from Klein-Nulend et al., 2005; González-Torres et al., 2010; the WT-model are analysed further. Gómez-Benito et al., 2011). It is noteworthy that in 8 Cécile Baron et al. Fig. 3: WSSIF estimated in 6 points over 3 cycles (3 ms) for (a) KC-model and (b) WT-model. tion frequency, 1 kHz. Both are much higher than the physiological loading frequencies (1-20 Hz). However, a c >2 Pa it has been demonstrated that a same load applied on bone at mesoscopic scale, at high frequency, could in- duce sufficient stress on the osteocytes to activate a bio- logical response when it would have had no effect at low frequencies (Weinbaum et al., 1994; Han et al., 2004). However the frequency range (0-100 Hz) of these obser- vations remains below the two frequencies considered 0 here. Following Bacabac et al. (2004), one of the de- terminant parameter for osteocytes response is the rate Fig. 4: WSSIF estimated for WT-model at a) 1 ms, b) 5 ms of fluid shear stress (the product of fluid shear stress and c) 10 ms (same color scale). amplitude in Pa and frequency loading in Hz). In their work on cell culture in a flow chamber, they assumed that low-magnitude, high-frequency mechanical stimuli may be as stimulatory as high-amplitude, low-frequency stimuli. In the present work, assuming that the trend remains similar for frequencies beyond 9 Hz (maximum frequency studied in Bacabac et al. (2004)), it would mean that considering the signal frequency of 1 MHz, a fluid shear stress magnitude of 50 µPa could be suf- ficient to activate an osteogenic response, and consid- ering the pulse repetition frequency of 1 kHz, the fluid Fig. 5: WSSIF estimated for WT-model in 12 points over 10 shear stress magnitude reliable to induce osteocyte re- cycles (10 ms). sponse should be around 0.05 Pa. In both cases, τWT is higher than these levels of fluid shear stress. Stimulation duration. One of the limitations of this study is that the maximum stimulation duration considered is 10 ms which is much shorter than the daily treatment time (20 min) defined in the fracture repair protocols (Poolman et al., 2017). Consequently, one of the points to check is that there are no cumulative effects during stimulation cycles. On Fig. 5, it can be seen that val- Fig. 6: WSSIF averaged over each of the ten first cycles. ues obtained on the chosen points are quite stable over time. This trend is confirmed on Fig. 6 which shows the value of τW T averaged per cycle of 1 ms for each con- our case, two frequencies can be considered and can sidered point. Nevertheless, further works are needed to play a role in the bone cell stimulation: the frequency resolve technical issues (simulation time and file size) of the ultrasonic signal, 1 MHz; and the pulse repeti- and to look at longer-time phenomena. Ultrasound and bone remodeling : a computational study 9 4 Conclusion P. Aspenberg. Comment to a BMJ EditorialIs LIPUS the baby in the bathwater? Acta Orthopaedica, 88 This work provides a numerical model of the interaction (1):1, 2017. of low pulsed intensity ultrasound and cortical bone as R. G. Bacabac, T. H. Smit, M. G. Mullender, S. J. a two-level porous medium. The main objective of this Dijcks, J. J. W. A. Van Loon, and J. Klein-Nulend. study is to gain insight into the process of mechan- Nitric oxide production by bone cells is fluid shear otransduction induced by ultrasound stimulation. Con- stress rate dependent. Biochemical and Biophysical sidering the osteocytes as the pilot of the bone remod- Research Communications, 315(4):823–829, 2004. elling under mechanical stimuli, the modelling of the A. D. Bakker, K. Soejima, J. Klein-Nulend, and LCN was essential. The reconstruction of this poros- E. H. Burger. The production of nitric oxide and ity network from images not being available over the prostaglandin E2 by primary bone cells is shear stress entire surface considered, the model uses the theory dependent. Journal of Biomechanics, 34(5):671–677, of poroelasticity to assess the physical quantities rel- 2001. ative to this porosity scale and estimate the effect of M. A. Biot. Theory of elasticity and consolidation for a ultrasonic waves at the cell scale. The results obtained porous anisotropic solid. Journal of Applied Physics, tend to prove the potential of the model detailed in this 26(2):182–185, 1955. study to represent the interaction of ultrasound stimu- M. A. Biot. The elastic coefficients of the theory of lation and cortical bone taking into account two levels consolidation. Journal of Applied Mechanics, 24:594– of porosity. Furthermore, they highlight how the value 601, 1957. of the LCN permeability could affect the WSSIF levels. L. F. Bonewald. The amazing osteocyte. Journal of According to literature, it tends to demonstrate that Bone and Mineral Research, 26(2):229–238, 2011. the KC-model is not relevant to represent the fluid flow S. Burra, D. P. Nicolella, W. L. Francis, C. J. Freitas, into the LCN and that the WT-model seems to be more N. J. Mueschke, K. Poole, and J. X. Jiang. Den- appropriate to evaluate the WSSIF levels potentially ex- dritic processes of osteocytes are mechanotransducers perienced by the osteocytes. that induce the opening of hemichannels. Proceedings Nevertheless, most of the assumptions have been inter- of the National Academy of Sciences of the United polated from previous works on bone mechanotransduc- States of America, 107(31):13648–13653, 2010. tion under physiological loading, more data are needed J. W. Busse, M. Bhandari, T. A. Einhorn, J. D. Heck- on osteocyte activation under higher frequency stimula- man, K-S. Leung, E. Schemitsch, P. Tornetta, S. D. tion to identify the physical and biological processes at Walter, and G. H. Guyatt. Trial to re-evaluate ultra- work in the context of pulsed ultrasound stimulation. sound in the treatment of tibial fractures (TRUST): Further studies are on going to investigate the influence a multicenter randomized pilot study. Trials, 15:206, of US parameters on the mechanical signal received by 2014. osteocytes and explore the effect of the multiscale bone J. A. Buza and T. Einhorn. Bone healing in 2016. Clin- anisotropy through a 3D-model including the absorp- ical Cases in Mineral and Bone Metabolism, 13:101– tion of the solid phase of the bone matrix. Investiga- 105, 2016. tions are under progress to experimentally validate the L. Cardoso, S. P. Fritton, G. Gailani, M. Benalla, and model. S. C. Cowin. Advances in assessment of bone poros- ity, permeability and interstitial fluid flow. Journal of Biomechanics, 46(2):253–265, 2013. C. W. Chan, L. Qin, K. M. Lee, M. Zhang, J. C. Y. Cheng, and Kwok S. Leung. Low intensity pulsed References ultrasound accelerated bone remodeling during con- solidation stage of distraction osteogenesis. Journal D. Ali and S. Sen. Permeability and fluid flow-induced of Orthopaedic Research: Official Publication of the wall shear stress of bone tissue scaffolds: computa- Orthopaedic Research Society, 24(2):263–270, 2006. tional fluid dynamic analysis using newtonian and L. E. Claes and C. A. Heigele. Magnitudes of local stress non-newtonian blood flow models. Computers in Bi- and strain along bony surfaces predict the course and ology and Medicine, 99:201–208, 2018. type of fracture healing. Journal of Biomechanics, 32 E. J. Anderson, S. Kaliyamoorthy, J. Iwan, D. Alexan- (3):255–266, 1999. der, and M. L. Knothe Tate. Nano-microscale models L. E. Claes and B. Willie. The enhancement of bone of periosteocytic flow show differences in stresses im- regeneration by ultrasound. Progress in Biophysics parted to cell body and processes. Annals of Biomed- and Molecular Biology, 93(1-3):384–398, 2007. ical Engineering, 33(1):52–62, 2005. 10 Cécile Baron et al. C. Corradi and A. Cozzolino. Effect of ultrasonics J. Klein-Nulend, R. G. Bacabac, and M. G. Mullender. on the development of osseous callus in fractures. Mechanobiology of bone tissue. Pathologie-Biologie, Archivio Di Ortopedia, 66(1):77–98, 1953. 53(10):576–580, 2005. O. Coussy. Thermoporoelasticity. In Poromechanics, J. Klein-Nulend, A. D. Bakker, R. G. Bacabac, pages 71–112. John Wiley & Sons, Ltd, 2005. A. Vatsa, and S. Weinbaum. Mechanosensation and S. C. Cowin. A recasting of anisotropic poroelasticity in transduction in osteocytes. Bone, 54(2):182–190, matrices of tensor components. Transport in Porous 2013. Media, 50(1):35–56, 2003. D. Lacroix and P. J. Prendergast. A mechano- S. C. Cowin, L. Moss-Salentijn, and M. L. Moss. Candi- regulation model for tissue differentiation during dates for the mechanosensory system in bone. Jour- fracture healing: analysis of gap size and loading. nal of Biomechanical Engineering, 113(2):191–197, Journal of Biomechanics, 35(9):1163–1171, 2002. 1991. T. Lemaire, S. Lemonnier, and S. Naili. On the para- L. R. Duarte. The stimulation of bone growth by doxical determinations of the lacuno-canalicular per- ultrasound. Archives of orthopaedic and traumatic meability of bone. Biomech. Model. Mechanobiol., surgery, 101(3):153–159, 1983. 11:933–946, 2012. L. Fan, S. Pei, X. Lucas L., and L. Wang. A multiscale P. Martinez de Albornoz, A. Khanna, U. G. Longo, 3d finite element analysis of fluid/solute transport in F. Forriol, and N. Maffulli. The evidence of low- mechanically loaded bone. Bone Research, 4:16032, intensity pulsed ultrasound for in vitro, animal and 2016. human fracture healing. British Medical Bulletin, 100 L. A. González-Torres, M. J. Gómez-Benito, (1):39–57, 2011. M. Doblaré, and J. M. García-Aznar. Influence S. M. J. Mortazavi, S. A. R. Mortazavi, and M. Pak- of the frequency of the external mechanical stimulus nahad. Mode & mechanism of low intensity pulsed on bone healing: a computational study. Medical ultrasound (LIPUS) in fracture repair. Ultrasonics, Engineering & Physics, 32(4):363–371, 2010. 71:142, 2016. G. C. Goulet, N. Hamilton, D. Cooper, D. Coombe, M. G. Mullender and R. Huiskes. Osteocytes and D. Tran, R. Martinuzzi, and R. F. Zernicke. Influ- bone lining cells: which are the best candidates for ence of vascular porosity on fluid flow and nutrient mechano-sensors in cancellous bone? Bone, 20(6): transport in loaded cortical bone. Journal of Biome- 527–532, 1997. chanics, 41(10):2169–2175, 2008. M. G. Mullender, A. J. El Haj, Y. Yang, M. A. van X. L. Griffin. Low intensity pulsed ultrasound for frac- Duin, E. H. Burger, and J. Klein-Nulend. Mechan- tures of the tibial shaft. BMJ, 355:i5652, 2016. otransduction of bone cells in vitro: mechanobiology M. J. Gómez-Benito, L. A. González-Torres, E. Reina- of bone tissue. Medical & Biological Engineering & Romo, J. Grasa, B. Seral, and J. M. García-Aznar. Computing, 42(1):14–21, 2004. Influence of high-frequency cyclical stimulation on V.-H. Nguyen and S. Naili. Simulation of ultrasonic the bone fracture-healing process: mathematical and wave propagation in anisotropic poroelastic bone experimental models. Philosophical Transactions of plate using hybrid spectral/finite element method. the Royal Society of London A: Mathematical, Phys- International Journal for Numerical Methods in ical and Engineering Sciences, 369(1954):4278–4294, Biomedical Engineering, 28(8):861–876, 2012. 2011. V.-H. Nguyen, T. Lemaire, and S. Naili. Numerical Y. Han, S. C. Cowin, M. B. Schaffler, and S. Weinbaum. study of deformation-induced fluid flows in periodic Mechanotransduction and strain amplification in os- osteonal matrix under harmonic axial loading. C. R. teocyte cell processes. Proceedings of the National Mecanique, 337:268–276, 2009. Academy of Sciences, 101(47):16689–16694, 2004. V.-H. Nguyen, T. Lemaire, and S. Naili. Poroelastic be- H. Isaksson, W. Wilson, C. C. van Donkelaar, haviour of cortical bone under harmonic axial load- R. Huiskes, and K. Ito. Comparison of biophysical ing: A finite element study at the osteonal scale. Med- stimuli for mechano-regulation of tissue differentia- ical Engineering & Physics, 32:384–390, 2010a. tion during fracture healing. Journal of Biomechan- V.-H. Nguyen, V. Sansalone, and S. Naili. Simulation of ics, 39(8):1507–1516, 2006. ultrasonic wave propagation in anisotropic cancellous J. Klein-Nulend, A. van der Plas, C. M. Semeins, N. E. bone immersed in fluid. Wave Motion, 47(2):117–129, Ajubi, J. A. Frangos, P. J. Nijweide, and E. H. 2010b. Burger. Sensitivity of osteocytes to biomechanical V.-H. Nguyen, T. Lemaire, and S. Naili. Influence of stress in vitro. The FASEB Journal, 9(5):441–445, interstitial bone microcracks on strain-induced fluid 1995. flow. Biomechanics and Modeling in Mechanobiology, Ultrasound and bone remodeling : a computational study 11 10:963–972, 2011. D. M. Wang and J. M. Tarbell. Modeling interstitial P. Ouyang, H. Dong, X. He, X. Cai, Y. Wang, J. Li, flow in an artery wall allows estimation of wall shear H. Li, and Z. Jin. Hydromechanical mechanism be- stress on smooth muscle cells. Journal of Biomechan- hind the effect of pore size of porous titanium scaf- ical Engineering, 117(3):358–363, 1995. folds on osteoblast response and bone ingrowth. Ma- S. Weinbaum, S. C. Cowin, and Yu Zeng. A model for terials & Design, 183:108151, 2019. the excitation of osteocytes by mechanical loading- F. Padilla, R. Puts, L. Vico, and K. Raum. Stimulation induced bone fluid shear stresses. Journal of Biome- of bone repair with ultrasound: a review of the pos- chanics, 27(3):339–360, 1994. sible mechanic effects. Ultrasonics, 54(5):1125–1145, L. Wu, L. Lin, and Y.-X. Qin. Enhancement of cell 2014. ingrowth, proliferation, and early differentiation in a R. Poolman, T. Agoritsas, R. Siemieniuk, I. Harris, three-dimensional silicon carbide scaffold using low- I. Schipper, B. Mollon, M. Smith, A. Albin, S. Nador, intensity pulsed ultrasound. Tissue Engineering. W. Sasges, S. Schandelmaier, L. Lytvyn, T. Kuijpers, Part A, 21(1-2):53–61, 2015. L. Van Beers, M. Verhofstad, and P. Vandvik. Low in- I. Yasuda. The classic: Fundamental aspects of fracture tensity pulsed ultrasound (LIPUS) for bone healing: treatment by Iwao Yasuda, reprinted from J. Kyoto A clinical practice guideline. BMJ, 356:j576, 2017. Med. Soc., 4:395-406, 1953. Clinical Orthopaedics and P. J. Prendergast, R. Huiskes, and K. Søballe. Bio- Related Research, 124:5–8, 1977. physical stimuli on cells during tissue differentiation D. Zhang, S. Weinbaum, and S. C. Cowin. Estimates at implant interfaces. Journal of Biomechanics, 30 of the peak pressures in bone pore water. Journal of (6):539–548, 1997. Biomechanical Engineering, 120(6):697–703, 1998. A. G. Robling and C. H. Turner. Mechanical Signaling for Bone Modeling and Remodeling. Critical reviews in eukaryotic gene expression, 19(4):319–338, 2009. G. Rosi, V.-H. Nguyen, and S. Naili. Numerical in- vestigations of ultrasound wave propagating in long bones using a poroelastic model. Mathematics and Mechanics of Solids, 21(1):119–133, 2016. S. Schandelmaier, A. Kaushal, L. Lytvyn, D. Heels- Ansdell, R. A. C. Siemieniuk, T. Agoritsas, G. H. Guyatt, P. O. Vandvik, R. Couban, B. Mollon, and J. W. Busse. Low intensity pulsed ultrasound for bone healing: systematic review of randomized con- trolled trials. BMJ, 356:j656, 2017. S. Scheiner, P. Pivonka, and C. Hellmich. Poromicrome- chanics reveals that physiological bone strains in- duce osteocyte-stimulating lacunar pressure. Biome- chanics and Modeling in Mechanobiology, 15(1):9–28, 2016. J. Schortinghuis, B. Stegenga, G. M. Raghoebar, and L. G. M. de Bont. Ultrasound stimulation of maxillo- facial bone healing. Critical Reviews in Oral Biology and Medicine, 14(1):63–74, 2003. T. H. Smit, J. M. Huyghe, and S. C. Cowin. Estima- tion of the poroelastic parameters of cortical bone. Journal of Biomechanics, 35(6):829–835, 2002. C. C. Swan, R. S. Lakes, R. A. Brand, and K. J. Stew- art. Micromechanically based poroelastic modeling of fluid flow in haversian bone. Journal of Biome- chanical Engineering, 125(1):25–37, 2003. M. Thompson and J. R. Willis. A reformation of the equations of anisotropic poroelasticity. Journal of Applied Mechanics, 58(3):612–616, 1991.