Evaluating the Role of Canalicular Morphology and Perilacunar Region Properties on Local Mechanical Environment of Lacunar–Canalicular Network Using Finite Element Modeling

Osteocytes are mechanosensory cells which are located in cavities called lacunae throughout the bone matrix. Long channels emanating from lacunae, called canaliculi, house the osteocyte cell processes forming the lacunar–canalicular network (LCN), which enables fluid transport, nutrient exchange, and signaling between osteocytes.

Osteocytes can alter their local mechanical environment through perilacunar/canalicular remodeling that may influence both the morphology of LCN and the composition of the tissue around lacunae [1,2]. Several previous studies demonstrated the change in osteocyte lacunar morphology in human bone due to physiological and pathological processes. These studies demonstrated (1) reduction in osteocyte lacunar density with age, diseases, and in fracture patients compared to healthy ones [3–9]; (2) reduction in osteocyte lacunae volume with age and disease [5,10]; and (3) increase in sphericity of lacunae with age [5]. Although not in human bone, osteocyte lacunae size was shown to increase with lactation, and decrease back to control levels after weaning in mice and rat bone [11–14]. There are limited data on the modifications in canalicular morphology due to physiological or pathological processes both in human and animal bone. Canalicular density was found to be lower in aged human bone compared to that of young [15]. Several other studies on mice and rats reported a reduction in canalicular density and/or length with aging, diabetes, perlecan deficiency, glucocorticoid-induced osteoporosis, defective transforming growth factor beta (TGF-β) signaling, and tumor growth in cancer metastases [16–22]. In addition, canalicular area increased during lactation and reduced with aging and diabetes in mice bone [12,21]. Impaired canalicular connectivity was also observed with aging, TGF-β signaling, and with X-linked hypophosphatemia in mice bone [19,23].

Besides morphological changes in LCN, modifications in perilacunar elastic modulus compared to the bulk modulus were also observed in mice and rat bones. A recent study demonstrated a gradation in perilacunar elastic modulus around lacunae in mice bone [24]. In this study, the perilacunar elastic modulus was lower than the bulk matrix modulus adjacent to the lacunae (∼0.2 μm from the lacunar wall), increased to a peak value (∼0.4 μm from the lacunar wall), and then decreased to the bulk matrix modulus value (∼2 μm from the lacunar wall) [24]. More pronounced changes were observed in perilacunar elastic modulus of mice bone that underwent glucocorticoid treatment [25] and with chronic kidney disease [26]. Modest differences between perilacunar modulus and bulk modulus were also reported in control, ovariectomized, and osteoporosis-treated mice bone [27]. Taken together, morphological changes in the LCN as well as modifications in the perilacunar region properties can impact the mechanical environment of osteocytes. Therefore, understanding the role of LCN morphology and perilacunar region properties on the local mechanical environment of osteocytes can provide additional insight into how osteocyte mechanosensitivity is modified with physiological or pathological processes.

External forces applied to the skeletal system deform the bone tissue around LCN inducing fluid flow through the LCN and applying shear stress on the osteocytes and their processes. To gain more insight into the mechanical environment around LCN local tissue strains were measured with advanced imaging technologies. These measurements revealed average strain amplification factors of 1.1 to 3.8 around osteocyte lacunae compared to macroscopic strains as well as local peak strain levels of 30,000 με under an applied load of 2000 με [28]. Another study reported a range of maximum strains in osteocytes around 19,000 to 30,000 in rat bone under a compressive loading of 3000 με [29]. These results indicate that in vivo strains (1000–3000 με) during moderate to vigorous physiological activities [30,31] due to externally applied loads can be magnified significantly around the lacunae.

Finite element (FE) models have also been utilized to overcome the experimental challenges of measuring the local deformation around LCN. Strains around single osteocyte lacuna FE models containing perilacunar region were amplified by 2.5 to 3 with canaliculi and 1.3 to 1.5 times without canaliculi compared to macroscopic strain levels [27,32]. In FE models of a single osteocyte lacuna including an osteocyte, osteocyte processes, extracellular matrix (ECM), and pericellular matrix (PCM) showed strain amplification of 3–4 times in idealized models [33–35] and 8–70 times in image-derived models [36,37]. The lacunar orientation and size as well as its distance to a Volkmann canal were also shown to influence the strain around lacunae in image-based and idealized models of multiple osteocyte lacunae [38–40]. Our previous FE study utilizing three-dimensional (3D) models of multiple osteocyte lacunae examined the influence of osteocyte lacunar morphology and perilacunar region properties on the mechanical environment around osteocyte lacunae [41]. The results showed that osteocyte lacunar orientation, volume, and density as well as perilacunar elastic modulus altered the extent of highly strained tissue around lacunae influencing the mechanical environment of osteocyte lacunae.

