**2. Bone remodeling simulation**

To develop the computer simulation model of bone remodeling considering bone metabolism macroscopically, this section describes the bone remodeling model referring to the strainbased mechanostat theory submitted by Harold Frost and refined more than 5 decades. Physiological, low, and high strain ranges were defined in mechanostat theory, and different mechanisms were considered to describe the mechanical surface remodeling of trabecular architecture resulting in the structural change (Kwon et al., 2010a).

The following was assumed as the result of a remodeling turnover for each strain range: 1) in the physiological strain range, the occurrence of bone resorption or formation depends on the degree of the local stress non-uniformity; 2) in the low strain range, only bone resorption occurs owing to osteocyte apoptosis and its frequency increases with decreasing strain; 3) in the high strain range, only bone formation occurs by targeted remodeling and its frequency increases with strain. Each window is distinguished by threshold values of equivalent strain *ε<sup>c</sup>* in the context of the mechanostat theory, and low, physiological, and high strain range correspond to the strain ranges of *εc* smaller than *εdu*, between *εpl* - *εpu* and larger than *εol* (Fig 1).

Trabecular architecture is discretized using voxel finite elements, as schematically shown in Fig. 2. The trabecular surface movement by remodeling is expressed by adding and removing voxel elements on the surface according to the mechanical stress/strain conditions determined by a finite element analysis. Each trabecular element was classified by equivalent strain into low, physiological, and high strain range. Bone surface remodeling was assumed to occur according to the mathematical models of each window described in the following sections in detail.

**Figure 1.** Strain ranges categorized by the equivalent strain and in-between transition ranges.

**Figure 2.** Discretized trabecular bone.

#### **2.1. Physiological strain range: Mechanical adaptation by remodeling**

In the physiological strain range (*εpl* < *εc* < *εpu*), the non-uniformity of the stress/strain distribution was taken as driving stimuli (Adachi et al., 1997). For the voxel representing a point *c* on bone surface, the magnitude of local stress non-uniformity *c* was quantified as

$$\Gamma\_c = \ln \left( \sigma\_c \sum\_{i}^{N} w(l\_i^c) \right) \left( \sum\_{i}^{N} w\left(l\_i^c \right) \sigma\_i \right) \tag{1}$$

Osteocyte Apoptosis-Induced Bone Resorption

(4)

(5)

in Mechanical Remodeling Simulation – Computational Model for Trabecular Bone Structure 29

*c l*

1 and 1. *LS <sup>c</sup> P f* (6)

max

min

(7)

*du*

 

*HS* 1 and 1. *<sup>c</sup> P f* (8)

*l c*

*u c*

*c u*

where Γ*u* and Γ*l* denote the lower and upper threshold values, respectively. More specifically, voxel elements are added around the surface point *c* when *f*(Γ*c*) = 1, and the voxel *c* is subtracted when *f*(Γ*c*)=-1. The case of *f*(Γ*c*)=0 represents local remodeling equilibrium and no change occurs at and around *c*. Moreover, in the cases of 0 < Γ*c* ≤ Γ*u* and Γ*l* ≤ Γ*c* < 0, bone formation and resorption occurs stochastically according to probability *f*(Γ*c*). In this physiological strain range, the activation frequency *AFf* or *AFr*, defined as a relative

In the low strain range (*εc* < *εdu*), only osteocyte apoptosis-accelerated bone resorption is assumed to occur (Noble et al. 2003; Gu et al., 2005; Li et al. 2005). The probability of bone

The degree of osteocyte apoptosis increases with the decrease of strain. Thus, the activation frequency for bone resorption *AFr* is also elevated with the decrease of strain and is

> (0 ) ln ( ) ( ), ln( ) *<sup>r</sup> <sup>c</sup> <sup>a</sup> <sup>r</sup>*

 

In the high strain range (*εol < εc*), only bone formation is assumed to occur. The probability of

Bone formation in high strain range undergoes the targeted remodeling mechanism (Burr, 2002; Da Costa Gómez et al., 2005), and thus the activation frequency for bone formation *AFf*

max min

*AF AF <sup>a</sup>* 

min

*r c c du du*

*c*

**2.3. High strain range: Targeted remodeling by over use** 

bone formation as the result of remodeling turnover is unity, i.e.,

is also elevated with the increase of strain and is described by

*AF*

where *AFr* reaches the maximum *AFrmax* at *εc= εmin*.

1 1 sin <sup>1</sup> <sup>0</sup>

1 1 sin 1 0

*c*

*c*

2 2

2 2

1

1

rate for bone formation or resorption, is set to unity, i.e., *AFf* = *AFr* = 1.

**2.2. Low strain range: Osteocyte apoptosis- induced bone resorption** 

