Osteocyte lacunar strain determination using multiscale finite element analysis.

Osteocytes are thought to be the primary mechanosensory cells within bone, regulating both osteoclasts and osteoblasts to control load induced changes in bone resorption and formation. Osteocytes initiate intracellular responses including activating the Wnt/β-catenin signaling pathway after experiencing mechanical forces. In response to changing mechanical loads (strain) the osteocytes signal to cells on the bone surface. However, this process of osteocyte activation appears heterogeneous since it occurs in sub-populations of osteocytes, even within regions predicted to be experiencing similar global strain magnitudes determined based on traditional finite element modeling approaches. Several studies have investigated the strain responses of osteocyte lacunae using finite element (FE) models, but many were limited by the use of idealized geometries (e.g., ellipsoids) and analysis of a single osteocyte. Finite element models by other groups included more details, such as canaliculi, but all were done on models consisting of a single osteocyte. We hypothesized that variation in size and orientation of the osteocyte lacunae within bone would give rise to micro heterogeneity in the strain fields that could better explain the observed patterns of osteocyte activation following load. The osteocytes in our microscale and nanoscale models have an idealized oval shape and some are based on confocal scans. However, all the FE models in this preliminary study consist of multiple osteocytes. The number of osteocytes in the 3D confocal scan models ranged from five to seventeen. In this study, a multi-scale computational approach was used to first create an osteocyte FE model at the microscale level to examine both the theoretical lacunar and perilacunar strain responses based on two parameters: 1) lacunar orientation and 2) lacunar size. A parametric analysis was performed by steadily increasing the perilacunar modulus (5, 10, 15, and 20 GPa). Secondly, a nanoscale FE model was built using known osteocyte dimensions to determine the predicted strains in the perilacunar matrix, fluid space, and cell body regions. Finally, 3-D lacunar models were created using confocal image stacks from mouse femurs to determine the theoretical strain in the lacunae represented by realistic geometries. Overall, lacunar strains decreased by 14% in the cell body, 15% in the fluid space region and 25% in the perilacunar space as the perilacunar modulus increased, indicating a stress shielding effect. Lacunar strains were lower for the osteocytes aligned along the loading axis compared to those aligned perpendicular to axis. Increases in lacuna size also led to increased lacunar strains. These finite element model findings suggest that orientation and lacunar size may contribute to the heterogeneous initial pattern of osteocyte strain response observed in bone following in vivo applied mechanical loads. A better understanding of how mechanical stimuli directly affect the lacunae and perilacunar tissue strains may ultimately lead to a better understanding of the process of osteocyte activation in response to mechanical loading.

Introduction

Osteoporosis and osteopenia are diseases characterized by low bone mass while osteopetrosis and osteosclerosis are characterized by increased bone mass. These diseases result in bone mineral density and architectural changes that affect the way bone responds to mechanical stimuli (van Hove et al., 2009). In the U.S. alone, osteoporosis leads to 1.5 million fractures per year (Black and Rosen, 2016). Bone is a dynamic organ that adapts its material and structural properties in response to its loading environment (Carter, 1984; Turner et al., 1992). According to Wolff's theory, bones adapt their structure and mass in response to mechanical stimuli to optimize their load bearing capacity (Santos et al., 2009). Mechanotransduction is a process in which the physical/mechanical forces applied to bone are transformed into biological responses. Bone tissue consists of three main cells, namely osteocytes, osteoclasts and osteoblasts, of which osteocytes comprise over 90% of the bone cells. Osteocytes are widely considered to be the primary mechanosensory cells in bone that are responsible for orchestrating bone formation and resorption (bone remodeling) through the regulation of osteoblasts and osteoclasts (Huang and Ogawa, 2010; Lanyon, 1993). Osteocytes occupy a fluid-filled space within their lacunae, and are interconnected by dendritic processes that pass through a network of small channels called canaliculi. In response to mechanical loading, osteocytes signal to osteoclasts and osteoblasts to control the function of these cells. The process of mechanotransduction, as is currently known, includes the response of the osteocyte to both direct mechanical deformation/strain as well as to fluid flow in the lacunocanalicular space. Mechanical strain is at least one driver of FFSS, which is commonly considered the more relevant stimulus for bone deposition. As a region is compressed or expands, fluid will move to equilibrate pressure throughout the network, resulting in fluid flow across cell membranes.

Previous studies (Forwood et al., 1996) have shown that dynamic loading has a stronger influence on bone remodeling compared to static loading. Mechanical loads are more effective at increasing bone formation through short durations of discrete loads followed by recovery periods rather than continuous application of the loading (Robling and Turner, 2002). Previous experiments on rat tibiae (Burr et al., 2002) suggest that a recovery period of 4–8 hours is sufficient for bone to reestablish its mechanically sensitive state. While dynamic loading is important to represent the full physiological behavior we have used quasi-static loading since we are considering only elastic behavior in this study.

Bone matrix is a heterogeneous material and structural components such as canaliculi, osteocyte lacunae and perilacunar matrix regions are considered to be potential stress concentrators as local microstructural strains at these locations are vastly different from the macro structural or global bone strains (Nicolella et al., 2006). In addition it has been hypothesized that osteopenic bone tends to have a larger fluid space between the osteocyte and the lacuna resulting in a lower stimulation due to mechanical loading, while the fluid space in an osteopetrotic bone is smaller, resulting in a greater stimulation (van Hove et al., 2009). Consequently, the overall motivation of this research is to study the local mechanical response (strain) on individual osteocytes due to mechanical loading of the bone and how this response is affected by differences in orientation of the cells. Future applications of the results of this study could be to study differences in the osteocyte's mechanical response due to the morphological changes that occur under normal and diseased conditions.

