Strain Amplification Analysis of an Osteocyte under Static and Cyclic Loading: A Finite Element Study

Osteocytes, the major type of bone cells which reside in their lacunar and canalicular system within the bone matrix, function as biomechanosensors and biomechanotransducers of the bone. Although biomechanical behaviour of the osteocyte-lacunar-canalicular system has been investigated in previous studies mostly using computational 2-dimensional (2D) geometric models, only a few studies have used the 3-dimensional (3D) finite element (FE) model. In the current study, a 3D FE model was used to predict the responses of strain distributions of osteocyte-lacunar-canalicular system analyzed under static and cyclic loads. The strain amplification factor was calculated for all simulations. Effects on the strain of the osteocyte system were investigated under 500, 1500, 2000, and 3000 microstrain loading magnitudes and 1, 5, 10, 40, and 100 Hz loading frequencies. The maximum strain was found to change with loading magnitude and frequency. It was observed that maximum strain under 3000-microstrain loading was higher than those under 500, 1500, and 2000 microstrains. When the loading strain reached the maximum magnitude, the strain amplification factor of 100 Hz was higher than those of the other frequencies. Data from this 3D FE model study suggests that the strain amplification factor of the osteocyte-lacunar-canalicular system increases with loading frequency and loading strain increasing.


Introduction
While bone-forming cells osteoblasts and resorbing cells osteoclasts make up only small proportions (<5 and <1%, resp.), osteocytes make up ∼90-95% of all bone cells in adult skeletons [1,2]. Osteocytes are the longest living bone cells, which can live up to decades within their mineralized environment and have the potential to live as long as the organism itself [3]. The osteocytes are believed to be the cells that sense the mechanical loads of the bone. The major known function of osteocytes is to translate the mechanical strain signals into biochemical signals between osteocytes and cells in the bone [4][5][6]. Osteocytes are embedded within the bone matrix, can remodel the perilacunar matrix (PCM), and are connected to each other by slender cell processes located within the small bony tubes of the canaliculi [1,4,7].
Osteocytes are surrounded by the PCM, which is a thin hyaluronan-rich coating. The PCM plays an important role in many cell-surface phenomena; it transfers mechanical signals between the cell body and cell processes and between PCM and the surrounding extracellular matrix (ECM) [8][9][10][11][12][13]. The PCM represents the entire space between the ECM and the osteocyte cell [14]. Due to these functional relations, the osteocyte-lacunar-canalicular system should be investigated based on the cell-PCM-ECM geometric structure ( Figure 1).
With the known biomechanosensory function of the osteocyte-lacunar-canalicular system in bone, dynamic mechanical stimulus like cyclic loading has been frequently used for the microscale environment biomechanics investigation [15]. Hsieh and Turner studied the anabolic effect of mechanical loading on bone tissue by modulating the loading frequency in the range of 1-10 Hz [16]. Gururaja et al. developed a 2D poroelastic model to investigate load-induced fluid flow in the lacunar-canalicular system of bone subjected to harmonic bending excitations in the range of 1-10 6 Hz [17]. Wang et al. investigated the net tracer transport in bone via the lacunar-canalicular porosity by using a mathematical model during cyclic mechanical loading at loading frequency 0.5-100 Hz [18]. To date, some researchers have investigated the osteocytelacunar-canalicular system by experimental and computational models [14,17,[19][20][21][22], and the unit of microstrain has been widely used in the loading and data analyses in the cell modelling. McCreadie and Hollister found the strain levels in and surrounding osteocytes were 1700 to 2700 microstrains [21]. Burr et al. measured the strain by surgically implanted strain gages and found that the strain was below 2000 microstrains even under conditions of strenuous activity [23]. Verbruggen et al. conducted a computational investigation of the strain amplification of osteocytes in vivo in response to mechanical stimulation under vigorous physiological activity (3000 loading levels) [14]. In addition, researchers have introduced methodologies to reconstruct biorealistic osteocyte models. A number of methods have been used to study 3D reconstruction of the osteocytelacunar-canalicular system. The methods include scanning electron microscopy (SEM) [24][25][26][27], serial focused ion beam/scanning electron microscopy (FIB/SEM) [28], transmission electron microscopy (TEM) [29,30], ultra high voltage electron microscopes (UHVEM) [31], light microscopy (LM) [32,33], confocal laser scanning microscopy (CLSM) [14,[34][35][36][37], atomic force microscopy (AFM) [38][39][40], computed tomographic (CT) scanning (such as nano-CT and synchrotron radiation micro-CT) [36,37,41,42], and X-ray (like X-ray phase nanotomography and X-ray microscope (TXM)) [43][44][45]. However, only a few of these models were used in finite element (FE) analysis of the osteocyte-lacunarcanalicular system. Verbruggen and his coworkers developed FE models of lacunar-canalicular network using confocal microscope imaging in rats [14,46].
While most previous computational geometries were 2D models, several researchers have proposed 3D models (including biorealistic and idealized FE models) for the computational investigation [14,21,22,[47][48][49]. For example, a 3D poroelastic FE model of rat tibia was developed to study the mechanical loads in modulating local flow distributions and concentration gradients within bone tissue. However, due to the simplicity of the geometry, only a thin layer of elements on the endosteal and periosteal surfaces was used in defining the boundary conditions [49]. As results of FE studies depend on many conditions (e.g., geometry, mesh density, element type and order, and material properties), more credible results can be obtained if more accurate geometries are used in the 3D models, when the other parameters have been confirmed.
Thus, although the osteocyte-lacunar-canalicular system has been investigated through the pure experimental tool, mathematics theory, or 2D FE models in previous studies, only a few researchers have used the 3D FE model. In the current study, a 3D FE model was proposed to predict the biomechanical behaviour of the osteocyte-lacunarcanalicular system. Firstly, a 3D geometry of an osteocytelacunar-canalicular system incorporating cell body, PCM, and ECM, which is similar to the real dimension, was developed. Secondly, the biomechanical behaviour of the system was investigated by using the 3D model subjected to the static and cyclic loading of different frequencies and magnitudes.

Geometric Modelling and
Meshing. An idealized model of osteocyte-lacunar-canalicular system was developed ( Figure 2). The triaxial lacunar osteocyte ellipsoid with the Cartesian equation for a general ellipsoid is shown as [50] ( ) where , , and are the osteocyte lacunar semiaxes in the , , and direction, respectively. is the minor axis, the major axis, and the intermediate axis of the osteocyte lacuna in the local coordinate system. In this study, , , and are 4.6, 9.45, and 2.4, respectively, with the values obtained from human femurs [42]. In this system, the PCM and osteocyte are embedded in the ECM block, whilst the osteocyte is in the center of the ECM block. While the mean length of each direction of ECM is 43 m for human [50], osteocyte cell bodies are surrounded by a pronounced ∼0.5-1 m thick layer of PCM [51]. The PCM thickness was assumed to be 0.75 m in this study. Beno et al. [50] estimated the numbers of canaliculi emanating from osteocyte lacunae for different species (six species: chick, rabbit, bovine, horse, dog, and human) by using the slicing method and surface area method based on the data from Remaggi et al. and Ferretti et al. [52,53]. Canaliculi were idealized as straight cylindrical channels and eighteen canaliculi were modelled as channels in the ECM [50] and using diameter of 0.25 m [4,54]. Furthermore, the average width of the canaliculi around the osteocyte process was ∼80 nm [4,55].
The FE analysis software ABAQUS 6.12 (SIMULIA, Providence, RI, USA) was used for simulations assuming fully saturated media. Because the model is symmetrical, only 1/8 symmetry model was used in all the simulations ( Figure 3). The reduced integration first order solid elements in Abaqus suffered from hourglassing when a mesh was coarse. In addition, the fully integration first order solid elements exhibited shear locking. To avoid the elements suffering from hourglassing and shear locking, the reduced integration with second order solid elements is recommended [56]. Twentynode hexahedral reduced integrated elements (C3D20R)  were used for all regions. The hexahedral FE mesh was mapped according to the high quality mesh geometries using twenty-node C3D20R elements in the whole FE analysis. Total number of nodes and elements were 453199 and 107868, respectively, in the 1/8 symmetry model (Table 1). Tie contact interfaces were used to ensure the PCM attached to the ECM and osteocyte and to prevent any relative movement during the simulation.

Materials.
All materials were assumed to be isotropic and linearly elastic. The PCM material properties are difficult to measure directly because of its size and connectivity to both cell body and ECM [10]. Meanwhile, it is hard to use the experimental data to define the material properties of the PCM surrounding the osteocyte [14]. An inverse FE approach was used to calculate the linear elastic properties of the PCM. The obtained PCM moduli ranged from 43 to 240 kPa based on this method [10]. Therefore, an elastic modulus of 43 kPa was assigned to the PCM and with Poisson's ratio of 0.4 [9,14,57]. The properties of the cells in the lacuna space were assigned a Young's modulus of 10 MPa and Poisson's ratio of 0.4 [21]. A modulus of 4471 ± 198 Pa and Poisson's ratio of 0.3 were attributed to the osteocyte cell body and processes [58].

Boundary and Loading Conditions.
To simulate the mechanical environment of the lacunocanalicular system, static and cyclic loads were applied and the responses of strain and stress distributions were analyzed. The ramped static compressive loads of 500, 1500, 2000, and 3000 microstrain were used in this study [14,22]. Generally, for most animals, the peak principal compressive strains are about 2000 to 3000 microstrains. However, peak principal strain in the human tibia at the location measured is slightly less than these values. In some cases such as uphill and downhill zigzag running, the principal compressive and shear strains are the greatest, reaching nearly 2000 microstrains, about three times higher than that recorded during walking [23].
Furthermore, a sinusoidal loading = | 0 [sin( )]| was used for compression [59] (Figure 4). The maximal loading magnitudes 0 were 500, 1500, 2000, and 3000 microstrains, respectively, in the simulations. The axial displacement 6 BioMed Research International boundary conditions were applied to the model faces to simulate the compression. The nodes on the opposing faces of applied displacement loading were constrained symmetrically to prevent rigid body motion [14,22,23]. To examine the influence of frequencies on strain, the loading frequencies of 1, 5, 10, 40, and 100 were investigated [16,17,59,60].

Static Loading.
To simulate the mechanical environment of the lacunocanalicular system, the global static compressive loads of 500, 1500, 2000, and 3000 microstrains were applied and the responses of strain and stress distributions were analyzed. The strain amplification factor was computed as the local maximum strain divided by the applied global strain as described [9,22]. The contour plots of half symmetry FE model (Figures 5(a)-5(d)) show the predicted strain distribution of the lacunar-canalicular FE model under 500, 1500, 2000, and 3000 microstrains, respectively, and the maximum principal strains were about 1705, 5129, 6848, and 10300 microstrains, respectively. When subjected to a nominal continuum strain level approximately equal to that measured in humans in vivo during rigorous activity like 2,000 microstrains, the osteocyte level strains can be as high as 12,000 to 15,000 microstrains (1.2% to 1.5%) [61]. Obviously, the maximum principal strain is located in the canaliculi near to osteocyte cell body and the interface between canaliculi and osteocyte. For the osteocyte cell, the principal strains in different areas followed the following trend: innermost center > center > outermost section. For PCM, the principal strains in the top and bottom fractions near to osteocyte cell body were higher than the rest, whilst those in the left and right side sections were lower than the other parts in the ECM.
The predicted results of strain amplification factors of the FE model were calculated for different loads (Figure 6), which were compared with the data in the literature [14,22]. The strain amplification factors were 3.41, 3.419, 3.24, and 3.43 for 500, 1500, 2000, and 3000 microstrain global loads, respectively. The strain amplification factors for loads of 500 and 3000 microstrains were comparable with the values of the idealized model and confocal image-derived model from Verbruggen et al. [14], while that for the loading of 1500 microstrains was comparable with the values of the idealized model from Verbruggen et al. [14] and Nicolella et al. [62]. Finally, the strain amplification factor for 2000-microstrain loading was comparable with those from the idealized model of Rath Bonivtch et al. [22] and Nicolella et al. [62].  The influence of loading frequency and loading magnitude on strain amplification was investigated, and the trends of strain amplification factors at different frequencies (1,5,10,40, and 100 Hz) and loading strains (500, 1500, 2000, and 3000 microstrains) were shown in Figure 12 for / = 0.5.

Discussion
In this work, a 3D FE model of osteocyte-lacunar-canalicular was developed to investigate the effects of cyclic loading at different frequencies on strain responses of the osteocytelacunar-canalicular system. The compressive loads of 500, 1500, 2000, and 3000 microstrains were applied in both static simulation and cyclic simulation. The strain amplification factors of all simulations were also calculated.
Some previous computational studies have examined idealized models of the osteocyte mechanical environment in 3D [21,22] and have studied both idealized models and biorealistic models [14,46]; these previous studies represented the first development of accurate 3D FE geometries of the lacuna-osteocyte-canaliculus-ECM system to predict osteocyte mechanobiology. In the previous studies, the shapes of idealized lacunar osteocyte were modelled as revolved  ellipsoid with the major and minor axes for idealized models. However, in our study, one triaxial lacunar osteocyte ellipsoid FE model was developed to predict the stain of osteocytelacunar-canalicular system for human.
In the static simulation, the strain amplification factors were compared with the existing data in literature ( Figure 6). The results of the simulation were found to be in good agreement with the results in literature [14,22,62]. The  [22]. Compared with the average strain amplification factors (varying from 1.1 to 3.8) from another study [62], the results of our FE analysis were consistent with these measurements. Thus, the simulation of static loading with our system was validated by the reference data in the literature, suggesting that our FE model of osteocyte-lacunar-canalicular is reliable.   For the cyclic loading simulations, the loading frequency and loading magnitude were investigated in this study, and we performed a total of 20 simulations with different loading strain ad 5 different frequencies (1,5,10,40,and 100 Hz) which are within the physiologically meaningful frequency range of 1-100 Hz [17,59,63]. The maximum strain increased with increasing loading frequency and increasing loading magnitude (loading strain). The maximum principal strain occurred at / = 0.5, which matched the maximum value of loading strain. Because the sinusoidal loading was applied in the cyclic loading simulations, the principal strain varied with time. When the loading strain reached the maximum magnitude value, the strain amplification factor of 100 Hz was higher than the other frequencies. When = 0.5 , the strain amplification factor increased with loading frequency and with loading strain. However, the strain amplification factor increased only slightly with the increasing loading magnitude for the same loading frequency and increased   slightly with the increasing loading frequency for the same loading magnitude.
There has been a lack of experimental data regarding the strain amplification factors for the osteocyte-lacunarcanalicular system. Due to the complex microstructural organization of osteocyte-lacunar-canalicular system, the osteocytes cannot be reliably and easily estimated from global strain measurements [36]. The results of static loading in our study are comparable with FE analysis estimates [14,22] and digital image measurements [36]. We have to point out that the cyclic loading data of our model is yet to be validated in future studies, as it is beyond the scope of present work.   Furthermore, since pore fluid pressure is an important loadinduced phenomenon of the osteocyte-lacunar-canalicular system as it in large part dictates the shear stress and degree of chemotransport experienced by the cells, future studies will be required to investigate the pore pressure, flow velocity, and poroelasticity in the 3D poroelastic FE model. Such data will facilitate the understanding of the biomechanical responses to loading and chemotransport around the osteocyte cell body in the osteocyte-lacunar-canalicular system.
Taken together, to understand the changes in strain of osteocytes with different loading frequencies and loading magnitudes, the current study has developed a 3D FE model of osteocyte-lacunar-canalicular system which was used to predict the responses of strain distributions of the osteocyte system under various static and cyclic loads. The study found that the strain amplification factor increased with loading frequency and increasing loading strain.