## Abstract

In bone functional adaptation by remodelling, osteocytes in the lacuno-canalicular system are believed to play important roles in the mechanosensory system. Under dynamic loading, bone matrix deformation generates an interstitial fluid flow in the lacuno-canalicular system; this flow induces shear stress on the osteocytic process membrane that is known to stimulate the osteocytes. In this sense, the osteocytes behave as mechanosensors and deliver mechanical information to neighbouring cells through the intercellular communication network. In this study, bone remodelling is assumed to be regulated by the mechanical signals collected by the osteocytes. From the viewpoint of multi-scale biomechanics, we propose a mathematical model of trabecular bone remodelling that takes into account the osteocytic mechanosensory network system. Based on this model, a computational simulation of trabecular bone remodelling was conducted for a single trabecula under cyclic uniaxial loading, demonstrating functional adaptation to the applied mechanical loading as a load-bearing construct.

## 1. Introduction

It is known that bone has the ability to adapt its outer shape and inner microstructure to a changing mechanical environment. This is called functional adaptation by remodelling (Martin *et al*. 1998). Known as Wolff's Law (Wolff 1869, 1892, 1986), the mechanism of bone functional adaptation has long been of interest in the field of orthopaedics biomechanics, and mathematical models for bone remodelling based on experimental studies have been proposed (Huiskes *etal*. 1987, 2000; Cowin 1993; Jang & Kim 2008; Coelho *et al*. 2009; Dunlop *et al*. 2009; Gerhard *et al*. 2009; Martinez-Reina *et al*. 2009). The authors proposed a model for trabecular surface remodelling, and conducted computational simulations to elucidate the mechanism of bone functional adaptation (Adachi *et al*. 1997, 2001; Tsubota *et al*. 2002, 2009; Tsubota & Adachi 2004, 2005, 2006). Those studies proposed a uniform stress hypothesis at the remodelling equilibrium (Adachi *et al*. 1998) as a phenomenological hypothesis, and investigated the relationship between trabecular morphology and local mechanical environment. This demonstrated that local stress non-uniformity could be a candidate for the driving force of trabecular surface remodelling, and the proposed remodelling rate equation could predict the trabecular patterns corresponding to the macroscopic stress trajectory, as stated in Wolff's Law (Tsubota *et al*. 2009). However, possible mechanisms at the cellular level were not considered, even though it is known that remodelling is accomplished by coupling cellular activities of osteoclasts, osteoblasts and osteocytes (Parfitt 1994).

One important system that could relate the microscopic cellular level activities to the macroscopic tissue level change in trabecular morphology by remodelling is the cellular network (Cowin *et al*. 1991) in bone matrix and on bone surfaces. In bone matrix, many osteocytes are embedded and interconnected in an intercellular network system (Donahue *et al*. 1995; Kamioka *et al*. 2001, 2009; Knothe Tate *et al*. 2004; Sugawara *et al*. 2005). The osteocytes are regarded as mechanosensory cells (Mullender & Huiskes 1997; Burger & Klein-Nulend 1999; Cowin 2007; Huo *et al*. 2008; Adachi *et al*. 2009*a*,*b*), and their network system acts as the pathway of mechanical signals from the osteocytes to the cells on the trabecular surfaces (Adachi *et al*. 2009*c*). When bone is subjected to mechanical loading, bone matrix deformation produces an interstitial fluid flow in the lacuno-canalicular system (Weinbaum *et al*. 1994; Knothe Tate *et al*. 1998; You *et al*. 2001; Wang *et al*. 2007). This fluid flow is believed to play important roles in providing nutrients and removing wastes and also in cellular mechanotransduction (Burger & Klein-Nulend 1999; Klein-Nulend *et al*. 2005; Bonewald & Johnson 2008; Fritton & Weinbaum 2009). Considering that osteocytes are sensitive to fluid flow (Klein-Nulend *et al*. 1995*a*,*b*), it is hypothesized that osteocytes respond to the fluid-induced shear stress (Weinbaum *et al*. 1994) and deliver its mechanical information to the surface cells by intercellular communication (Adachi *et al*. 2009*c*); as a result, bone formation and resorption are regulated.