Finite element (FE) models are a good tool for estimating strain responses at the microstructural level. Therefore, researchers have used FE models to understand both the mechanical and biological mechanisms of increased bone formation (Cowin et al., 1991). One cellular response initiated in osteocytes after experiencing mechanical forces is activation of the Wnt/β-catenin signaling pathway (Bonewald and Johnson, 2008). However, even within regions in which traditional FE models have predicted similar global strains, the process of activation is not uniform among all osteocytes. Our previous studies have shown that the initial response to in vivo loading occurs within sub-populations of osteocytes (Lara-Castillo et al., 2015a). This suggests that there is a “strain threshold” for activation of osteocytes. In addition, osteocyte lacunae vary substantially in geometry, orientation and other factors within a single mouse bone.

Several studies have investigated the strain responses of osteocyte lacunae using FE models, but many were limited by the use of 2-D models or idealized geometries (e.g., ellipsoids) (Stern, 2013; McCreadie and Hollister, 1997). In this study, a multi-scale computational approach is used to study the strain response at the osteocyte cellular level by creating FE models at the micro and the nano scale level. At the microscale level, where the element sizes were on the order of a few micrometers, numerical lacunar and perilacunar strain responses were examined based on lacunar orientation and lacunar size. The biological motivation to study lacunar size, elongated shape, and orientation (the preferred alignment being parallel to the long axis of the bone) of the cells in the cortical region of the bone is based on observations in various normal or diseased states in both cancellous and calvarial bone (Sugawara et al., 2013).

The perilacunar region surrounding the lacuna is known to have varying material properties and is altered during lactation (Qing et al., 2012; Nicolella and Lankford, 2002). Hence, a parametric analysis was performed by steadily increasing the perilacunar modulus (5, 10, 15, and 20 GPa). Subsequently, the osteocyte cell, fluid space, lacunar and perilacunar region was resolved at a finer level in a nanoscale FE model containing a few lacunae, based on known osteocyte dimensions in order to determine the predicted strains in the perilacunar matrix, fluid space, and cell body regions. Finally, to study strains in an in situ situation, we created three dimensional lacuna models using confocal image stacks from mouse femurs to determine the strain in lacunae represented by realistic geometries derived from empirical data.

Finite element modeling methods

This section describes the methodologies used to develop the finite element 2-D microscale and nanoscale models and 3-D osteocyte lacunar models. The 2-D models presented in this paper are simplified representations of the 3D structure. The following sections describe the node and element details, along with the associated loading and boundary conditions. Table 1 outlines the material properties used for the models.

Table 1

Material properties for bone models.

Material	Elastic modulus	Poisson's ratio
Bone matrix	25 GPa	0.3
Cell body & fluid space	0.1 MPa	0.45
Perilacunar matrix	5–20 GPa	0.3

FEBio Software Suite (Maas et al., 2011; Maas et al., 2012) was used to create the microscale and nanoscale FE models of the osteocyte. The pre-processing and post-processing steps in FEBio were completed using Preview and PostView, respectively. Microscale (Micro Model) and nanoscale (Nano Model) versions of the model consisted of the bone matrix, osteocyte lacunae, and perilacunar region. The 3-D lacunar models were developed from confocal image stacks of mouse femurs using Mimics® Innovation Suite (Materialise P, MI, n.d.). A linear elastic FE analysis was used for all the models.

Microscale model development

Based on previous work (Bonivtch et al., 2007), bone matrix tissue was modeled as an isotropic elastic material with three material regions in the Micro model namely bone matrix, lacuna and perilacunar matrix. The perilacunar matrix region was defined as the region 10 μm around the edge of the lacuna (Bonivtch et al., 2007). The size of the bone matrix used for the Micro model was 1 mm × 2 mm × 0.01 mm. A hexahedral mesh was used with element size of 10 μm × 10 μm × 10 μm with the element size chosen to approximately represent one lacuna. This size was chosen based on the study of osteocyte cell/lacunae sizes in human female tibia in osteoarthritic, osteopenic and osteopetrotic conditions (van Hove et al., 2009). Fig. 1 shows the Micro model with several lacunae labeled as follows: L1, L2, L3 and L4 are four identical lacunae at different locations, AD and BD are two lacunae oriented diagonal to the longitudinal loading direction, V represents lacunae aligned in the vertical or along the loading direction and H are lacunae perpendicular to the loading direction. The Micro model consisted of 40,602 nodes and 20,000 elements. An elastic modulus value of 25 GPa was selected for the bone matrix and 0.1 MPa was selected for the lacunar region (Bonivtch et al., 2007).

Fig. 1

