Effect of osteoporosis treatment agents on the cortical bone osteocyte microenvironment in adult estrogen-deficient, osteopenic rats.

Though osteoporosis is a significant cause of disability worldwide, treatment with pharmacologic agents decreases risk of fragility fracture. Though these treatments act through the bone remodeling system to improve bone mass, it is unclear if they alter the response of bone to mechanical loading at the level of the osteocyte. This pre-clinical study determined the relationship between microstructural bone tissue properties and osteocyte lacunar size and density to strain around osteocytes with standard osteoporosis treatment or sequential therapies. Six-month-old female ovariectomized (OVX) Sprague-Dawley rats were cycled through various sequences of pharmacological treatments [alendronate (Aln), raloxifene (Ral) and human parathyroid hormone-1,34 (PTH)] for three month intervals, over nine months. Linear nanoindentation mapping was used to determine Young's modulus in perilacunar and bone matrix regions around cortical bone osteocyte lacunae. Measurements of lacunar diameter and density were completed. Treatment-related differences in Young's modulus in the perilacunar and bone matrix regions were not observed. We confirmed previous data that showed that the bone matrix region was stiffer than the perilacunar matrix region. Whole bone material properties were correlated to perilacunar matrix stiffness. Finite element models predicted a range of mechanical strain amplification factors estimated at the osteocyte across treatment groups. In summary, though the perilacunar matrix near cortical osteocyte lacuna is not as stiff as bone matrix further away, osteoporosis treatment agents do not affect the stiffness of bone tissue near osteocyte lacunae.

Introduction

Osteoporosis is a prevalent musculoskeletal disease that causes significant disability and death throughout the world (Salomon et al., 2012; Wang et al., 2012). Approved medicinal treatments for osteoporosis currently include anti-resorptive agents (e.g., bisphosphonates [BPs], Denosumab, and Raloxifene) and an anabolic agent (hPTH (1-34). Anti-resorptive agents generally improve bone mineral density (BMD) and reduce fracture risk by both reducing bone turnover and increasing the duration of the secondary phase of bone mineralization. Preclinical studies of all these agents generally show that their greatest effects are on trabecular bone (Evans et al., 1996; Ominsky et al., 2008; Seedor et al., 1991; Allen et al., 2006). Similarly, clinical studies with anti-resorptive agents find up to 70% fracture risk reduction in trabecular bone sites like the spine and 50% reduction in cortical bone sites like the hip (Black et al., 1996; Black et al., 2007; Cummings et al., 2002; Cummings et al., 2009; Liberman et al., 1995; Neer et al., 2001). While hPTH (1-34), the anabolic agent, also reduces risk of spine fracture by 65%, it increases trabecular and cortical bone turnover rate, with no significant effect on hip fracture risk and only modest positive effects on cortical bone mass (Neer et al., 2001; Zanchetta et al., 2003; Ma et al., 2014; Borggrefe et al., 2010; Eriksen et al., 2014).

Bone active agents may influence bone modeling and remodeling indirectly by altering the ability of osteocytes to sense the mechanical load signals that are transmitted to bone. Structural properties of the bone tissue surrounding osteocyte lacunae may influence the degree of mechanical stimulation to which an osteocyte is exposed (Eswaran et al., 2006; Kinney et al., 2000). Past work suggests that the composition of the bone tissue in the perilacunar region adjacent to osteocyte lacunae differs from that further from osteocyte lacunae (Meers et al., 2013; Meers, 2013). This difference may not only affect the mechanical loading signals received by the osteocyte (Stern and Nicolella, 2013), but may also be a sign of the osteocyte's known ability to modify its nearby surrounding matrix (Parfitt, 1977; Qing et al., 2012; Sharma et al., 2012). For example, stiffer perilacunar tissue may attenuate mechanical signals received by osteocytes (Bonivtch et al., 2007). Thus, under normal skeletal loading conditions, osteocytes surrounded by perilacunar bone tissue that has been stiffened by a condition such as estrogen deficiency, would receive lower intensity, and perhaps fewer, mechanical signals that alters bone modeling and remodeling activity to favor increased bone resorption and bone loss (Frost, 1997). Rats treated continuously with hPTH (1-34) had larger osteocyte lacunae due to perilacunar bone tissue resorption by osteocytes, with evidence of regeneration of perilacunar tissue (Tazawa et al., 2004). Osteocytes respond quickly to low doses of calcitonin by depositing mineral aggregates near the lacunar wall (Matthews and Talmage, 1981; Norimatsu et al., 1979). Others have shown that osteocyte lacunae change with loading (Bozal et al., 2001; Ferreyra et al., 2000) and that matrix metalloproteinase (e.g., MMP-13) concentration in the osteocyte lacuna increases during lactation (Qing et al., 2012). The osteocyte appears to alter perilacunar bone in response to metabolic demands (Matthews et al., 1978) and mechanical signals.

Osteoporosis treatment agents may alter cortical bone matrix. This can be inferred from measurements of ultimate stress (Ardawi et al., 2016; Wang et al., 2011; Bahar et al., 2016; Takeda et al., 2017), a whole bone biomechanical strength endpoint that indirectly reflects bone material properties including those of bone matrix, and more direct methods such as nanoindentation testing (Yang et al., 2013; Brennan et al., 2009a) and FTIRI (Fourier transform infrared imaging) (Gamsjaeger et al., 2011). Separate studies report that osteoporosis treatment agents can also alter the bone response to mechanical loading (Roberts et al., 2009; Hagino et al., 2001). Since patients being treated clinically for osteoporosis are now routinely exposed to sequential treatment periods that involve anti-resorptives, cessation of treatment, and anabolic agents, additional knowledge of how such sequences influence the perilacunar tissue properties and thereby the localized strain, would be helpful (Yao et al., 2007; Shahnazari et al., 2010). Therefore, the purpose of this study was to assess: a) treatment-related differences in Young's modulus of the perilacunar and bone matrix regions near osteocyte lacunae and lacunar major and minor diameters; b) the relationship between perilacunar and bone matrix tissue properties and whole bone strength and material properties; in estrogen-deficient, osteopenic adult rats given three sequential three month treatments with common, approved osteoporosis treatment medications; and c) the implications of any differences in perilacunar tissue properties on the mechanical signals received by osteocytes, using finite element modeling (FEM). We hypothesized that estrogen deficient rats have higher Young's modulus in the perilacunar matrix around osteocytes with lower diameter and density of osteocyte lacunae than Sham rats. We also hypothesized that rats given sequential administration of medications that improve bone mass and decrease fracture risk in humans have lower Young's modulus in the perilacunar matrix and greater diameter and density of osteocyte lacunae than untreated OVX rats. It is well known that osteoporosis treatment medications influence bone remodeling activity, and that bone remodeling rate is a major determinant of tissue age. In order to limit the influence of bone remodeling activity and bone age on local bone material properties, we analyzed the middle third of the long bone cortex, a region whose bone tissue lacks Haversian remodeling in rats (Frost and Jee, 1992).