*r c l*

*P*

*P*

resorption is unity, i.e.,

described by

*f c u*

where *σc* denotes the stress at the point *c* and *σ<sup>i</sup>* denotes the representative stress at a point *i* with the distance *<sup>c</sup> i l* from *c*. The function *w*( *<sup>c</sup> i l* ) is described as

$$w\left(l\_i^c\right) = \begin{cases} 1 - l\_i^c \int l\_L & \left(0 \le l\_i^c < l\_L\right) \\ 0 & \left(l\_L \le l\_i^c\right) \end{cases} \tag{2}$$

defining a positive weight for the neighboring point *i* within the limit of sensing distance *lL* from the point *c* under consideration. This represents the sensory integration by osteocyte network.

The probability *f*(Γc) of bone resorption or formation at the point *c* is described by

$$f\left(\Gamma\_c\right) = \begin{cases} P\_f\left(\Gamma\_c\right) & \left(0 < \Gamma\_c\right) \\ 0 & \left(\Gamma\_c = 0\right) \\ -P\_r\left(\Gamma\_c\right) & \left(\Gamma\_c < 0\right) \end{cases} \tag{3}$$

Osteocyte Apoptosis-Induced Bone Resorption

in Mechanical Remodeling Simulation – Computational Model for Trabecular Bone Structure 29

$$P\_f\left(\Gamma\_c\right) = \begin{cases} \frac{1}{2} \left\{ \sin \pi \left( \frac{\Gamma\_c}{\Gamma\_u} - \frac{1}{2} \right) + 1 \right\} & \left( 0 < \Gamma\_c \le \Gamma\_u \right) \\\\ 1 & \left( \Gamma\_u < \Gamma\_c \right) \end{cases} \tag{4}$$

$$P\_r\left(\Gamma\_c\right) = \begin{cases} \frac{1}{2} \left\{ \sin \pi \left( \frac{\Gamma\_c}{\Gamma\_l} - \frac{1}{2} \right) + 1 \right\} & \left( \Gamma\_l \le \Gamma\_c < 0 \right) \\ 1 & \left( \Gamma\_c < \Gamma\_l \right) \end{cases} \tag{5}$$

where Γ*u* and Γ*l* denote the lower and upper threshold values, respectively. More specifically, voxel elements are added around the surface point *c* when *f*(Γ*c*) = 1, and the voxel *c* is subtracted when *f*(Γ*c*)=-1. The case of *f*(Γ*c*)=0 represents local remodeling equilibrium and no change occurs at and around *c*. Moreover, in the cases of 0 < Γ*c* ≤ Γ*u* and Γ*l* ≤ Γ*c* < 0, bone formation and resorption occurs stochastically according to probability *f*(Γ*c*). In this physiological strain range, the activation frequency *AFf* or *AFr*, defined as a relative rate for bone formation or resorption, is set to unity, i.e., *AFf* = *AFr* = 1.

#### **2.2. Low strain range: Osteocyte apoptosis- induced bone resorption**

In the low strain range (*εc* < *εdu*), only osteocyte apoptosis-accelerated bone resorption is assumed to occur (Noble et al. 2003; Gu et al., 2005; Li et al. 2005). The probability of bone resorption is unity, i.e.,

$$P\_{LS} = 1 \quad \text{and} \quad f\left(\Gamma\_c\right) = -1. \tag{6}$$

The degree of osteocyte apoptosis increases with the decrease of strain. Thus, the activation frequency for bone resorption *AFr* is also elevated with the decrease of strain and is described by

$$AF\_r(\mathcal{E}\_c) = \begin{cases} AF\_{r\max} & \text{( $0 < \mathcal{E}\_c < \mathcal{E}\_{\min}$ )}\\ \left(\frac{\mathcal{E}\_{du}}{\mathcal{E}\_c}\right)^a & \text{( $\mathcal{E}\_{\min} < \mathcal{E}\_c < \mathcal{E}\_{du}$ )}\\ \end{cases} \\ \text{or} \\ \begin{aligned} \mathcal{E}\_{r\max} & \text{( $a <\_{du} \mathcal{E}\_{r\min}$ )} \\ \end{aligned} \\ \text{and} \\ \begin{aligned} AF\_{r\max} & \text{( $\mathcal{E}\_{d} < \mathcal{E}\_{r\min}$ )} \\ \end{aligned} \tag{7}$$

where *AFr* reaches the maximum *AFrmax* at *εc= εmin*.

28 Apoptosis and Medicine

**Figure 2.** Discretized trabecular bone.

*i*

with the distance *<sup>c</sup>*

**Figure 1.** Strain ranges categorized by the equivalent strain and in-between transition ranges.