In the current study, we aim to extend our previous study and evaluate the influence of canalicular morphology in combination with lacunar density, and the perilacunar region properties on the local mechanical environment around LCN. To address this aim, we utilized 3D FE models of a representative volume of bone with several lacunae including their canaliculi and varied canalicular morphology, lacunar density, and perilacunar region properties in the models based on previously reported experimental measurements. Our models differ from previous models in the literature as they incorporate multiple lacunae and their canaliculi explicitly. As a result, they capture the interaction of multiple lacunae with canaliculi and how they influence the local mechanical environment of the LCN. In addition, this study presents a thorough and systematic analysis of how variation in critical canalicular morphological characteristics and perilacunar region properties that were observed to vary with physiological and pathological processes influence the mechanical environment around LCN. This information has not been presented as a whole in previous studies. These results can enhance the understanding of how changes in canalicular morphology and perilacunar region properties impact the mechanical environment around LCN.

The lacunar–canalicular morphology and perilacunar region properties that were investigated in the current study included (i) canalicular density (Ca.D), (ii) canalicular length (Ca.L), (iii) canalicular diameter (Ca.Dm), (iv) lacunar density (Lc.D), (v) perilacunar region size (Pl.S), and (vi) perilacunar elastic modulus (Pl.E). This resulted in seven model groups (Table 1). In each group, one of the properties was varied based on previous experimental measurements while the others remained constant (Table 1). Due to the limitations in generation of FE meshes with large number of canaliculi, lacunar density of group 1–3 was set at the lowest limit of 10,000 #/mm³ reported in the literature [43]. To facilitate the model generation, a slightly smaller canalicular length and a larger canalicular diameter was used in group 4 models compared to group 1–3. The models that varied lacunar density contained a group of lacunae that were present in all models to track the change in the mechanical environment around LCN in a consistent manner as the lacunar density increased. These lacunae are referred to as “common lacunae.” Additional simulations were performed introducing perilacunar regions in the models. The first two groups of these simulations (group 5 and 6) were derived from three models in group 2 with canalicular length of 3, 6, 9 μm, respectively, and the third group (group 7) was derived from all models in group 4.

The surrounding bone matrix around the LCN was represented as a cubic volume with a size of 5 × 10⁻⁴ mm³ (0.079 × 0.079 × 0.079 mm³) for groups 1–3, 5, 6 (five lacunae/model) and 1 × 10⁻³ mm³ (0.1 × 0.1 × 0.1 mm³) in group 4 and 7 (10–25 lacunae/model). In the 3D models, lacunae were represented as ellipsoidal voids and canaliculi as cylindrical tunnels emanating normal to the lacunar surface (Fig. 1(a)). Individual osteocyte lacuna within the 3D models with different canalicular morphology and perilacunar region are shown in Figs. 1(b)–1(e).

The coordinates of lacunar centroids were generated randomly via an in-house developed MATLAB (version R2019a, MathWorks, Natick, MA) script from our previous work [41]. The lacunae size (460 μm³) and orientation (longest axis of lacuna at 30 deg from loading direction) were held constant in all models. These values were chosen based on experimental data reported in the literature [36]. Three ellipsoidal axes of the lacunae were chosen based on previous experimental measurements as outlined in our previous study [41]. Canaliculi were generated as straight hollow cylinders emanating outward from lacunar surfaces. Coordinates of the points where canaliculi emanate from the lacunae were selected such that the distance between two canaliculi was larger than $Distmin=1/Ca.D·π$⁠. This value was calculated by first dividing the lacunar surface area with the number of canaliculi to identify a circular area associated with each canaliculus. The radius of this circular area was defined as the minimum distance between two canaliculi. This approach allows overlap between circular areas associated with each canaliculus. However, it ensures that there is only one canaliculus in each circular area leading to an even distribution of canaliculi over the lacunar surface. This morphological information was then translated into a SOLIDWORKS (version 2019, Dassault Systèmes SolidWorks Corporation) macrofile for generating lacunar assembly, and a Python script for generating canaliculi assembly in abaqus (version 6.18, Simulia, Providence, RI). All parts were then assembled in abaqus through Boolean operations, which cut out the lacunar and canalicular parts from the cubic volume (Fig. 1(a)). Perilacunar regions were generated by scaling up the original lacunae in all dimensions (Fig. 1(b)). Any features protruding out of the boundaries of the cube were truncated.