The purpose of this study is to propose a theoretical framework of trabecular bone remodelling that interconnects the microscopic cellular activities to the macroscopic morphological changes through the mechanical hierarchy by constructing a mathematical model of this complex adaptation mechanism. First, the mechanical behaviour of a trabecula with a lacuno-canalicular system is modelled as a poroelastic material (Cowin 1999; Kameo *et al*. 2008, 2009) to evaluate the interstitial fluid flow in bone matrix and the fluid shear stress on the osteocyte process membrane under dynamic mechanical loading. Second, a mathematical model of trabecular surface remodelling that considers mechanosensing by osteocytes (Mullender *et al*. 1994; Mullender & Huiskes 1995; Adachi *et al*. 2009*a*,*b*) and intercellular communication between osteocytes (Adachi *et al*. 2009*c*) is proposed. Then, a remodelling simulation is conducted by applying the finite-element method for a single trabecula under cyclic uniaxial loading to investigate the ability of functional adaptation by remodelling.

## 2. Quantification of the interstitial fluid flow

### (a) Equilibrium equation

When bone is subjected to mechanical loading, the bone matrix deforms and produces interstitial fluid flow in the lacuno-canalicular system. The mechanical behaviour of bone, which contains both solid and fluid phases, can be modelled as a poroelastic material (Cowin 1999). Bone matrix generally possesses anisotropic mechanical properties owing to microscopic structures, such as lamellae and canalicular systems. For simplicity, it is assumed here that bone matrix behaves as an isotropic elastic material. That is, the stress and strain tensors, *σ*_{ij} and *ε*_{ij}, are related through the following constitutive relation (Detournay & Cheng 1993; Wang 2000; Coussy 2004):
2.1
with the Einstein summation convention for the repeated subscript *k*, where *δ*_{ij} denotes Kronecker's delta. The parameters in equation (2.1) are the Biot–Willis coefficient *α*, the drained shear modulus *G* and the drained Poisson's ratio *ν*. The fluid (pore) pressure is denoted by *p*.

Substituting equation (2.1) into the equilibrium equation and using the displacement–strain relation, the equilibrium equation is derived as
2.2
in terms of the displacement *u*_{i} and pressure *p*, where *F*_{i} is the body force and ignored in this study.

### (b) Diffusion equation

Darcy's Law, which relates the pressure gradient *p*_{,j} to the fluid flux *q*_{i} in the direction of *x*_{i}, is generally expressed as
2.3
using the permeability tensor *k*_{ij} and the fluid viscosity *μ*. Substituting equation (2.3) into the law of conservation of mass leads to
2.4
where *ζ* is the variation of fluid content, defined as the variation of fluid volume per unit volume of porous material, and *Q* is the fluid source density.

Substituting the relations among the volumetric strain *ε* (=*ε*_{kk}), fluid pressure *p* and the variation of fluid content *ζ*:
2.5
into equation (2.4), in the absence of fluid source and using the displacement–strain relation, the diffusion equation is derived as
2.6
in terms of the displacement *u*_{i} and the pressure *p*, where *S*_{ε} is the constrained specific storage. Assuming the isotropic interstitial fluid flow, we set *k*_{ij} = *kδ*_{ij} here. By applying equations (2.2) and (2.6) to the finite-element method, the change in interstitial fluid pressure owing to changes in mechanical loading can be numerically analysed.

## 3. Mathematical model of trabecular bone remodelling

### (a) A modelling framework

A proposed modelling framework of the trabecular bone remodelling that considers the osteocytes' network system is as described below. Mechanical loading applied to bone induces bone matrix deformation and an interstitial fluid pressure gradient. This pressure gradient drives the interstitial fluid flow that produces fluid shear stress on the surface of the osteocyte processes. The osteocytes respond to the shear stress, and deliver its mechanical signal to the neighbouring osteocytes and the cells on the trabecular surface, such as osteoblasts and osteoclasts, through the intercellular network system. This framework is illustrated in figure 1.

A mathematical model of the trabecular bone remodelling was constructed based on this framework. First, the interstitial fluid pressure *p* and its spatial gradient ∇*p* under mechanical loading can be analysed using poroelasticity. Second, the interstitial fluid velocity *v*_{z}(*R*) in canaliculi and the fluid-induced shear stress *τ*_{p}(** n**) acting on the osteocyte processes can be estimated. Here, the normal vector

**indicates the direction of the cell process at point**