2-D FE Micro model (1 mm × 2 mm × 0.01 mm) showing lacunae of different idealized shapes, sizes and orientations. The loading direction for this model is vertical on the page, corresponding to the longitudinal axis of the bone. The fixed surface is indicated by red X markers. Green represents bone, and purple represents lacunae, with darker pixels representing the perilacunar regions for each shape. Top left - L1, L2, L3 and L4 are four lacunae with sides of equal size. Top right – lacunae with same size but different orientations: 2 V is vertical, 2H is horizontal, 1 CE is diagonal from upper left to lower right, and 1BD is diagonal from upper right to lower left. Bottom – magnified image showing greater detail of the perilacunar regions for each idealized cell.

Boundary conditions and calculations of strain

A sinusoidal load was applied to the top surface of the model and all the nodes on the bottom surface were fixed to prevent motion in the three directions. A sinusoidal is typically applied in all our in vivo animal loading experiments. However, the load value is low and the frequency value essentially results in a quasi-static loading. Since only elastic material behavior is considered, the loading type will not have an effect on the results. Previous work on the strain response in osteocyte lacunae (Bonivtch et al., 2007) proposed a global strain of 2000 microstrain which was considered to be the maximum physiological strain experienced by human bone under active conditions. For this study, the applied theoretical load was adjusted until a global strain of 2000 microstrain was achieved. The global strain was calculated by averaging the strains from nine elements in the region of bone at the center of the model. All theoretical strain values reported in this study corresponded to the longitudinal Lagrangian strain in the axial direction.

Nano model development

The nanoscale FE model (i.e., Nano model) was developed based on the actual size of the osteocyte, which is approximately 20 μm × 8 μm × 4 μm as shown in one study (Rochefort et al., 2010). Other studies (Mader et al., n.d.; Heveran et al., 2018) have developed in silico models of the osteocyte lacunar systems with varying spatial distribution and alignment of the cells. Each lacuna was modeled as an oval shape in order to closely resemble an actual lacuna. A solid, rectangular geometry with dimensions of 80 μm × 160 μm × 4 μm was created for the bone matrix, and the hexahedral meshing method was used to adjust the size of each element to 0.8 μm × 0.8 μm × 4 μm. In the Nano model, there were four material regions which included the fluid flow space and cell body regions (See Table 1 for the material properties of each region). The Nano model also consisted of 40,602 nodes and 20,000 elements.

As shown in Fig. 2A, six osteocytes were created with two osteocytes oriented perpendicular to the loading direction (H1 and H2), two aligned with the loading direction (V1 and V2), and two oriented at forty-five degrees to the loading direction (D1 and D2). The osteocyte cell to cell distance has been shown to be in the range of 35–40 μm, which was used as a guideline for this model (Sugawara et al., 2008). Theoretical strains were analyzed in three osteocyte sub regions, namely the cell body, fluid space, and perilacunar matrix region (Fig. 2B). The perilacunar region was chosen to be two elements wide representing a region about 1.6 μm wide, in this model based on an earlier work (Bonivtch et al., 2007). Strain ratios, computed as the ratio of the maximum strain in each region to the average bone strain determined at the center of the model, were also calculated.

Fig. 2

2-D FE Nano model with six osteocytes. The loading direction for this model is vertical from top to bottom and the fixed surface is indicated by red X markers . (A) The six osteocytes represent each possible orientation in the matrix. (B) A Representative osteocyte showing the bone matrix region (light green), perilacunar matrix (blue), fluid space (red), and cell body (purple).

Parametric analyses of Micro and Nano models

The parameters considered for analysis included the orientation and size of the osteocyte lacuna. While examining each parameter separately, we also investigated the effect of changing the perilacunar elastic modulus on the strain response. The tissue strain amplification was studied by adjusting the perilacunar modulus with values based on earlier work (Bonivtch et al., 2007), which suggested that the perilacunar tissue may alter its properties to modulate the strain environment around the embedded osteocyte.

Strain variation based on orientation in Micro model

This analysis considered the orientation of the osteocyte as another contributor to the differences in strain magnification because osteocytes in vivo can be aligned in a heterogeneous manner. Hence, in the Micro model we considered four lacunae that were oriented differently with respect to the loading surface. To minimize the effect of position, we only considered lacunae located at similar distances from the loading surface. In Fig. 1, the lacunae were named 1AD CE, 1BD, 2H and 2V with 1 representing diagonal orientations and 2 representing non-diagonal orientations. D stands for Diagonal while H and V stand for Horizontal and Vertical, respectively. The notations A and B were added to distinguish the two diagonally oriented lacunae.

Strain variation based on size in Micro model

A third parameter analyzed was the variation in osteocyte lacunar size. Three lacunae of different sizes were selected with all three possessing square-shaped geometries to minimize the effect of their respective orientations. The lacunar sizes were: S1 for 1 × 1, S2 for 2 × 2 and S3 for 3 × 3 element size lacuna. Since the size of each element was 10 μm × 10 μm, the size of lacuna S1 was 10 μm × 10 μm, lacuna S2 was 20 μm × 20 μm and lacuna S3 was 30 μm × 30 μm. The effect of position was minimized by selecting all three lacunae that were in close proximity to each other.

3-D lacuna model