The elastic modulus and the Poisson's ratio of the bone matrix were set as 20 GPa and 0.3, respectively, based on previous studies [45,46]. Symmetric boundary conditions were implemented in the models that show similar results to periodic boundary conditions [47] (Fig. 2(a)). A tensile displacement loading was applied to the top face (Fig. 2(a)), resulting in a total of 3000 microstrains representing bone matrix strain during vigorous activities [31]. Tensile loading was chosen in this study to investigate the mechanical environment around LCN in the sections of long bones experiencing tension particularly due to bending such as in femur and tibia.

Models were meshed using four-node tetrahedral elements with a global element size of 4.0 μm and a minimum element size of 0.11 μm (Figs. 2(b) and 2(c)). The minimum element size was chosen to be smaller than the smallest feature in the models (smallest canaliculus diameter = 0.3 μm) to minimize the dependency of stress and strain on the mesh density [48,49]. A curvature control factor of 0.09 was applied to ensure at least eight divisions across canalicular radial sections (Fig. 2(c)). The element sizes (global element size = 4 μm and minimum element size = 0.11 μm) were selected through a mesh sensitivity study. This study showed a small variation in the apparent elastic modulus (±0.02%) with the global and minimum element size (Fig. S1(a) available in the Supplemental Materials on the ASME Digital Collection). The maximum principal strain also demonstrated a small variation as the minimum mesh size was varied (±3.4%) and this difference was 0.13% between the two smallest minimum element sizes (0.10 and 0.11 μm) (Fig. S1(b) available in the Supplemental Materials on the ASME Digital Collection). The selected mesh size resulted in a total number of 3 to 18 × 10⁶ elements depending on the canalicular morphology. The simulation time for each model ranged from 47 to 180 min based on the element number using 64 cores in a parallel computing environment.

A total of 28 simulations were performed using models with parameters described in Table 1 (five models for group 1 and 2, varying canalicular density and length, respectively, four models for group 3, 4, 7, varying canalicular diameter and lacunar density, respectively, and three models for group 5 and 6, varying perilacunar size and modulus, respectively). The apparent elastic modulus of the models was calculated by using the average stress on the loading surface (total load/top surface area) and the average strain (total displacement/original height of the model). In addition, the strain amplification factor (Lc.SAF) was calculated as the ratio of the maximum principal strain in the model, which occurs on the surface of LCN, to the total applied strain (3000 με). Furthermore, in order to compare the extent of highly strained regions around LCN, the volume of elements (E.Vol) that have strain over 1.5 times of the applied strain (4500 με) were extracted from each model. This value was selected based on the strain distribution in the models (details in Sec. 3.2). E.Vol was normalized with the total element volume in the entire model (M.E.Vol) and with the element volume associated only with LCN surface (Lc.E.Vol). M.E.Vol represents the deformation around LCN in the entire model while Lc.E.Vol corresponds to the deformation of LCN walls.

For all models, the change in apparent elastic modulus of the bone with canalicular morphology and lacunar density was within 1% due to the small change in total porosity (Fig. 3(a)). The apparent elastic modulus of the bone showed strong linear dependence on canalicular length (R² = 0.99), diameter (R² = 0.98), and lacunar density (R² = 0.99) but not on canalicular density (R² = 0.45) (Fig. S2 available in the Supplemental Materials on the ASME Digital Collection). There was a direct correlation between apparent elastic modulus and perilacunar elastic modulus (Fig. 3(b)) whereas perilacunar region size had a very small impact on the apparent modulus (Fig. 3(c)).

The strain concentration on the lacunar surface appeared as an annular region along the longest axis direction in all the models and showed no correlation with canalicular density (Fig. 4). Higher strain values were observed near the lacuna-canaliculi junctions within high-strain areas on the lacunar surface (Fig. 4). The strain levels at the junctions were around two times higher compared to models with no canaliculi (Fig. S3 available in the Supplemental Materials).

In the models, a very small number of elements had strains below 1000 με, while the majority of the elements had strains within 1000–4500 με (Fig. 5(a)). Elements strained over 4500 με were directly associated with LCN (Fig. 5(b)). Thus, the strain level of 4500 με was selected as the threshold for defining the volume of highly strained elements. This value was 1.5 times the applied strain to the model.