**2.1. Physiological strain range: Mechanical adaptation by remodeling** 

*l* from *c*. The function *w*( *<sup>c</sup>*

*w l*

In the physiological strain range (*εpl* < *εc* < *εpu*), the non-uniformity of the stress/strain distribution was taken as driving stimuli (Adachi et al., 1997). For the voxel representing a point *c* on bone surface, the magnitude of local stress non-uniformity *c* was quantified as

> ln *N N c c c c i ii i i*

where *σc* denotes the stress at the point *c* and *σ<sup>i</sup>* denotes the representative stress at a point *i*

 

*c c iL i L c i c*

*ll l l*

1 0

defining a positive weight for the neighboring point *i* within the limit of sensing distance *lL* from the point *c* under consideration. This represents the sensory integration by osteocyte network.

*c c*

*P*

*i*

0

The probability *f*(Γc) of bone resorption or formation at the point *c* is described by

*P*

*f*

*wl wl*

 

*l* ) is described as

*l l*

0 0

*rc c*

*f c c*

0

0

*L i*

(1)

(2)

(3)

#### **2.3. High strain range: Targeted remodeling by over use**

In the high strain range (*εol < εc*), only bone formation is assumed to occur. The probability of bone formation as the result of remodeling turnover is unity, i.e.,

$$P\_{HS} = 1 \quad \text{and} \quad f\left(\Gamma\_c\right) = \ 1. \tag{8}$$

Bone formation in high strain range undergoes the targeted remodeling mechanism (Burr, 2002; Da Costa Gómez et al., 2005), and thus the activation frequency for bone formation *AFf* is also elevated with the increase of strain and is described by

$$AF\_f\left(\mathcal{E}\_c\right) = \left\{\frac{AF\_{f\text{max}} - 1}{2} \left[\sin\pi\left(\frac{\mathcal{E}\_c - \mathcal{E}\_{ol}}{\mathcal{E}\_{\text{max}} - \mathcal{E}\_{ol}} - \frac{1}{2}\right) + 1\right] \quad \left(\mathcal{E}\_{ol} \le \mathcal{E}\_c < \mathcal{E}\_{\text{max}}\right)\right\} \tag{9}$$
 
$$AF\_{f\text{max}} \tag{9}$$

Osteocyte Apoptosis-Induced Bone Resorption

in Mechanical Remodeling Simulation – Computational Model for Trabecular Bone Structure 31

6. For voxels with activation frequency *AFr* or *AFf* larger than one, the same surface

**3. Simulation of trabecular bone structures for human proximal femur** 

The three-dimensional profile of human proximal femur is reconstructed using computed tomography (female, 70 years, hip osteoarthritis), the resolution of the images was 0.7 mm so that the entire volume could be captured with slice image of 512512 cubic voxels. As the random initial trabecular structure in proximal femur, an artificial trabecular structure was given by the random-placement of hollow spheres with inside and outside diameters of 2100 and 4200 μm, respectively (Kwon et al., 2010a). One voxel was transferred to one finite element of eight node brick element. The isotropic elastic body with Young's modulus of 20.0 GPa and Poisson's ratio of 0.3 was assumed for trabecular bone material. The cortical outer surface was fixed throughout the remodeling simulation, and the intra-structure, from the head to mid-shaft, was remodeled in accordance with the same remodeling procedure. The fixed boundary conditions are imposed on the distal end. In the daily loading condition, we considered the one-legged stance, abduction, and adduction, which accounted for 50, 25, and 25% of this condition as illustrated in Fig. 3. The loading magnitudes and frequencies of these three stances were shown in Table 1 (Beaupré et al., 1990; Adachi et al., 1997). The extreme ranges of motion of abduction and adduction were assumed. The trabecular structure obtained under this daily loading condition was referred to as the normal structure and also used for the initial structure for the remodeling simulation under reduced

movements in procedure (4) are repeated *AFr* or *AFf* times. 7. Procedures (2)-(6) are repeated until equilibrium is achieved.

**3.1. Voxel FE model of human proximal femur** 

weight-bearing conditions (Kwon et al., 2010b).

**Figure 3.** Standard weight-bearing cases at each stance (daily-loading condition).

The simulation results from initial random structure to an equilibrium state are shown in Fig. 4. In the diaphyseal region, random inner structures disappeared whereas the thickness of cortical bone increased with simulation step to the similar level observed in CT image (Fig.

**3.2. Healthy trabecular structures in proximal femur** 

where *AFf* reaches the maximum *AFfmax* at *εc* = *εmax*. It is noted that the pathological overuse is not considered here and the upper limit *εmax* is enforced for *εc* in high strain range.