The 3-D lacunar models were developed by importing thresholded confocal microscope Z stacks of osteocytes from sections of mouse femur into MIMICS from Materialise® (Materialise P, MI, n.d.). The confocal images were collected from thick 50 μm decalcified transverse bone sections from the femoral midshaft of 5 month old C57BL/6 female mice which were generated during a previous study, and the fixation, decalcification, staining and mounting of these bone specimens is described elsewhere (Tiede-Lewis et al., 2017). The osteocyte cell bodies were stained using the membrane stain 3,3′-dioctadecyloxacarbocyanine perchlorate (DiO), the nuclei were stained with 4,6-diamidino-2-phenylindole dihydrochloride (DAPI) and the lacunocanalicular systems had been labeled by injecting the mice intravenously with a Texas Red conjugated lysine fixable 10 kDa dextran (32 mg/kg in PBS), 4–5 min prior to euthanization, which distributes into the lacunae and canaliculi. These staining methods were described previously (Tiede-Lewis et al., 2017; Kamel-ElSayed et al., 2015) and the fluorescent dyes were obtained from ThermoFisher Scientific, Waltham, MA. For the current FE modeling study, the stained sections were imaged on a Leica TCS Sp5 II confocal microscope (Leica Microsystems, Wetzlar, Germany) in resonant scanner mode to collect detailed Z-stacks of 250–350 Z-planes with a field size of 91.18 × 91.18 μm (100× oil objective, NA 1.44, zoom 1.7 with a 0.09 μm step size). Therefore, these Z-stacks represented bone volumes of 187,060 to 261,884 μm3 (0.187 to 0.262 mm3). DAPI signal was collected in a separate scan using 405 nm excitation and a collection window of 410–480 nm. For DiO/Texas Red-Dextran imaging the DiO was excited at 488 nm with a collection window of 494–564 nm and Texas Red was excited in the same scan at 594 nm with a 604-684 collection window. To determine the lacunar volumes, RGB image stacks of triple merged Texas Red-Dextran, DiO and DAPI images were converted to 8 bit greyscale images to create a “filled in” image of the lacuna using Image J software (Rasband, W.S., ImageJ, U. S. National Institutes of Health, Bethesda, MD, http://imagej.nih.gov/ij/, 1997–2015). The 8 bit stacks were normalized to correct for signal loss with tissue depth and then thresholded in the Object Counter 3D plugin (Bolte and Cordelieres, 2006) so that just the lacunae were included. A lower volume limit of 10 μm3 was used.

In MIMICS®, segmentation of each image was histogram based, using values in Hounsfield Units (HU). Then using operations such as region-growing, gap-filling, pixel closing, wrapping, and Boolean subtraction, a working 3D 17 lacunae femur model was generated composed of bone and lacunar spaces. Pixel closing was required because regions represented by the pixels must be continuous with the next slice to avoid discontinuity errors in 3D volume meshing. Wrapping was used sparingly to smooth the surfaces of the lacunae reducing effects of pixilation and developing a less complicated model for computation. Finally, two solid model STL files (one bone matrix and one lacunae, with the bone matrix created by subtracting the lacunae from a solid block) were generated for each sample, and were imported to another software program – 3-Matic® - for meshing and finite element purposes. The STL files were imported from MIMICS to 3-Matic, repaired, refined, and re-meshed to generate a non-manifold assembly, and a volume mesh was finally generated using tetrahedral elements. A 3-D lacuna model, from a mouse specimen, containing seventeen lacunar bodies of various sizes is shown in Fig. 3A. Theoretical strain ratios (i.e. lacunae/bone) were calculated based on a load oriented longitudinally to the bone for ten lacunar bodies shown in Fig. 3B. Only fully intact lacunae were used for the analysis, therefore six lacunar bodies (O1, O4, O8, O10, O13, and O17) were not included because they were cropped off in the image field and thus were incomplete. A seventh lacunar body (O14) was also eliminated due to its sharp geometry. The volume mesh generated by 3-Matic comprised of 137,736 nodes and 530,205 4-node tetrahedral elements. Although bone is heterogeneous with varying material properties, it was treated as a homogeneous elastic material in this study. An elastic modulus of 25 GPa was assigned to the bone region based on values reported in the literature and previous work (Guess et al., 2010; Rho et al., 1998). Poisson's ratio was set at 0.3.

Fig. 3

3-D FE nano-model with seventeen lacunae created from a confocal image stack of mouse femoral bone. The load is applied in the longitudinal direction, which is also the longitudinal axis of the bone. (A) Front surface view showing mesh. One incomplete lacuna is lying in the cut plane (purple). (B) With outer surfaces hidden, the realistic lacunar shapes are revealed.

In addition, one model was created from the mid-shaft of a femur of a 5-month old female mouse to analyze for any unanticipated effects on results due to model mesh sizes. This model is termed the 7 lacunae femur model. Normalized lacunae/bone strain measurements, representing the strain magnification ratios, were calculated using the average strain in the lacunae and the bone, taking care to avoid elements influenced by edge effects. Four different mesh sizes from very coarse to fine, were created using procedures described in an earlier section. Loading, boundary conditions and material properties used were also similar.

Results

This section presents the results from the analyses on the different FE models described in this study. Data from the Micro model followed by the Nano model and finally the 3-D Lacuna model are presented.

Micro model

Strain amplification in lacunae

The first parametric analysis with the Micro model was performed to determine the theoretical strain magnification in various lacunae throughout the bone matrix. Table 2 provides an analysis summary and also includes other information such as the Young's modulus of the perilacunar tissue and the ratio of the lacunar strain to the global bone strain (i.e., L/B Ratio).

Table 2

Comparison of theoretical lacunar strains (με) and strain magnification (lacuna/bone or L/B ratio) in Micro model.

	Lacuna L1		Lacuna L2		Lacuna L3		Lacuna L4
Young's Modulus (Perilacunar Tissue)	Lacunar Strain	L/B Ratio	Lacunar Strain	L/B Ratio	Lacunar Strain	L/B Ratio	Lacunar Strain	L/B Ratio
5 GPa	9325	4.66	9167	4.58	9012	4.51	9000	4.50
10 GPa	8392	4.20	8168	4.08	7763	3.88	7713	3.86
15 GPa	7709	3.85	7390	3.70	6869	3.43	6798	3.40
20 GPa	7130	3.57	6760	3.38	6187	3.09	6102	3.05

A key observation from this iteration of the model was that in all cases the lacunar strain was amplified at least 3-fold with respect to the global bone strain, confirming that osteocyte lacunae act as stress concentrators. In all lacunae examined, the L/B strain ratio decreased when the perilacunar modulus increased. For example, in Lacuna L1, the L/B strain ratio decreased from 4.66 to 3.57 as the modulus increased from 5GPa to 20GPa. These results are similar to a previous study by Bonivtch et al. (2007), who varied the perilacunar modulus in a single osteocyte lacunar model, to investigate the magnification in strains due to mechanical loading. We have extended that work by including multiple lacunae in the bone matrix. For example, with a perilacunar modulus of 5 GPa, which may correspond to a perilacunar region that has more calcium removed, the strain magnification ratio ranged from 4.66 to 4.50 for lacuna (L1) and lacuna (L4). The strains were also lower in the osteocytes farther away from the loading surface.

Strain variation based on orientation

In this analysis, the four lacunae under analysis were named according to their orientation as Lacuna 1AD CE, 1BD, 2V and 2H. Fig. 4 shows the relationship between the elastic modulus of the perilacunar region (x-axis) and the maximum Lagrangian strain of the lacuna (y-axis). The perilacunar modulus was varied from 5 GPa to 20 GPa and the theoretical lacunar strains in each osteocyte were computed. From Fig. 4, it is evident that an increase in the elastic modulus of the perilacunar region from 5 GPa to 20 GPa resulted in decreased maximum lacunar strains for all the lacunae, regardless of their orientation. For instance, the plot for Lacuna 2H exhibited a decrease in lacunar strain from 17,027 με to 11,421 με. Similar trends were observed in the other lacunae.

Fig. 4

Comparison of theoretical lacunar strains as a function of perilacunar modulus in osteocyte lacunae with different orientations.

Lacuna 2H, which was oriented horizontally parallel to the loading surface, experienced the highest strains when compared to all other lacunae whereas Lacuna 2V, which was oriented vertically perpendicular to the loading surface experienced the lowest strains. Depending on the perilacunar modulus assigned, the horizontally aligned lacuna (i.e., 2H) generated strains that were higher by factors of 2.42 to 2.57 compared to the strains from the vertically aligned lacuna (i.e., 2V). As expected, the two diagonally-oriented lacunae exhibited similar strains although a slightly larger difference (6.3%) was observed for a perilacunar modulus of 5 GPa versus 20 GPa.

Strain variation based on size

The second parameteric analysis of the Micro model investigated the lacunar strain based on variations in the lacunar size. The three lacunae selected for this analysis were labeled as lacuna S1, S2 and S3 from smallest to largest. Fig. 5 shows the lacunar strain based on their size and the perilacunar modulus assigned.

Fig. 5

Comparison of lacunar strains as a function of perilacunar modulus in lacunae of different size.

At a given perilacunar modulus (e.g., 5 GPa), lacuna S3 generated the highest strain followed by lacuna S2 and finally lacuna S1. This data suggests that an increase in size of the lacuna resulted in an increase in the lacunar strain value. Although lacuna S3 consistently displayed the highest strain, a relatively small difference (1.46%) was observed between lacuna S3 and S2 at the lowest perilacunar modulus of 5 GPa. Meanwhile, the strain differences between lacuna S3 and S2 varied by 9.99%, 17.0%, and 23.1% for higher perilacunar moduli of 10, 15, and 20 GPa, respectively. In addition, the lacunar strain decreased as the perilacunar modulus increased from 5 GPa to 20 GPa. For lacuna S1, the strain value decreased from 7319 με at 5 GPa to 4021 με at 20 GPa.

2-D Nano model

Results from the Nano model analysis are presented as the theoretical strain ratios of the cell body, fluid space, and perilacunar region to the surrounding bone and are summarized in Table 3, Table 4, Table 5, respectively. For all six osteocytes and all three osteocyte subregions, the strain ratios decreased as the perilacunar modulus increased. Changing the lacunar orientation from horizontal to vertical (e.g. Osteocyte H1 vs Osteocyte V1) was found to decrease strains in the cell body (~82%) and fluid space (~76%) regions but not the perilacunar tissue. It can also be observed that the diagonally oriented (D1 and D2) cells behave similarly to the horizontal cells (H1 and H2) in terms of the cell body strains, they are “in between” the horizontal and vertical cells for the fluid space strains and there is not much effect of the cell orientation on the perilacunar strains. For the osteocytes oriented horizontally and diagonally the strain magnifications in the cell bodies were higher compared to those in the fluid space. In the case of vertical osteocytes, strains measured in the fluid space were higher than the cell body.

Table 3

Maximum cell body region strain to global bone strain (2000 με) for each osteocyte in the Nano model.

Perilacunar modulus	Osteocyte ID
	H1	H2	V1	V2	D1	D2
5 GPa	10.31	9.25	1.92	1.70	9.98	10.86
10 GPa	9.77	8.86	1.80	1.60	9.39	10.16
15 GPa	9.35	8.56	1.69	1.51	8.94	9.64
20 GPa	9.02	8.31	1.59	1.43	8.57	9.21

Table 4

Maximum fluid space region strain to global bone strain (2000 με) for each osteocyte in the Nano model.

Perilacunar modulus	Osteocyte ID
	H1	H2	V1	V2	D1	D2
5 GPa	8.67	7.78	2.11	1.98	6.18	6.66
10 GPa	8.08	7.32	1.96	1.85	5.70	6.11
15 GPa	7.63	6.98	1.85	1.74	5.34	5.71
20 GPa	7.28	6.70	1.75	1.65	5.06	5.40

Table 5

Maximum perilacunar strain to global bone strain (2000 με) for each osteocyte in the Nano model.

Perilacunar modulus	Osteocyte ID
	H1	H2	V1	V2	D1	D2
5 GPa	1.27	1.11	1.53	1.42	1.10	1.14
10 GPa	1.10	0.99	1.41	1.30	0.96	1.30
15 GPa	0.98	0.87	1.31	1.22	0.86	1.22
20 GPa	0.89	0.79	1.22	1.14	0.78	1.14

From the 17 lacunae 3D lacuna model developed by using confocal imaging data from a 5-month old female mouse femur it was observed that all the lacunae were oriented similarly along a direction slightly diagonal to the longitudinal axis of the bone. Hence the analysis of the 3-D lacuna model was designed to observe the strain differences due to variation in lacunar size. Table 6 shows the strain responses of the ten complete lacunar bodies selected for analysis and is organized by decreasing lacunar size. The 3-D modeling using geometries from real lacunae again confirmed that there was strain amplification at the lacunae, with the lacunar strains 1.7–3.5 fold higher than the global bone strain.

Table 6

Lacunar strains in 3-D lacuna model arranged by decreasing lacunar size.

Lacuna ID	No. of elements	Max strain (μƐ)	Strain magnification ratio
O16	876	7093	3.55
O5	705	5821	2.91
O6	467	5319	2.66
O12	458	4691	2.35
O2	447	5681	2.84
O15	425	4221	2.11
O11	389	4124	2.06
O9	354	4266	2.13
O7	310	4042	2.02
O3	278	3481	1.74

From Table 6, it can be observed that the maximum strain values and the strain amplification ratios for the 3-D lacunar model are similar to the values calculated in the nano and Micro models with a vertical lacunar orientation. The strain amplification in the lacunae compared to global strain are similar. There is also generally an increase in lacunar strain with increased lacunar size (any variations to this could be due to differences in the lacunar orientation or position).

Mesh size effects

Effects of mesh size on model results were evaluated by analyzing another 3-D seven lacuna femur model, shown in Fig. 6 which maps the third Principal Lagrange Strain in the lacunar regions. The figure shows the bone as a boxed region and the lacunae. The load is applied in the negative z region at the top and the bottom is fixed. It can be seen that the strain in the lacunae are non-uniform with the maximum negative strains toward the bottom.

Fig. 6

Color maps of the third Principal Lagrange Strain values for lacunae in the mid-shaft of the 7 lacunae female murine femur model. The loads were applied to the top surface of the model along the z-axis, in the negative z-direction, which corresponds to the longitudinal axis. The bottom surface of the model was fixed.

Four different mesh sizes were analyzed and Table 7 shows the details of the node and element numbers.

Table 7

Node and element numbers for various mesh sizes.

	Very Fine	Fine	Medium	Coarse
Nodes	985,764	538,168	22,517	7714
Elements	3,958,585	2,117,321	91,069	35,343

Fig. 7, compares the normalized strain ratios, defined as the ratio of the average strain in the lacuna to the strain in the central bone region, in each osteocyte for each mesh size. The figure shows that variations in mesh density does have a small effect on the strain magnification ratio with the ratio increasing as the mesh size becomes smaller. However, the magnification ratios for fine and very fine meshes were similar indicating that fine meshes are required for this type of simulation. The L/B ratios for these osteocytes were in the same range as results previously discussed.

Fig. 7

Comparison of normalized strain values in individual osteocytes in Model shown in Fig. 6 shown from coarse to a very fine mesh.

Discussion

Our previous study (Lara-Castillo et al., 2015a) demonstrated that in vivo loading of the mouse bone resulted in a very rapid, yet heterogeneous activation of Wnt/β-catenin signaling in a subset of osteocytes. In many instances, activated osteocytes resided adjacent to non-activated osteocytes. This suggests there is a “strain threshold” for osteocyte activation. The heterogeneity of the activation response suggests that other factors contribute to the determination of which osteocytes activate and that the strain microenvironment varies on a microscale that is not accounted for in traditional finite element models of global strain fields in bone (Lu et al., 2012; Fan et al., 2016; Moustafa et al., 2012). The exact mechanism by which mechanical load is perceived by the osteocyte remains unknown and is widely debated in the bone field. Our findings illustrate that the material properties of the bone surrounding the osteocyte lacunae as well as the size and orientation of the lacunae relative to the direction of applied load can have significant effects on the magnitude of strain experienced by the osteocyte.

As a first attempt to better determine strains experienced by osteocytes, we built finite element models that incorporate multiple osteocyte lacunae rather than assuming a uniform bone tissue. These models were used to test the effects of osteocyte lacunar orientation and size on the strain experienced by osteocytes. We built three models with increasing degrees of sophistication. First, we built a Micro model, which tested the effects of different idealized sizes and orientations on strain fields. Then we constructed a Nano-model in which appropriate osteocyte dimensions were used and also incorporated the osteocyte cell body and lacunar fluid compartments. Finally we built 3-D Lacuna models generated using empirical data from confocal image stacks from mouse femurs to determine the strain in lacunae represented by realistic in situ geometries. All three models gave internally consistent results and they confirmed the previous results of Bonivtch et al. (2007) showing that osteocyte lacunae act as strain amplifiers and that the lacunar strain is affected by the modulus of the perilacunar matrix. Using models that included multiple osteocytes rather than a single osteocyte, our study also extended those findings by demonstrating that the lacunar (osteocyte) orientation and size also affect strain amplification.

The perilacunar matrix is a region extending out from the edge of the lacunar wall about 5–10 μm. The material properties and hence the modulus of this region differs from the rest of the bone matrix (Nicolella et al., 2008). Verbruggen et al. (2012) found that the absence of this region in an idealized (i.e., ellipsoid) osteocyte FE model leads to poor strain magnification. The perilacunar matrix region therefore appears to be vital for mechanosensing and mechanotransduction since it modulates the strain environment experienced by the osteocyte cell membrane. Our FE results indicate that strains in the osteocyte lacunae decrease in a predominantly linear fashion with respect to an increasing perilacunar modulus and that this holds true regardless of the lacunar size and orientation. This also implies that small changes in the perilacunar matrix properties, for example due to mineral removal during lactation (Qing et al., 2012), could lead to large changes in the strain amplification.

Interestingly, of all the parameters tested, the orientation of the lacuna with respect to the applied load had the most significant effect on strain magnification. The horizontally aligned lacunae had the highest strain magnification ratio while vertically aligned lacunae had the lowest strain magnification ratio. Regardless of lacunar orientation, the strain magnification ratio still decreased as the perilacunar modulus increased from 5 GPa to 20 GPa. Based on the data in Table 3, Table 4, Table 5 the average decrease was 14% in the cell body, 15% in the fluid space region and 25% in the perilacunar space. The decreased strain magnification ratio indicates a stress shielding effect, wherein an increased perilacunar modulus (hence increased stiffness) results in a decreased lacunar strain. While the strain fields present in the bone material are routinely measured, it is actually the strain experienced by the osteocyte within its lacuna that will ultimately determine if its function is altered. As a first step in attempting to measure strain at the level of the osteocyte itself, which is critical for determining a strain threshold for activation, we created a nanoscale model that incorporated known osteocyte lacuna dimensions. The strain magnifications were analyzed at two locations, namely the fluid space and cell body. For the osteocytes oriented horizontally and diagonally the strain magnifications in the cell bodies were higher compared to those in the fluid space. In the case of vertical osteocytes, strains measured in the fluid space were higher than the cell body. This observation is particularly interesting since osteocytes in cortical bone of the mouse femur tend to have their long axis oriented parallel with the long axis of the bone (Sugawara et al., 2013; Vatsa et al., 2008), which would tend to reduce the increased strain perceived by the cell and result in strains that are higher in the fluid space compared to the cell body. It would also potentially allow the actin rich cytoskeleton within the cell body to align along the longitudinal direction and help to transmit the strains experienced by the cell membrane to the nucleus to affect pathway activation. In both the fluid space and cell body regions, the strain magnification ratio was highest for horizontally oriented osteocytes and lowest for vertically oriented osteocytes. Diagonally oriented osteocytes had strain ratios ranging between those of the horizontally and vertically aligned osteocytes.

Interestingly, we observed that most of the lacunae in the geometrically realistic 3-D Lacuna model that was based on empirical data were oriented vertically and/or diagonally. This is possibly to reduce the strain concentrations as observed in our theoretical models. The strain ratio range of 1.74–2.91 for nine of the ten osteocytes within the 3-D lacuna model coincided with the strain ratios observed for vertically oriented lacunae in our Micro and Nano models. Studies have provided experimental evidence that the alignment of osteocytes parallel to the long axis of the bone may be an adaptive response to mechanical loading to minimize the strains the osteocytes experience. Using adult mice, Vatsa et al. (2008) showed that osteocytes harvested from the fibula were aligned parallel to the primary axis of loading whereas the osteocytes from the calvaria, which is loaded at lower intensities, showed no preferential alignment. Furthermore, Sugawara et al. (2013) showed that newborn mice subjected to hind limb unloading for 6 weeks due to sciatic neurectomy, which is a model for unloaded long bone, formed cortical osteocytes with a much more random alignment and a more rounded morphology compared to normal ambulatory mice.

Another key observation from this study was that the lacunar strain value increased when the size of the lacunae increased. Regardless of lacunar size, the strain magnification ratio decreased with increasing values for the perilacunar modulus. Our finding that lacunar size alters the strain experienced by the osteocyte for a given load input may be important because lacunar size has been linked previously with bone pathology by van Hove et al. (2009). These investigators found that osteocyte morphologies differed among osteopenic, osteopetrotic, and osteoarthritic bone. Osteoarthritic bone contained small and elongated lacunae while osteopetrotic bone contained small and disc-shaped lacunae. Compared to the two previous bone types, osteopenic bone possessed the largest lacunae which may have been formed to counteract the decreased mechanosensitivity and progressive bone loss.

Our ultimate goal for developing a 3-D lacuna model is to have the capability to simulate in vivo conditions where realistic osteocyte geometries and heterogeneous bone properties are likely to produce varying strains in the lacunae. Verbruggen et al. (2012) demonstrated that FE models created with realistic dimensions from confocal images generated osteocyte strain magnifications that were 3.5 to 4 times larger than with 3-D FE models with idealized ellipsoid geometries. Furthermore, the strain generated from confocal image-derived models typically exceeded 10,000 με. Our 3-D lacuna model exhibited an increase in strain magnification with an increase in lacunar size which was consistent with both of the idealized models. Interestingly, the strain distribution in each lacuna was non uniform as seen in Fig. 6.

One limitation of this work is the simulation of bone as an isotropic material versus an anisotropic material although previous investigators have also followed this assumption in their FE studies (Bonivtch et al., 2007; Verbruggen et al., 2012; Verbruggen et al., 2014; Joukar et al., 2016). Measurement of local material anisotropic properties is a challenge and can influence the strain magnification factors shown in this study. Another limitation worth noting is that we did not vary the size of the lacunae in the Nano model. Previous studies have shown the interaction among the morphology, mechanical properties, and function of a cell morphology (Bacabac et al., 2008). Furthermore, the Nano model incorporated the lacunae as ellipsoids even though they may possess an irregular geometry as evidenced by the 3-D Lacuna model. Also, the current version of the 3-D Lacuna model does not include dendritic processes or heterogeneous bone properties, which would likely influence the strain magnifications already observed, and does not include vascular channels that exist within bone. Nonetheless, we have begun to provide a framework for building finite element models that can be merged with biological data from our TOPGAL mouse (a transgenic mice with a reporter line to visualize activation of the Wnt canonical pathway by beta-galactosidase activity) loading studies (Lara-Castillo et al., 2015a). This will enable us to better understand how load and generated strain can alter the biological activity of osteocytes in bone. Holguin et al. (2016) recently demonstrated that in the TOPGAL mouse model, which carries a lacZ reporter to indicate β-catenin activation, aged mice do not activate the β-catenin pathway as effectively compared to young adult mice.

Conclusions

This work is the first step in the development of a multiscale FE modeling approach for determining the strain threshold for activation of osteocytes in response to load. The size, orientation, and shape of the osteocyte lacunae were analyzed to investigate other potential factors that could explain the differential Wnt/β-catenin activation observed in previous mouse forearm loading experiments by Lara-Castillo et al., 2015a, Lara-Castillo et al., 2015b. Results from the nano FE models, which included multiple osteocytes, confirmed that the maximum strains in the osteocyte cell body, fluid space, and perilacunar tissue all increased due to the increased size of the lacunae. In addition, strains in all three osteocyte subregions were influenced by the orientation of the lacunae with the horizontal lacunae having the highest strain magnification ratio and vertical orientations the least. This study confirms that osteocyte size, orientation, and geometry all affect the strain magnification; however, the extent of these effects as well as the interaction of these parameters warrant further investigation through a combination of FE models (i.e., both idealized and realistic geometries). A better understanding of how mechanical stimuli directly affects the cell body, fluid space, and perilacunar tissue may provide evidence that can be used to identify the minimum strain threshold required for osteocyte activation. Finally, by constructing more models and looking at more regions of interest within bone, we can study changes in bone properties and their effect on bone responsiveness to loading. For example, it is well established that as we age the response to loading is attenuated (Burr, 1997; Pearson and Lieberman, 2004). Understanding how aging affects bone properties and incorporating those changes into our models will enable us to better predict/understand why the aging process results in diminished bone formation responses to loading (e.g., exercise).

Author contributions

GT and MLJ were responsible for the conceptual design of the studies. SKK developed the microscale and nanoscale models and performed analyses under the guidance of GT. MTB and LEL developed the 3-D lacuna models using confocal images provided by LMTL and SLD. All authors contributed toward the writing of this manuscript.

CRediT authorship contribution statement

Sravan K. Kola: Conceptualization, Funding acquisition, Supervision, Writing - original draft. Mark T. Begonia: Conceptualization, Funding acquisition, Supervision, Writing - original draft. LeAnn M. Tiede-Lewis: Conceptualization, Funding acquisition, Supervision, Writing - original draft. Loretta E. Laughrey: Conceptualization, Funding acquisition, Supervision, Writing - original draft. Sarah L. Dallas: Conceptualization, Funding acquisition, Supervision, Writing - original draft. Mark L. Johnson: Conceptualization, Funding acquisition, Supervision, Writing - original draft. Thiagarajan Ganesh: Conceptualization, Funding acquisition, Supervision, Writing - original draft.