Normalized volume of highly strained elements in the model (M.E.Vol) increased linearly with increasing values of canalicular morphological properties (R² > 0.99 for all fits) (Fig. 6). Normalized volume of highly strained elements on LCN surfaces (Lc.E.Vol) increased linearly with increasing canalicular density (R² = 0.97, Fig. 6(a)), but showed a quadratic relationship with canalicular length (R² = 0.93, Fig. 6(b)) and an exponential relationship with canalicular diameter (R² > 0.99, Fig. 6(c)). Lacunar strain amplification factor varied between 2.3 and 5 in all models but did not show a clear relationship with any of the canalicular parameters (Fig. S4 available in the Supplemental Materials on the ASME Digital Collection).

Comparison of the common lacunae, (Cm. Lc, lacunae that were present in all models in group 4 as described in Sec. 2.1) in models with varying lacunar density showed that the normalized volume of highly strained elements on LCN surfaces associated with common lacunae (Cm.Lc.E.Vol) only varied slightly with increasing lacunar density (Fig. 7(a)). Normalized volume of highly strained elements in the model (M.E.Vol) increased with an increasing lacunar density whereas normalized volume of highly strained elements on LCN surfaces (Lc.E.Vol) did not show a dependence on lacunar density (Fig. 7(b)). In addition, the lacunar strain amplification factor in common lacunae (Cm.Lc.SAF) did not vary significantly with lacunar density (Figs. S5 and S6 available in the Supplemental Materials).

Normalized volume of highly strained elements in the model (M.E.Vol) increased by around one fold compared to models without perilacunar region when the perilacunar elastic modulus was reduced to 75% of the matrix elastic modulus (Fig. 8(a)). This increase in the normalized volume of highly strained elements in the model (M.E.Vol) was independent of lacunar density (Fig. S7 available in the Supplemental Materials on the ASME Digital Collection). The reduction in perilacunar elastic modulus to 75% of the matrix elastic modulus also led to a ∼50% increase in normalized volume of highly strained elements on LCN surfaces (Lc.E.Vol) (Fig. 8(b)).

Lacunar strain amplification factor decreased from a mean value of 3.9 to 3.2 as the perilacunar elastic modulus increased from 10 GPa to 20 GPa (Fig. 9(a)). Normalized volume of highly strained elements in the model (M.E.Vol) and on LCN surfaces (Lc.E.Vol) exhibited descending exponential (R² > 0.99) and linear (R² > 0.99) trends with increasing perilacunar elastic modulus, respectively (Fig. 9(b)). Lacunar strain amplification factor did not show a trend with increasing perilacunar region size (Fig. S8(a) available in the Supplemental Materials). Normalized volume of highly strained elements in the model (M.E.Vol) and on LCN surfaces (Lc.E.Vol) varied only slightly with an increase in perilacunar region size (Fig. S8(b) available in the Supplemental Materials).

This study utilized a 3D FE modeling framework that incorporated modifications in canalicular morphology, lacunar density, and perilacunar region properties based on experimental measurements in the literature. The modeling approach has the capability to quantify the changes in local mechanical environment of LCN as a result of the local morphological and tissue property changes around LCN.

The primary mechanism of mechanosensation for the osteocytes is fluid flow in the LCN [35,50] and the fluid flow associated drag/shear force on the osteocytes and their processes [51]. In this study, we focused on evaluating the local strain environment around LCN, particularly the extent of highly strained tissue, because fluid flow is altered by tissue deformation. The simulation outcomes showed that the extent of highly deformed regions around lacunae is more sensitive to changes in canalicular morphology and perilacunar region properties than the strain amplification factor. Therefore, the variation in the extent of highly strained tissue may have implications for fluid flow in LCN and in turn on osteocyte mechanosensation.

The bone apparent elastic modulus decreased with modifications in LCN morphology that increased porosity. This reduction was small because of the small increase in lacunar–canalicular porosity in the models. If a more substantial change in porosity occurs due to modifications in LCN morphology, a considerable change in the bone apparent elastic properties can occur such as seen in experimental measurements [12]. In addition, our simulation results suggest that changes in the elastic modulus of the perilacunar region may lead to a detectable change in the apparent elastic properties of the bone matrix. This indicates that the local tissue property changes around LCN may have impact on the mechanical properties at larger length scales.