Materials and methods

Animals and experimental procedures

The bones studied here are from a subset of rats that were treated for nine months in a previously-reported experiment (Amugongo et al., 2014a; Amugongo et al., 2014b). Six-month-old female ovariectomized (OVX) or sham-OVX (Sham) Sprague-Dawley rats were purchased from and operated on at Harlan Laboratories (Livermore, CA; USA). After a one week acclimation period to our laboratory, pair-feeding of OVX to Sham rats was initiated. All OVX rats were allowed to develop trabecular osteopenia of the distal femur and lumbar vertebral body for eight weeks post-surgery (Amugongo et al., 2014a), then randomized by body weight into seven groups that represent both currently-applied and potential sequences of anti-osteoporosis medications (see below). OVX rats were treated for three 90 day periods with vehicle [Veh, normal saline, 1 ml/kg/dose, three times per week by subcutaneous (SC) injection] (Life Technologies, Cat# 10010; Grand Island, NY, USA); parathyroid hormone, [hPTH (1–34) (human) acetate 25 μg/kg/dose, 5×/wk. SC] (Bachem Biosciences Inc., Cat# H-4835; King of Prussia, PA, USA); alendronate [Aln, 25 μg/kg, 2×/wk. SC] (Sigma, Cat# A-4978; St. Louis, MO, USA); or raloxifene [Ral, 5 mg/kg/dose, 3×/wk. by oral gavage] (Sigma, Cat# R-1402), in sequences described below. Each treatment period involved only monotherapy. No two treatments were ever given at the same time, meaning that combination treatment never occurred in this experiment. The study protocol was approved by the University of California Davis Institutional Animal Care and Use Committee.

After 90 days (end of Period 1), the animals switched to their second treatment regimen. After a total of 180 days (end of Period 2), the animals switched to their third treatment regimen. After a total of 270 days (end of Period 3), all rats were necropsied. The sequential treatments applied were Sham, Veh-Veh-Veh, Aln-Aln-Aln, Ral-Ral-Ral, PTH-Veh-Veh, Aln-Veh-Aln, Aln-PTH-Aln, and Ral-PTH-Ral (Amugongo et al., 2014a; Amugongo et al., 2014b). The central right tibia from these rats was previously tested in three-point bending (Amugongo et al., 2014b). Though cortical bone was unchanged at two months post-OVX, Veh-Veh-Veh rats had lower cortical area and cortical thickness at the central femur than Sham rats by the end of the nine month treatment period (Amugongo et al., 2014b). Aln-Aln-Aln, Aln-Veh-Aln, Aln-PTH-Aln, and Ral-PTH-Ral groups had significantly greater distal femoral BV/TV than Veh-Veh-Veh rats after nine months treatment (Amugongo et al., 2014a). Aln-Aln-Aln, Aln-Veh-Aln, Aln-PTH-Aln, and Ral-PTH-Ral groups had significantly greater central femur cortical area than OVX rats after nine months.

At necropsy, left tibiae were dissected free, wrapped in saline-soaked gauze, and frozen at −20 °C. With eight groups (seven primary comparisons to the OVX group) and N = 5/group, and reducing the alpha to 0.0071 to account for the seven comparisons, we have 80% power to detect a difference of 2.24 standard deviations between groups. Although a large effect size, given the precision with which Young's modulus is estimated, such a difference corresponds to ~10.5% difference in means. A recent article found a 24% difference in Young's modulus between groups when testing drug effects on Young's modulus (Nowak et al., 2016). Since we have power to detect a difference less than half their observed difference and a 10.5% difference in Young's modulus would be clinically meaningful, our study design with eight groups and N = 5/group was justified. Thus, we randomly selected left tibiae from five rats of each group. The cortical bone of the tibia was chosen because any effects of estrogen deficiency and the treatment agents on its bone matrix are highly likely to be mediated only through their actions on osteocytes, rather than on bone remodeling activity.

Polymethyl methacrylate (PMMA) embedding

All remaining soft tissues were removed from the thawed left tibiae. They were halved longitudinally in the parasagittal plane with a water-cooled low-speed diamond saw (Buehler, Lake Bluff, IL: USA). The resulting lateral and medial halves were placed in 70% ethanol for 24 h, and then dehydrated in 95% (two solution changes) and 100% (four changes) ethanol each for at least 8 h in each change, a procedure that does not alter bone material properties (Bushby et al., 2004). Each specimen was then placed in a 20 ml glass scintillation vial (Wheaton; Millville, NJ) in an infiltration solution that contained 8.4 ml of unpolymerized methyl methacrylate (MMA, Sigma-Aldrich, Cat.# M55909, St. Louis, MO, USA), 1.4 ml of dibutylphthalate (Sigma-Aldrich, Cat.# D2270, St. Louis, MO, USA), 100 μL of polyethylene Glycol 400 (PEG400, Hampton Research, Aliso Viejo, CA), and 70 mg benzoyl peroxide (Polysciences, Inc., Warrington, PA, USA) at room temperature for 24 h. This solution was changed three times. On the fourth day, sufficient N,N-Dimethyl-p-toluidine (MP Biomedicals, Solon, OH), was added to make a 33 μM concentration of it in the embedding solution. Both medial and lateral tibial specimens were oriented with the trimmed aspect of the bone face down in the glass vial in the embedding solution at 4 °C for a week, until the PMMA block was completely hardened (Qing et al., 2012). After embedding, the exposed face of the bone was sequentially polished on POLIMET I polishers with diamond powder suspensions (Buehler; Lake Bluff, IL) of decreasing particle sizes starting at 6 μm and ending with a final size of 0.25 μm, to create a PMMA block with a smooth mirror-like face on the exposed bone surface for use in nanoindentation tests (lateral half) and scanning electron microscopy (medial half) (see below).

Nanoindentation

A nanoindenter (Triboscope, Hysitron Inc.; Minneapolis, MN, USA) attached to a Nanoscope IIIa atomic force microscope (AFM, Digital Instruments Inc.; Santa Barbara, CA, USA) with a calibrated diamond-tipped Berkovich indenter (tip diameter = 150 nm) was used. The loading direction was perpendicular to the polished surface of the exposed bone face in the PMMA block. The indenter tip was aligned in x-, y-, and z-axes by viewing with a microscope (Nikon SMZ800, Model P-IBSS, Japan) mounted to the indentation apparatus. The polished block containing the exposed face of the lateral half of the tibia was glued to a steel stub with cyanoacrylate adhesive and mounted on the microscope stage with polished face up. 100 μm × 100 μm AFM scanning area images were generated to identify osteocyte lacunae for testing. The lacunae were randomly selected from the middle one-third region (relative to the periosteal and endocortical surfaces) of the medial and lateral quadrants of cortex ~2 mm proximal to the tibia-fibula junction (TFJ).

A linear sequence of nanoindentation test sites located perpendicular to a lacunar wall was designated, starting at ~1 μm from the lacunar wall and proceeding in 1 μm increments away from the lacunar wall for 20 μm (Fig. 2B). Six such lacunar tracks were identified and studied per bone. Thus, on each track, five nanoindentation tests were performed at 1–5 μm from the lacunar wall and five tests were performed at 16–20 μm from the lacunar wall (Fig. 2B). Young's modulus is higher in bone matrix located >4–7 μm from the osteocyte lacunar wall than in bone matrix located adjacent to the lacunar wall (Meers, 2013; Stern and Nicolella, 2013). Each individual nanoindentation test was conducted at a peak force load of 1000 μN with a loading/unloading rate of 250 μN/s and holding segment time of 3 s. Based on the generated force-displacement curves from each individual test in the sequence of ten, the initial part of the unloading curve was analyzed using the Oliver–Pharr method to provide Young's modulus (E) values at each nanoindentation test site (Xu et al., 2012). The mean of all 1–5 μm nanoindentations for a bone (N = 30) was used for its perilacunar matrix Young's modulus value and the mean of all 16-20 μm nanoindentations (N = 30) was used for its bone matrix Young's modulus value.

Fig. 2

A. SEM image of an osteocyte lacuna from acid-etched resin-embedded cortical bone from the lateral quadrant of the tibial midshaft. At left, resin-filled lacunar void is completely visible, with multiple resin-filled canaliculi extending from the lacuna. Only lacunae completely exposed during the acid-etch process, as shown on the left, were used for measurement of lacunar diameters.

At the right, long (horizontal) and short (vertical) axis diameter measurement bars are superimposed on the same lacuna. Scale bars (10 μm) are in the bottom center of each image.

B. A 20 μm long linear nanoindentation path is shown in the upper right. The first crossing line on the path is 5 μm from the lacunar wall. The second crossing line is 15 μm from the lacunar wall. The end of the path is 20 μm from the lacunar wall. Five evenly-spaced measurements of Young's modulus were made between 0 and 5 μm from the lacunar wall and five more identical measurements were made between 16 and 20 μm from the lacunar wall. Scale bar (10 μm) is in the bottom center.

Scanning electron microscopy (SEM)

The PMMA-embedded medial half of the tibia was polished as above to expose the bone face, and then immersed for etching in 37% phosphoric acid once for 30 s, and then rinsed three times with distilled water for three minutes each time. The etched specimens were then immersed in 5% sodium hypochlorite for five minutes and rinsed three times with distilled water for three minutes each time. The specimens were dried overnight, mounted on aluminum stubs, and sputter coated with 20 nm of gold‑palladium (Lu et al., 2011). Specimens were examined using SEM (Philips XL30 ESEM-FEG 515; Philips Electron Optics Inc.; Hillsboro, OR) at 15 kV (Bozal et al., 2001).

The SEM was used to image completely surface-exposed osteocyte lacunae. Eight SEM images were taken from the central cortical bone 2 mm proximal to the TFJ, in both the medial and lateral region of each tibia at 400× magnification for quantitation of lacunar density (lacunae/mm2). Vascular channels and resorption cavities, that were excluded from all lacunar analysis procedures, were recognized by the depth with which they extended into the bone, that was always much greater than that of osteocyte lacunae, and their diameter, that was usually much larger than lacunae in both axes. Ten to eleven lacunae per specimen were randomly selected and images were captured at 3000× magnification for measurement of the length of the minor and major lacunar diameters (Fig. 2) (Bonivtch et al., 2007). ImageJ software was used to analyze the images (ImageJ 1.45 k, NIH, USA).

Finite element analysis (FEA)

Eight models, one for each treatment group, were created for finite element analysis (FEA), based on material properties and lacunar geometry data (Table 1 and Fig. 4).

Table 1

Data used to generate finite element models.

Variable/group	Lacunar diameter long (LDL) axis (μm)	Lacunar diameter short (LDS) axis (μm)	Young's modulus- perilacunar region (GPa)	Young's modulus- bone matrix region (GPa)	Lacunar density (mm−2)
Sham	12.81	6.06	27.33	31.44	589
OVX	13.24	6.10	27.55	29.23	785
Aln-Aln-Aln	14.39	5.68	28.62	30.23	788
Ral-Ral-Ral	13.12	5.73	27.93	31.23	839
PTH-Veh-Veh	14.27	6.45	28.74	29.72	691
Aln-Veh-Aln	13.82	5.84	26.35	27.51	818
Aln-PTH-Aln	13.46	6.20	27.31	29.67	777
Ral-PTH-Ral	13.79	6.03	28.80	30.09	872

Fig. 4

Fringe plots illustrating the strain concentrations that occur around the osteocyte lacuna. Note that each model is presented on the same scale, so differences in the above sizes are representations of varying lacunar density and size. The resulting first principle strain is plotted for an applied global physiological strain of 2000 microstrain.

The models were created using our parametric finite element model of the osteocyte lacuna (Bonivtch et al., 2007). Briefly, the osteocyte lacuna was modeled as an ellipsoid revolved around the major axis. The dimensions of the surrounding bone matrix region were calculated based on the lacunar density. The bone tissue within 20 μm of the lacunar wall was idealized as an isotropic linearly elastic material with two adjacent regions (Fig. 1A). The thickness of the perilacunar region (PT) was held constant at 5 μm. A Poisson ratio of 0.3 was used for each model (Bonivtch et al., 2007).

Fig. 1

(A) Cross-section of the parametric finite element model of an osteocyte lacuna and its two types of surrounding bone matrix (LDL- Lacunar Diameter Long; LDS- Lacunar Diameter Short; PT- Perilacunar Thickness). (B) An eighth of the model was modeled, meshed, and then analyzed during the simulations. (C) The prescribed boundary conditions effectively reflect the eighth of a lacuna model in three dimensions infinitely. Spacing is based upon the measured lacunar density.

The boundary conditions for the analysis included a forced axial displacement perpendicular to the long axis of the lacuna (vertically in Fig. 1A). The prescribed displacement for each model was individually calculated to result in an overall structural strain of 2000 microstrain. This was previously described as the upper bound of physiological strain in an investigation of in vitro strain during human rigorous physical activity and has been shown to scale to bones the size of rodents (Burr and Allen, 2013). The nodes on the model face opposing the prescribed displacement were fixed in all degrees of freedom, while the nodes on the perpendicular faces were constrained to allow in-plane motion. These constraints resulted in an aligned and uniform distribution of infinite osteocyte lacunae in all three directions (Fig. 1C).

A mesh density study was conducted to ensure that the numerical solution had indeed converged. Briefly, lacunar strain was measured as the maximum strain (first principle) in the model. Mesh convergence was reached when the maximum lacunar strain differed by <2% (Bonivtch et al., 2007). This same hexahedral mesh density was utilized in each of the eight FE models. To decrease processing time, the actual analysis was performed using one-eighth of each model and assuming a symmetrical shape (Fig. 1B).

The geometries and meshes for each model were created in and exported from TrueGrid (XYZ Scientific Applications, Livermore, CA, USA). The quasi-static simulation and post-processing were all executed using an LS-Dyna implicit solver followed by LS-Prepost (Livermore Software Technology Corporation, Livermore, CA, USA). Each simulation took approximately ten minutes on a single high performance CPU (OptixPlex 790 with 3.3GHz Dual Processors). The average maximum lacunar strain (first principle) was calculated for each model as well as the resulting strain amplification ratios. The strain amplification ratios were calculated as the average maximum lacunar strain experienced by each osteocyte lacuna modeled, divided by the applied macroscopic strain (2000 microstrain). Fringe plots were also prepared to visualize the strain contours and concentrations.

Statistical analysis

Two-factor ANOVA (treatment X location [perilacunar or bone matrix]) was completed for Young's modulus data. The difference between bone matrix and perilacunar matrix Young's modulus in each rat was calculated. Young's modulus, inter-region difference, and lacunar diameter datasets were subjected to one-way ANOVA to identify treatment effects with a Tukey's post-hoc test when the ANOVA was significant. Pearson's correlation coefficient was calculated for three point bending strength of the right central tibia (Amugongo et al., 2014b) vs. perilacunar and bone matrix Young's modulus for the specific bones measured, that represent all seven treatment groups. All values are expressed as mean ± SD. A P-value < 0.05 was considered significant. All statistical analyses were performed using SPSS statistical software (version 23.0; SPSS Inc., Chicago, IL, USA).

Results

Young's modulus (Nanoindentation)

There were no significant differences among treatment groups in the perilacunar matrix, the bone matrix, or the difference between the two regions, of cortical bone for Young's modulus (Table 1, Table 2). Two-factor ANOVA showed that Young's modulus of the perilacunar region was significantly less than that of the bone matrix region (P = 0.0004) with no difference across the different treatments. There was no interaction between matrix region and treatment (Table 2). Ultimate stress, work to failure, yield stress, and Young's modulus of the whole bone were significantly and negatively correlated to Young's modulus of the perilacunar region, but not to Young's modulus of the bone matrix region (Table 3).

Table 2

Young's modulus (GPa) of perilacunar matrix and bone matrix regions.

Group/variable	Region		Inter-region difference
	Perilacunar	Bone matrix	Bone-perilacunar
Sham	27.33 ± 1.50	31.44 ± 3.41	4.09 ± 2.09
OVX	27.55 ± 1.29	29.23 ± 1.11	1.68 ± 1.26
Aln-Aln-Aln	28.62 ± 1.56	30.23 ± 2.91	1.61 ± 2.63
Ral-Ral-Ral	27.93 ± 1.74	31.23 ± 1.77	3.30 ± 3.07
PTH-Veh-Veh	28.74 ± 1.65	29.72 ± 1.49	0.98 ± 2.31
Aln-Veh-Aln	26.35 ± 1.25	27.51 ± 1.31	1.16 ± 2.30
Aln-PTH-Aln	27.31 ± 2.89	29.67 ± 2.32	2.36 ± 1.56
Ral-PTH-Ral	28.80 ± 0.87	30.09 ± 2.99	1.30 ± 2.65

Mean ± SD.

Two factor ANOVA (P-values):

Region (P = 0.0004)

Treatment (P = 0.7682)

Interaction of Region X Treatment (P = 0.8464).

Inter-Region Differences ANOVA (P-value):

Treatment (P = 0.3209).

Table 3

Pearson's correlation coefficients (R) of whole bone strength endpoints in the right central tibia from rats of all seven treatment groups to perilacunar and bone matrix Young's modulus (by nanoindentation).

Variable	Young's modulus (perilacunar)		Young's modulus (bone matrix)
	R	P=	R	P=
Ultimate load (N)	−0.0704	0.7069	−0.0234	0.9007
Ultimate stress (N/mm2)	−0.3767	0.0367	−0.0243	0.8966
Work to failure (J)	−0.3761	0.0370	−0.2620	0.1544
Yield stress (N/mm2)	−0.5183	0.0028	−0.2088	0.2596
Young's modulus (whole bone)	−0.3728	0.0389	−0.0891	0.6335

Lacunar density and diameters (SEM)

Representative lacunar images of the Aln-Veh-Aln (left), Sham (middle), and Veh-Veh-Veh (right) groups are presented (Fig. 2, Fig. 3). There were no significant intergroup differences in either lacunar diameters of the two axes (Fig. 1) or lacunar density (Fig. 3).

Fig. 3

Osteocyte lacunar density (A). Mean ± SEM.

Mean values for Young's modulus of the perilacunar matrix and bone matrix region (Table 1) were used in the FE models to predict strain distributions around osteocytes (Fig. 4). A range of maximum lacunar strains was predicted by the FEA models across the treatment groups (Table 4). The highest predicted maximum lacunar strains were found in the Sham model, while the lowest strains were seen in the Aln-Aln-Aln and Aln-Veh-Aln groups. Of all the treatment groups, Aln-PTH-Aln had the highest estimated maximum osteocyte lacunar strains. Varying patterns of strain distribution across the osteocyte lacunae of the treatment groups were also observed (Fig. 4). As expected due to the geometry of the osteocyte lacuna, for all treatment groups, the area of highest strain was located at the center of the lacuna. In some treatment groups, the areas of high strain did not reach as far as in the Sham model (e.g., Aln-Aln-Aln) (Fig. 4).

Table 4

Maximum Strain amplification at the osteocyte lacunar wall determined by FEA.

Group	Maximum lacunar strain (microstrain)	Strain amplification factor (max strain/imposed strain)
Sham	2871	1.44
OVX	2762	1.38
Aln-Aln-Aln	2616	1.31
Ral-Ral-Ral	2750	1.38
PTH-Veh-Veh	2736	1.37
Aln-Veh-Aln	2665	1.33
Aln-PTH-Aln	2781	1.39
Ral-PTH-Ral	2695	1.35

Discussion

Bone active agents such as BPs, SERMs, and hPTH (1–34) used to treat osteoporosis have been evaluated for their efficacy in reducing risk of fragility fracture and preventing loss of and/or increasing hip and spine bone mineral density (BMD). BMD is currently the best measurable assessment of future fracture risk available for humans, usually showing relative risks of 2–3 (Stone et al., 2003; Black et al., 2017). However, a positive history for prevalent spine or hip fracture is the best single observable assessment of future fracture risk available for humans, usually showing relative risk of 5 or more (Wasnich et al., 1994; Ross et al., 1991; Lindsay et al., 2001; Black et al., 1999). This signals a gap in what can be learned from BMD about risk of future fracture and the actual risk of future fracture. Thus, most now feel that there is ample room for measurements that go beyond BMD to improve the understanding of how osteoporosis treatments affect bone quality, clarifying how such treatments reduce fracture risk. Bone quality, as generally defined, incorporates a number of endpoints such as bone microarchitecture/geometry, bone turnover rate, collagen properties, and tissue level material properties (Compston, 2006). At present it is not possible to perform many of these assessments in living humans. Animal models have been employed to bridge this gap.

Many studies have used animal models to evaluate the effects of these bone active agents on whole bone strength, stiffness, and architectural properties. BMD, biomechanical, and architectural properties as currently assessed reflect only macro-level properties of bone. At the tissue material or micro/nano-level, fewer studies have been published (Yang et al., 2013; Brennan et al., 2009a; Gamsjaeger et al., 2011; Aruwajoye et al., 2017; Mandair and Morris, 2015; Roschger et al., 2014; Buckley et al., 2014; Gamsjaeger et al., 2014; Donnelly et al., 2010). We found that as bone tissue in the perilacunar region of cortical bone stiffens, decreases in bone material properties that can be considered to reflect bone quality, such as ultimate stress, work to failure, and yield stress, occur. A current theory suggests that stiffer perilacunar matrix reduces the level of mechanical loading signals that are transmitted, thereby causing osteocytes to experience reduced loading and accommodate by increasing bone resorption in the modeling mode, reducing whole bone strength (Frost, 1997; Frost, 1987). We found that bone material properties, including ultimate stress, work to failure, yield stress, and Young's modulus of the whole bone, were significantly and negatively correlated to Young's modulus of the perilacunar region, but not to Young's modulus of the bone matrix region (Table 3). These data suggest that this theory has merit and should be tested by additional experimentation, perhaps involving osteocyte lacunae in trabecular regions, such as the long bone metaphysis or vertebral body, that have higher bone remodeling rates than the cortical region studied here. Additional studies of this type could yield a better understanding of how these bone properties reflect bone quality. Furthermore, a better understanding of how alterations in these various components alter strain distribution in bone at the local tissue level is needed, as fracture is ultimately a failure of bone tissue to manage local strain overload.

We examined the effects of three conventional treatment regimens and three sequential treatment regimens with those same bone active agents, on both cortical bone around osteocyte lacunae and dimensions of the lacunae themselves. We found no significant treatment-related differences in Young's modulus in cortical bone tissue around osteocyte lacunae. Two previous studies made randomly-placed nanoindentation measurements over a wide surface area of trabecular or cortical bone tissue from OVX rats given PTH, bisphosphonate, or raloxifene. Increased elastic modulus was observed in six month old rats that had been treated for two months with PTH or bisphosphonate (Yang et al., 2013). Increased elastic modulus was also seen in 12-month-old rats that had been treated for four months with raloxifene, but not with bisphosphonate (Brennan et al., 2009a). Rats in the current study were treated more than twice as long as either previous study. In addition, our sampling method was confined to a small region of bone matrix with a specific relationship to nearby osteocyte lacunae. Treatment duration and methodologic differences appear to make comparing our data to existing published work (Yang et al., 2013; Brennan et al., 2009a) difficult.

Even when no treatment-related differences in these variables are present, it is the dynamic interaction of these variables, and not their individual values, that determine strain distributions. Therefore, we used these data to construct finite element (FE) models of the strain distributions around the osteocyte that would occur with each combination of cortical bone properties and lacunar dimensions. As these osteoporosis treatment agents have different mechanisms of action, we hypothesized that they would have different effects on the strain to which osteocytes are exposed.

We observed no significant differences among the treatments in Young's modulus of either the perilacunar or bone matrix region around osteocytes in cortical bone. We also saw no significant effect of treatments on the difference in Young's modulus between the bone matrix and perilacunar region around osteocytes in cortical bone. These data indicate that osteoporosis drugs administered sequentially for up to three months in rats do not cause osteocytes to influence the local mechanical properties of bone tissue around them in cortical bone, in either a relative or absolute sense. We chose a region of rat cortical bone, because rats lack secondary osteonal remodeling of cortical bone (Frost and Jee, 1992), minimizing the chance that local bone properties measured here would be influenced by agent effects on bone remodeling activity, thereby providing us the most straightforward chance to examine local actions of osteocytes on nearby bone tissue. Furthermore, the positive effects of PTH and anti-resorptive agents on cortical bone in small animals and osteoporotic patients are less robust than those in trabecular bone (Arita et al., 2004; Wronski et al., 1993; Brouwers et al., 2009; Fox et al., 2006; Yano et al., 2014). A study of perilacunar matrix and bone matrix around osteocyte lacunae in trabecular bone would yield important complementary data to our study. The small differences in osteocyte lacunar shape and density among treatment groups also contributed to the altered strain fields in our models and are reflected in the size of the box in the strain models, a larger box indicating a lower lacunar density.

To date, there are few published studies with which to compare directly our findings. A single 16 week treatment regimen in 8-month-old OVX rats has been reported (Brennan et al., 2009b). Using nanoindentation they observed lower Young's modulus in pamidronate-treated rats than in Sham and higher Young's modulus in raloxifene-treated OVX rats compared to OVX-veh. Though we did not observe significant treatment effects on either perilacunar or bone matrix Young's modulus, we did observe a significant difference between the perilacunar and bone matrix Young's moduli. Our results support those of previous studies that showed a significant transition of mechanical properties from the perilacunar region to the surrounding bone matrix, with the bone matrix being stiffer than the perilacunar matrix (Meers et al., 2013; Meers, 2013), and previous findings showing that the perilacunar matrix is an active site of mineral exchange (Norimatsu et al., 1979; Matthews et al., 1978).

In a companion study, all these treatment sequences had higher maximum load, yield stress, and work to failure than the OVX group (Amugongo et al., 2014b). However, significant differences among the treatment groups for local mechanical properties measured by the Biodent system (Active Life Scientific; Santa Barbara, CA), did not exist (Amugongo et al., 2014b). Biodent measures of local mechanical properties can be considered to account for osteocyte lacunar density because Biodent's indentation tip (150 μm diameter) is 1000× larger than that of the Berkovich indenter used here, allowing multiple lacunae to be present under the Biodent tip. Perilacunar tissue is also grouped into the same Biodent material property measurement. In effect, lacunar density, lacunar size, and local material properties are combined in a single Biodent measurement. In comparison, the FE model in the present study takes into account differences in lacunar density, lacunar size, and local material properties. However, each variable is considered independently, and the resulting strain fringe plots in Fig. 4 can be calculated and considered. The microstructural and material changes at the lacunar level result in material properties that are similar at the 150 μm diameter interrogation level measured by Biodent. Differences in bone density, size, and/or geometry of the tibia likely contribute to the differences in whole bone strength at the right central tibia (Amugongo et al., 2014b). Though bone mass and geometry are surely determinants of whole bone strength, endpoints measured here are meaningful contributors to bone material properties that may reveal a more mechanistic view of how the skeleton adjusts its strength.

A weakness of our study is that we do not know if all the lacunae observed by nanoindentation contained viable osteocytes at necropsy; further studies that consider this are warranted. Some data suggest that aged bone does not initiate new bone formation at strain levels known to induce bone formation in younger bone (Rubin et al., 1992; Turner et al., 1995). It has been previously shown that raising or lowering the stiffness of the perilacunar region in this FE model amplifies or diminishes the strain transmitted to the osteocyte lacuna from a uniformly imposed global strain (Stern and Nicolella, 2013).

The use of linear elastic material models in the FE models was necessary as Young's modulus for the perilacunar and bone matrix regions for each of the models was determined through the nanoindentation study. The nanoindentation data was collected only in compression and along a single plane of the exposed bone surface. Using a viscoelastic anisotropic material model would be more representative of bone. However, multidimensional material properties were not determined for the bone samples.

The FE models employed an idealized ellipsoid shape for the lacuna. As seen in Fig. 2 and a number of recent publications (Carter et al., 2013; Mader et al., 2013; Hesse et al., 2014; Dong et al., 2014; Carter et al., 2014; Akhter et al., 2017), the lacuna is actually more of a flattened ellipsoid with a rough, irregular surface. In order for the model to remain parametric, canaliculi in the perilacunar and bone matrix were not added. FE models that incorporate an ability to vary lacunar shape in the third dimension have not yet been validated. The incorporation of the aforementioned features in the model would certainly affect the estimated strains, likely increasing the average strain values. However, these changes would likely have a similar effect across each of the models. The amount of strain actually transmitted to the embedded osteocyte is not known, but would also likely have similar effects across each of the models. Meanwhile, until such FE models become available, the state-of-the-art FE model applied here has been validated and thus supplies meaningful new information about the current topic (Ruffoni and van Lenthe, 2017).

Conclusion

In clinical practice today, many patients are cycled through a number of bone active medications for the treatment of osteoporosis over many years. We utilized detailed high resolution measurements of bone biomechanical properties in cortical bone of adult osteopenic, estrogen deficient rats, in combination with SEM and FE modeling to explore the relationship of structural differences in hard tissues to determine the effects on osteocyte strain. Using a cortical bone region in the rat, where bone remodeling is absent, we determined that Young's modulus was lower in the perilacunar region than in the bone matrix region, and that Young's modulus in the perilacunar region was negatively related to whole bone material properties. We found no significant treatment-related effects on Young's modulus in the perilacunar and bone matrix regions around osteocyte lacunae. Additional studies are now needed, particularly in trabecular bone regions in the rat where bone remodeling plays a prominent role (Erben, 1996; Wronski et al., 1999), to understand whether the same relationships hold in bone regions with a higher remodeling rate. In addition, data that would indicate morphologic and/or densitometric causes for the lower Young's modulus in the perilacunar region should be sought, perhaps with either non-destructive ultrahigh resolution (<50 nm) nanoCT imaging or destructive focused ion beam/scanning electron micropsy (FIB/SEM) (Taiwo et al., 2016; Schneider et al., 2010), both of which are capable of resolving individual canaliculi. Ultimately, the goal is to understand completely how these agents, individually, sequentially, or even in combination, alter bone quality in relation to fracture prevention, which is really the desired outcome of any treatments.

Disclosures

Authors report no conflicts of interest.

Authors' contributions

ARS co-drafted manuscript and executed FE models. XY, YW, AB, MD, and MLJ conducted the nanoindentation portion. WY directed the sequential treatment study. NEL and MLJ also co-drafted the manuscript. NEL is responsible for the final version.