*n***, and it is assumed that the wall of the canaliculus and the surface of the process membrane were modelled as an annular pipe with a cylindrical coordinate system (**

*x**R*,

*Θ*,

*Z*), as illustrated in figure 1

*c*. Third, the mechanical stimulus

*S*

_{oc}(

**) that the osteocyte at point**

*x***can sense is estimated by integrating the shear stress**

*x**τ*

_{p}(

**) on the processes, and is communicated to the cell at point**

*n*

*x*_{sf}summarized as

*S*

_{sf}(

*x*_{sf}). Finally, the rate of trabecular surface remodelling , which describes the formation and resorption of the trabecular surface, is regulated by

*S*

_{sf}(

*x*_{sf}). In the following subsections, each process of this mathematical model is explained in detail.

### (b) Fluid shear stress owing to the interstitial fluid flow

The pressure gradient owing to bone matrix deformation induces an interstitial fluid flow in the lacuno-canalicular system. Referring to the microstructure model of Weinbaum *et al*. (1994), shown in figure 1*c*, it is assumed that interstitial fluid flow through a fibre matrix, such as proteoglycan, in an annular canaliculus can be governed by the Brinkman equation. Solving the Brinkman equation under no-slip boundary conditions on the walls, the fluid velocity *v*_{z}(*R*) in the axial direction is expressed as a function of the pressure gradient ∂*p*/∂*Z*, the fluid viscosity *μ*, the radius of the osteocytic process *r*_{p}, the radius of the canaliculus *r*_{c} and the permeability of the fibre matrix in the fluid annulus *k*_{p}. The permeability of the fibre matrix *k*_{p} is approximated by 0.0572(*Δ*/*a*_{0})^{2.377} using the fibre radius *a*_{0} and the spacing of the fibre matrix *Δ* (Tsay & Weinbaum 1991). Thus, the shear stress *τ*_{p}(** n**) acting on the osteocyte processes that align in direction

**, as indicated by an arrow in figure 1**

*n**b*, can be derived as 3.1 where constants

*A*

_{1}and

*B*

_{1}are given by 3.2

In these equations, *q* is the ratio of the radius of the canaliculus *r*_{c} to that of the process *r*_{p} (*q* = *r*_{c}/*r*_{p}), *γ* is the dimensionless parameter defined by and *I*_{0}, *I*_{1}, *K*_{0} and *K*_{1} are modified Bessel functions (Weinbaum *et al*. 1994).

### (c) Mechanosensing by osteocytes

It is known that the osteocytes respond sensitively to flow shear stresses acting on the processes (Klein-Nulend *et al*. 1995*a*,*b*). Thus, the shear stress on the process was assumed as the mechanical signal input to the osteocytes. Considering the volume fraction of canaliculi oriented to ** n**,

*ρ*

_{c}(

**), at position**

*n***, the model can take into account the local canalicular anisotropy. That is, by integrating the fluid-induced shear stress on the cell processes, the mechanical stimulus**

*x**S*

_{oc}(

**) sensed by the osteocytes per unit of volume can be defined as 3.3 where**

*x**θ*is the angle between the

*x*

_{3}-axis in the Cartesian coordinate system (

*x*

_{1},

*x*

_{2},

*x*

_{3}) and vector

**, and**

*n**φ*is the angle between the

*x*

_{1}-axis and the projection of

**onto the**

*n**x*

_{1}

*x*

_{2}-plane, measured anticlockwise. The mechanical stimulus is considered as its average over a sufficiently longer time scale than the repetitive daily loading history. For simplicity, the cell processes are assumed to extend isotropically in the radial direction, i.e. the canalicular volume fraction

*ρ*

_{c}does not depend on the orientation

**.**

*n*### (d) Intercellular communication of mechanical signals

The intercellular communication of mechanical signals through the network system, as illustrated in figure 2, is considered in this model. Assuming that cellular communication depends on the distance between cells and that there exists a maximum distance for intercellular communication, the mechanical signals are delivered to the bone-forming and -resorbing cells on the trabecular surface within distance *l*_{L}. Thus, using a monotonically decreasing function *w*(*l*) = 1 − *l*/*l*_{L} (*l* ≤ *l*_{L}) as a weight function, the mechanical signal *S*_{sf}(*x*_{sf}) that the cell at position *x*_{sf} on the trabecular surface can sense is defined as
3.4

This is a positive scalar function, indicating the cellular activity on the trabecular surface. The function integrates signals in the region *Ω* that can also cover the trabecular surface when the cells on the trabecular surface, such as bone lining cells and osteoblasts, are involved in the mechanosensory system. Therefore, we can consider the difference in mechanosensitivity in the different types of cells by changing weight functions depending on the location of the cells.

### (e) Rate equation of trabecular surface remodelling

Self-regulatory mechanism of bone remodelling is qualitatively considered. Bone formation and resorption are regulated by the increase and decrease in the mechanical signal from the physiological equilibrium condition. Therefore, by considering *S*_{sf}(*x*_{sf}) as the remodelling driving force, the rate of trabecular surface remodelling can be modelled as
3.5
where the mechanical signal at the remodelling equilibrium , the lazy zone and the maximum formation/resorption rate and are considered in the rate equation as shown in figure 3.

## 4. Simulation of trabecular bone remodelling

### (a) Single trabecular model

A two-dimensional single trabecula was modelled as shown in figure 4. A 200 µm wide trabecula, a skew angle (30°) and two parallel plates, 100 µm wide, were subjected to cyclic uniaxial loading. A plane strain condition was applied to this model. The region analysed was 1.0 × 1.0 mm and discretized by quadrilateral finite elements with an edge size of 10 µm. The permeability coefficient *k* was set as 1.1 × 10^{−21} m^{2} (Beno *et al*. 2006), and the other material parameters were: fluid viscosity *µ* = 1.0 × 10^{−3} Pa s, the drained shear modulus *G* = 5.94 GPa, the drained Poisson's ratio *ν* = 0.325, the bulk modulus for solid *K*_{s} = 17.66 GPa and for fluid *K*_{f} = 2.3 GPa and a porosity *φ* = 0.05 (Smit *et al*. 2002).

For boundary conditions, a shear-free condition was applied to the lower edge, and the inflow and outflow of interstitial fluid from all trabecular surfaces were assumed. On the upper trabecula, a cyclic loading *F*(*t*) = −*a *sin 2*πft* (*a* = 0.4 MPa and *f* = 1 Hz) was applied in the direction of the *x*_{2}-axis. The time-averaged shear stress in equation (3.3) was considered as the average value during one cycle.

The settings of the physiological parameters explained in the previous section were specified as follows. The radius of the osteocytic process *r*_{p} and the radius of the canaliculus *r*_{c} were determined based on the average values measured by You *et al*. (2004): *r*_{p} = 52 nm and *r*_{c} = 129.5 nm. The fibre radius *a*_{0} and the fibre spacing Δ were taken as *a*_{0} = 0.6 nm and *Δ* = 7 nm (Weinbaum *et al*. 1994). The volume fraction of canaliculi was set as *ρ*_{c}(** n**) = 1/

*π*, assuming their structural isotropy. The maximum remodelling rate = 40 µm d

^{−1}(Jaworski & Lok 1972) was based on the resorption rate of osteoclasts. The maximum distance for intercellular communication,

*l*

_{L}, and the parameters associated with the mechanical signal, , , and , were arbitrarily set as:

*l*

_{L}= 200 µm, =22.5 nN µm

^{−1}, =7.5 nN µm

^{−1}, =15 nN µm

^{−1}, and the parameter was set as 0 nN µm

^{−1}, for simplicity, i.e. the lazy zone was not considered. To express trabecular surface movement in the simulation, the level set method (Osher & Sethian 1988), which is a numerical technique for tracking interfaces and shapes of materials, was employed.

### (b) Results

Morphological change in a single trabecula by remodelling and the distribution of the interstitial fluid-induced shear stress acting on the osteocytic processes are illustrated in figure 5. In the initial state, as shown in figure 5*a*, a higher shear stress was observed at the trabecular surfaces around both ends with an acute angle. This generated large mechanical stimuli to the osteocytes and rapid bone formation was promoted on the surfaces, while at the middle of the trabecula, the original morphology was preserved. At 5 days, the surfaces with an obtuse angle had a lower shear stress, and bone resorption dominated in the region. After bone remodelling for 10 days, the concentration of fluid-induced shear stress was relieved, as shown in figure 5*c*. As a result of successive remodelling, the single trabecula aligned along the loading direction and shear stress in the neighbourhood of the trabecular surfaces became uniform at 50 days.

To quantify the morphological change in the trabecula, we focused on the bone volume fraction and apparent stiffness. Changes in the bone volume fraction and the apparent stiffness during remodelling are plotted in figure 6. The bone volume fraction increased slightly at early stages, and then gradually decreased and converged at a constant value owing to remodelling. In contrast, the apparent stiffness of the trabecula for the applied loading increased monotonically and reached a uniform plateau value.

## 5. Discussion

We proposed a mathematical model of trabecular surface remodelling that considers the osteocyte mechanosensory network system. Based on the model, we simulated the morphological changes of a single trabecula under cyclic uniaxial loading. As a result of bone remodelling, the trabecula aligned along the loading direction and the stress distribution on the trabecular surface became uniform. The apparent stiffness of the trabecula increased and the bone volume fraction decreased with the progress of bone remodelling.

One of the first mathematical models relating osteocyte mechanosensing to trabecular remodelling was presented by Mullender *et al*. (1994). They and their co-workers simulated the development of trabecular architecture by assuming the strain-energy density as a signal sensed by the osteocytes (Mullender & Huiskes 1995; Huiskes *et al*. 2000; Ruimerman *et al.* 2005). A bone remodelling algorithm that combined strain-adaptive remodelling and microdamage-induced resorption was presented under the hypothesis that both mechanical strain and microdamage accumulation drive cellular responses (McNamara & Prendergast 2007; Mulvihill & Prendergast 2008). The bone remodelling simulation based on these models provides reasonable results from a physiological point of view. Previously, we also developed a mathematical model for trabecular surface remodelling, which assumed that non-uniformity of local stress on trabecular surface drives bone remodelling (Adachi *et al*. 1997, 2001; Tsubota *et al*. 2002; Tsubota & Adachi 2004, 2005, 2006). Applying this algorithm to the trabecular remodelling simulation for an entire human proximal femur using large-scale voxel finite-element models, we demonstrated the creation of complex three-dimensional trabecular structure similar to that in the actual femur (Tsubota *et al*. 2009).

Although these mathematical models are based on the phenomenological hypothesis that the cellular responses are affected by the mechanical quantities at the tissue-scale level, the cellular mechanosensing will depend on the mechanical state at the microscopic level, such as lacuno-canalicular-scale level. In this study, we quantified the mechanical behaviour of interstitial fluid in lacuno-canalicular porosity. Assuming the fluid-induced shear stresses acting on the osteocyte processes as a driving force for bone remodelling, we presented the theoretical framework for remodelling across the mechanical hierarchy from the macroscopic changes in trabecular morphology to the microscopic cellular activities.

The trabecular orientation obtained in this simulation coincides with the observation of cancellous bone remodelling *in vivo* (Goldstein *et al*. 1991; Guldberg *et al*. 1997). This suggests the validity of our theoretical framework for trabecular bone remodelling. In general, the loss of volume causes a reduction in the stiffness in structural materials. However, the apparent stiffness of the trabecula increased despite the decrease in the bone volume fraction owing to bone remodelling. From a mechanical point of view, the trabecular morphology at 50 days is appropriate for its function of bearing an external load. These results indicate that the trabecula actively adapted to the mechanical environment by changing its own morphology.

Thus, the proposed model, which considers the osteocyte network system, can predict bone functional adaptation by remodelling. As observed in actual trabeculae, the lacuno-canalicular network system has a strong anisotropy that is perpendicular to the trabecular surfaces (Kameo *et al*. 2010). This anisotropy can affect the cellular mechanosensing and the intercellular communication of mechanical signals, as well as the interstitial fluid flow. To understand functional adaptation by bone remodelling, it is important to consider the contribution of this anisotropy to the mathematical model of bone remodelling. This development of a computational simulation for trabecular bone remodelling will enable us to gain new insights into its mechanism.

## Acknowledgements

This study was partially supported by a Grant-in-Aid for Specially Promoted Research (20001007) from the Ministry of Education, Culture, Sports, Science and Technology of Japan, and by a Grant-in-Aid from the Japan Society for the Promotion of Science Fellows (201475).

## Footnotes

One contribution of 13 to a Theme Issue ‘The virtual physiological human: computer simulation for integrative biomedicine I’.

- © 2010 The Royal Society