The simulation results showed that the increase in canalicular density, length, and diameter did not have a large influence on the strain amplification factor but they increased the amount of highly strained bone matrix. In addition, increasing canalicular density and diameter increased the amount of highly strained elements that are associated with the lacunar–canalicular surfaces. Normalized volume of highly strained elements on LCN surfaces demonstrated a nonlinear increase with canalicular diameter indicating that a slight increase in canalicular diameter when the canalicular diameter is small may lead to a steep increase in the amount of highly deformed tissue on LCN surfaces. However, this effect diminished as the canalicular diameter became larger. In addition, a quadratic relationship was observed between the amount of normalized volume of highly strained elements on LCN surfaces and the canalicular length when averaged over the entire model. However, when the trends for individual lacunae were evaluated, some of the lacunae demonstrated a steady linear decrease with increasing canalicular length whereas others showed an initial increase followed by a decrease (Fig. S9 available in the Supplemental Materials). With the current data, it is difficult to conclude whether the highly deformed element volume near lacunar–canalicular surfaces reach a peak value and decrease. Further work with more lacunae included in each model is required to answer this question.

These results indicate that changes in the canalicular morphology influence the amount of bone tissue that is highly strained and in turn may impact deformation driven interstitial fluid flow in the LCN. Based on our results, the experimentally observed reduction in the canalicular density and/or length with aging, diabetes, perlecan deficiency, glucocorticoid-induced osteoporosis, defective TGF-β signaling, and tumor growth in cancer metastases [15–22] may reduce the amount of highly strained bone matrix and potentially reduce the fluid exchange rate in the LCN. This may adversely impact the local mechanical environment of osteocytes. Previous studies showed that the change in the canalicular morphology can impact fluid flow characteristics in LCN. In particular, reduction in canalicular density and increase in canalicular length decreased the solute transport rate and perilacunar diffusion [19,52,53]. Our results complement these results and indicate that the degradation in canalicular morphology, e.g., canalicular density and diameter, reduce the amount of highly strained bone matrix around LCN that can potentially deteriorate fluid flow.

In addition, the lacunar-canaliculi junctions in our models had higher strains compared to the lacunar surfaces [33,36]. This indicates that canaliculi may further amplify the strain transferred to the osteocyte dendrites that are connected to the canalicular wall. This can have an impact on the mechanosensation of the osteocytes because of the critical role of the osteocyte dendrites in mechanotransduction [54].

The increase in lacunar density did not have an influence on the strain amplification (Lc.SAF) and highly strained regions on lacunae and canaliculi surfaces (Lc.E.Vol). However, it increased the total highly strained element volume (M.E.Vol) because of an increase in number of lacunae in the models. This suggests that a higher Lc.D can significantly increase the amount of bone tissue under high deformation, but has minor influence on the extent of deformation on individual lacunar–canalicular surfaces. This may indicate that the adverse impact of reduction in lacunar density, as observed with aging and osteoporosis [43,55], may be due to its negative influence on the connectivity of the osteocytes and in turn on fluid flow, signaling, and nutrient exchange but not due to its impact on local deformation of existing lacunae.

Our simulation results showed that a reduction in perilacunar elastic modulus increased both the strain amplification and the volume of highly strained tissue around and on the surface of LCN (Fig. 9). The change in perilacunar elastic modulus had an impact on apparent and local mechanical properties even with a relatively small perilacunar region size, which is consistent with the outcomes of our previous study [41]. These findings indicate that reduced perilacunar elastic modulus may enhance fluid flow as a result of the increased strain amplification and volume of highly strained tissue on and around LCN and lead to elevated strain sensitivity and mechanotransduction for osteocytes. Reduced perilacunar elastic modulus was observed in young and old mice bone [24], in mice bone under various osteoporosis treatments [27], and in more pronounced amounts in mice bone under glucocorticoid treatment, and with chronic kidney disease [25,26]. Our simulation results combined with experimental observations may indicate that changes in the perilacunar modulus may be a mechanisms to restore osteocyte mechanosensitivity to compensate for adverse effects of physiological and pathological processes in bone. Conversely, an increase in perilacunar elastic modulus may lower the strains and the extent of highly deformed tissue around LCN and reduce fluid flow mediated mechanosensitivity of osteocytes adversely impacting bone quality. This can occur due to higher mineralization around the lacunae such as in the elevated mineralization around the edge of the lacunae in estrogen deficient (ovariectomized) mice bone [56].

The main strength of our FE models is that they included the explicit representation of multiple lacunae including their canaliculi which is different from existing FE models in the literature that mainly focused on a single lacuna and its associated canaliculi. Limited number of existing multiple lacunae models either did not include canaliculi or vary canalicular morphology. By incorporating multiple lacunae and their canaliculi explicitly in our models, we were able to capture the interaction of multiple lacunae/canaliculi and how they influence the local mechanical environment of the LCN. In addition, this study presented a thorough and systematic analysis of how variation in critical canalicular morphological characteristics and perilacunar region properties that were observed to vary with physiological and pathological processes influence the mechanical environment around LCN. This information has not been presented as a whole in previous studies. In addition, the developed model makes it possible to evaluate the influence of individual variations in the lacunar–canalicular morphology and perilacunar region properties on the local mechanical environment of lacunar–canalicular network. In real specimens, concurrent changes in these parameters may occur and may confound the individual influence of each parameter. The simulation results provide additional insights into the relationship between canalicular morphology, perilacunar region properties, and the local mechanical environment of lacunar–canalicular network which has not been determined experimentally.

The main limitation in our models is the idealized geometries of the lacunae and canaliculi which may underestimate the strain amplifications that can occur due to the actual irregular geometry of these features. The strain amplification observed around LCN in our models were similar to other idealized and microcomputed tomography derived models [32,34,35,40] as well as experimental measurements [28]. However, our values were smaller than the ones observed in image-derived FE models that incorporate osteocytes and the surrounding PCM and ECM [36,37]. Previous computational studies showed that ECM strains were magnified when they were transferred to the PCM and the osteocytes [33,36,40]. Our models also did not include PCM, cells, and the extracellular fluid in the LCN and did not directly measure the strains that are transferred from the bone matrix to these features. The reason for this choice was to isolate the relationship between canalicular morphology and the tissue deformation around LCN. Despite the potential underestimation of the strain amplification factors around LCN, the extent of highly deformed tissue, which was the primary outcome of our simulations, is not expected to be influenced by the idealized LCN in the models. In our models, the canaliculi network was not represented as a connected network with branching and tortuosity due to the complexity of model generation. Steep variations in strain can be observed due to canalicular tortuosity. The representation of branching and tortuosity in the models may further elevate the impact of canalicular morphology on the extent of highly deformed tissue around LCN compared to values found in the current models. In addition, the variation in lacunar orientation and size can impact the strain amplification in the models as we have demonstrated in our previous study that evaluated the influence of lacunar morphology on mechanical environment of lacunae independent of canalicular morphology [41]. In the current models, these values were held constant to isolate the influence of canalicular morphology on local mechanical environment of LCN. Another limitation is that the models assumed a uniform low or high perilacunar elastic modulus compared to the bulk modulus whereas perilacunar modulus has been shown to be graded as a function of distance from the lacunar wall [24]. Challenges remain for generating meshes with graded material properties in models incorporating canaliculi because of additional partitions needed which will lead to many small intersecting surfaces complicating the meshing process. Similarly, elastic modulus may be graded not only around lacunae but adjacent to the canaliculi as well based on mass density and mineralization measurements in the literature [57,58]. Incorporation of such material gradation will also lead to challenges in meshing due to multiple layers of partitions that will be required. A specific constitutive material model that incorporates graded material property variation will need to be developed to overcome these meshing challenges. The model developed in this study was not directly validated because experimental data that correspond to the model outcomes are not available due to the difficulty of carrying out mechanical tests at this small length scale. However, the simulation results are consistent with experimental results as outlined above and support the modeling approach developed in this study.

In summary, this study presented a 3D FE modeling framework for evaluating the influence of modifications in LCN morphology and perilacunar region properties on the local mechanical environment of LCN as well as the apparent mechanical properties of the bone. The results provided additional insights on how the local mechanical environment around LCN are modified by changes in the canalicular morphology, lacunar density as well as perilacunar region properties. Specifically, the variation in canalicular morphology did not influence the level of strain amplification. However, increased canalicular density and diameter resulted in larger highly strained regions around and on the surface of LCN. Longer canaliculi increased the highly strained volume around LCN but reduced the highly strained regions on the surface of lacunae and canaliculi. A decrease in perilacunar modulus increased the lacunar strain amplification and the extent of highly strained regions around LCN and on LCN surfaces while perilacunar size did not influence these outcomes. The FE model presented here can be used to further understand how the mechanosensitivity of osteocytes may be impacted from the changes in LCN morphology and perilacunar region properties due to physiological and pathological processes such as aging, disease, lactation, and treatments.

This work used the Extreme Science and Engineering Discovery Environment (XSEDE) which is supported by National Science Foundation Grant No. ACI-1548562 through the allocation MSS-160001.

• National Science Foundation (Grant No. ACI-1548562 through the allocation MSS-160001; Funder ID: 10.13039/100000001).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.