Study on the biomechanical responses of the loaded bone in macroscale and mesoscale by multiscale poroelastic FE analysis

Background Bone is a hierarchically structured composite material, and different hierarchical levels exhibit diverse material properties and functions. The stress and strain distribution and fluid flow in bone play an important role in the realization of mechanotransduction and bone remodeling. Methods To investigate the mechanotransduction and fluid behaviors in loaded bone, a multiscale method was developed. Based on poroelastic theory, we established the theoretical and FE model of a segment bone to provide basis for researching more complex bone model. The COMSOL Multiphysics software was used to establish different scales of bone models, and the properties of mechanical and fluid behaviors in each scale were investigated. Results FE results correlated very well with analytical in macroscopic scale, and the results for the mesoscopic models were about less than 2% different compared to that in the macro–mesoscale models, verifying the correctness of the modeling. In macro–mesoscale, results demonstrated that variations in fluid pressure (FP), fluid velocity (FV), von Mises stress (VMS), and maximum principal strain (MPS) in the position of endosteum, periosteum, osteon, and interstitial bone and these variations can be considerable (up to 10, 8, 4 and 3.5 times difference in maximum FP, FV, VMS, and MPS between the highest and the lowest regions, respectively). With the changing of Young’s modulus (E) in each osteon lamella, the strain and stress concentration occurred in different positions and given rise to microscale spatial variations in the fluid pressure field. The heterogeneous distribution of lacunar–canalicular permeability (klcp) in each osteon lamella had various influence on the FP and FV, but had little effect on VMS and MPS. Conclusion Based on the idealized model presented in this article, the presence of endosteum and periosteum has an important influence on the fluid flow in bone. With the hypothetical parameter values in osteon lamellae, the bone material parameters have effect on the propagation of stress and fluid flow in bone. The model can also incorporate alternative material parameters obtained from different individuals. The suggested method is expected to provide dependable biological information for better understanding the bone mechanotransduction and signal transduction.

a detailed study on the propagation of physiological loading and fluid flow across the length scales of different functional units has not been made yet.
The objective of this research was to develop a multiscale model that included various functional units at each hierarchical level to evaluate the response of poroelastic behaviors under axial compressive cyclic loading. FE analysis is performed three times, at macroscale, macro-mesoscale, and mesoscale structural levels ( Fig. 1, Tables 1, 2, 3 and 4), with each analysis at a greater level of refinement, using COMSOL Multiphysics software. This model neglected the bone marrow cavity and trabeculae, only considers the tissue from the endosteum to periosteum, and the mechanical property and the flowing discipline of various functional units was observed. At the mesoscale, an osteon which cut from the whole model was refined, and the effects of E and k lcp on the fluid flow and stress and strain field were investigated. All properties and parameters were based on literature reports, and the solid structure and interstitial fluid were assumed as transverse isotropic poroelastic material and compressible liquid, respectively. This paper provides a deeper understanding of the mechanotransduction and stimulation by fluid flow which induces bone remodeling and bone metabolism. Macro-mesoscopic model including endosteum, periosteum, interstitial bone, and hollow osteon. c The hollow osteon in macro-mesoscopic model. d One-eight symmetry model. Regional M and line L are selected to analyze the biomechanical responses along the radius of bone tissue. e Region M is divided to 8 small regions (Regions [1][2][3][4][5][6][7][8], and each small region contains one osteon (Osteon 1-8). f Some points are selected at the internal wall, the middle and the external wall, and the interstitial bone. C Mesoscopic model. g FE mesoscopic model including macro-meso interface (including osteon lamellae structure). The critical part near the periosteum in macro-mesoscopic model is refinement. h A cube (340 µm) that cut from the macro-mesoscopic model near the periosteum as mesoscale model. Line H is a part of L along radius direction of done tissue. i Transversal cross-section of the mesoscale model with 10 osteon lamellae. j The cement line

Comparison of finite element method and numerical simulation
The comparison of FP and FV along the radial direction of macroscale model (Fig. 1A) was shown in Fig. 2, and the computed results of the macroscale FE model showed good agreement with the results of numerical simulation, and the result error was acceptable, which verified the validity of the model.

The analyses of macro-mesoscopic FE model
Because of the symmetry of geometric structure (Fig. 1B, C), in order to reduce the computation, 1/8 model was established (Fig. 1d). A region M was taken along the radius of the bone tissue (Fig. 1d), and some points were get at different positions of the endosteum, interstitial bone, osteon, and periosteum within the region M. The mean values of VMS ( Fig. 3A Fig. 4. The maximum VMS at interstitial bone was 2.24 × 10 7 Pa, which was    significantly higher than other locations (Fig. 4B). The maximum FP around the endosteum (14,786 Pa) and periosteum (5264 Pa) was significantly lower than other locations (Fig. 4C). In the same small region of bone tissue, the maximum FP of interstitial bone (5.9 × 10 4 Pa) was greater than osteon (4.8 × 10 4 Pa) (Fig. 4C), but the FV of interstitial tissue (1 × 10 −8 m/s) was less than osteon (8 × 10 −8 m/s) (Fig. 4D). The FV in the area where the periosteum and bone tissue contact had significant changes (Fig. 4D).

The analyses of osteon level in macro-mesoscopic FE model
Some points were selected at the internal wall, the middle, the external wall, and the interstitial bone in Region 1-Region 8 (Fig. 1e, f ) to analyze and compared the response of FP and FV, respectively. The peak values of Region 2-Region 7 have no significant difference; however, the peak values of Region 1 and Region 8 were significantly smaller than other regions. The peak FP decreases gradually from interstitial bone to the inner wall of osteon (Fig. 5A). The FV at the inner wall of osteon had the maximum value, and the range of peak FV in osteon was about 5 × 10 −8 m/s-8 × 10 −8 m/s (Fig. 5B). The value of FV in both interstitial bone and outer wall of osteon was less than 2 × 10 −8 m/s ( Fig. 5B), which meant that the fluid stimulation generated by fluid flow in these locations was too small to cause the response of osteocyte's mechanoreceptors and was called "dead zone" [14]. The FP (Fig. 5C) and FV ( Fig. 5D) of Osteon 2-7 showed the symmetrical distribution around the Haversian canal, but because of the influences of the boundary conditions, the FP and FV of Osteon 1 and Osteon 8 were asymmetric distribution. The FP and FV on the side near endosteum and periosteum both had a trend of decreasing markedly.

Effects of Young's modulus (E) and permeability (k lcp ) on mesoscale FE model
As a verification, the VMP and FP states at the cuts where the mesoscale models had displacements and pressure prescribed by the results in the whole model ( Fig. 1g) were checked. Such a comparison was shown in Fig. 6. For the range of parameters investigated, the computed peak VMS of the whole macro-mesoscopic models (Whole 1, 2, and 3) and mesoscale models (Case 1, 2, and 3) was about 2.6 × 10 7 Pa and 2.62 × 10 7 Pa, respectively, and the peak value of FP was about 4.8 × 10 4 Pa and 4.7 × 10 4 Pa, respectively. The peak values of VMS and FP for the whole models were about 0.1% ~ 1% and 1% ~ 2% different compared to that in the mesoscale models, which were expected given the coarse mesh in the whole model. The results strongly indicate that the mesoscale models had been set up correctly. The cut through the whole models (Whole 1, 2, and 3) in Fig. 6 displayed that the stresses and pressures were not well resolved. The stress and pressure fields are not smooth, especially around Haversian canal. Due to the large size span of the whole model, sometimes it is not feasible to have a mesh that at the same time captures the global behavior and resolves the detail structures with high accuracy. For thin bone lamella structures, it is common that the regions with different properties are small. The mesoscale models can simulate the detail structures more accurately. In  With the different E of each osteon lamella, the VMS of the bone lamella showed gradient change (Fig. 7A). As shown in Fig. 7B, the local MPS at medial lamellae position increased obviously in three cases, and the local MPS in Case 1 was significantly larger than Case 2 and Case 3. As shown in Fig. 7A, B, the peak VMS in Case 1 was at lateral lamellae, while the peak VMS occurred at medial lamellae in Case 2 and Case 3.

Discussion
This paper presented a multiscale poroelastic FE model under axial compressive cyclic loading to analyze the stress and strain field and fluid flow of bone. The distributions of the VMS, MPS, FP, and FV in endosteum, interstitial bone, osteon, and periosteum were analyzed on macro-mesoscopic model, and the effects of E and k lcp on stress and strain field and fluid flow in mesoscopic model were investigated.
The osteogenic function of periosteum is more important than the bone marrow and endosteum, and the periosteum has abundant pluripotent stem cells and molecular factors of regulating cell behavior play an important role in providing cells resource for bone repairing and bone healing [15,16]. The material properties of periosteum selected in our model, which was considered to have abundant osteoprogenitor cells   [17], made it more compliant. This had important effects as periosteum deform elastically at small strains. As shown in Fig. 4B, the MPS on the periosteum and endosteum was much larger than the interstitial bone and osteon, which results in a large strain gradient that may lead to greater strain stimulation on the bone cells on the endosteum and periosteum. However, the use of a transversely isotropic model for the periosteum and endosteum with a small elastic modulus of 4.41 MPa in radial direction and 25.67 MPa in axial direction resulted in a much lower VMS (Fig. 4A). The large VMS (Fig. 4A) in interstitial bone makes it easy to produce crack and fatigue damage, and the properties of high VMS and crack-prone of interstitial bone can provide a protection to the osteon, which helped to protect the osteocyte and maintained normal metabolic activity in the bone. This result was reasonable in the functions of maintaining structural integrity and resistance to fatigue damage in the process of biology evolution.
The maximum FP and FV of the whole bone tissue tend to be in a distribution of great in the middle and small at both ends (Fig. 5). The higher the FP and FV in osteon, the stronger the fluid stimulation is perceived by the osteocytes, which may lead to the physiological activities such as bone remodeling and bone metabolism more active. The distributions of FP and FV of Osteon 2-7 in Fig. 5 had little difference, which suggested that the fluid stimulation that the osteocytes felt in these osteons was basically the same. However, the asymmetrical distributions of FP and FV in Osteon 1 and Osteon 8 were found. The difference of FP and FV in Osteon 1, Osteon 8, and Osteon 2-7 may lead to different mechanosensation in osteocytes. The difference may relate to endosteum and periosteum. Studies had shown that the endosteum and the periosteum had important material exchange and signal transduction with bone tissue [18]. The periosteum and endosteum are closely related to the activities of bone formation and bone resorption, which together maintain the stability of the thickness of cortical bone [19]. The FV at the inner wall of the osteons was much larger than other positions, with a peak value about 8 × 10 −8 m/s (Fig. 5B). The maximum FV in both interstitial bone and outer wall of osteon was less than 2 × 10 −8 m/s (Fig. 5B), which means that the fluid stimulation generated by fluid flow in these locations was too small to cause the response of osteocyte's mechanoreceptors and was called "dead zone" [14].
The variations of E in 3 cases in the mesoscale model are all derived from the experimental observations [20][21][22][23]. As shown in Fig. 6, there were obvious high stress at the medial lamella in Case 2 and Case 3. High stress can lead to initiation of cracking by fatigue. Therefore, it is easy to appear microcrack damage near the Haversian canals in Case 2 and Case 3. Under physiological load, the bone formation of osteoblasts interacted with bone resorption of osteoclasts to complete the growth of new bone and the removing of old bone, and therefore, there would be more new bone formations in the large stress region to form strong bone lamellae. According to the geometry of osteon, the redistribution of stress in the existing structure will be hindered by Haversian canal, which should lead to the presence of stress concentration. The existence of stress concentration around the Haversian canal will lead to the increase of E (Case 2 and Case 3), which conformed to principle of bone adaptability. In Case 1 (Fig. 6), the lateral bone lamellae had higher stress and medial bone lamellae had lower stress, and the high stress region of lateral bone lamellae could cause stress shielding to medial bone lamellae. This stress shielding effect can prevent the stress concentration around Haversian canal to damage the osteon and provide protection to osteon [24]. Therefore, the results in Case 1 also conformed to the functions of maintaining structural integrity and resistance to fatigue damage in biological evolution. As shown in Fig. 7C, D, the different distributions of E had influence on the FP and FV, which indicated that the variations of E had effect on the fluid flow of osteon. The microcracks could affect the fluid flow behaviors of osteon [25]. The osteocytes around the Haversian canal could feel larger fluid stimulation (Figs. 5B, 7D and 8D), which may be related to the high stress around the Haversian canal and the susceptibility to initiation of microcracks. The existence of stress concentration and microcracks could make the bone formation more frequent to increase the strength of bone lamellae and repair microcracks.
Permeability can be regarded as a macro-index to describe the microscopic flow behavior in bone tissue [26,27]. According to previous theoretical and experimental results, the range of k lcp was 10 −17 m 2 to 10 −25 m 2 [28]. The variations of k lcp in bone lamellae might be related to the distribution of the bone canaliculi and lacuna, the low-permeability regions might have less bone canaliculi and lacuna [29]. Although researchers had done extensive research on the permeability of bone, the effect of gradient distribution of permeability of osteon lamellae on the biological responses of osteocyte under fluid flow was still unclear. Our study found that the heterogeneous distribution of permeability in osteon lamella had a significant effect on FP and FV (Fig. 8). The FV was found to fluctuate at bone lamellae and cement line that was different from the findings of [10,26]. In reality, the micropore of bone tissue may pass through the cement line [30], so the cement line may be permeable. Although the permeability value of the cement line is still unclear, this permeability of the cement line may change the mechanical behavior of cortical tissue. Therefore, we propose that the osteon in the ideal model is covered by a cement line. In addition, because the permeability of the cement line is lower than that of the osteon, the flow velocity at the cement line increases obviously, which will increase the fluid stimulation of the osteocytes near the outer wall of osteon. This mesoscale model was investigated by multiscale method, and the boundary conditions were more approach to real mechanical environment.
This paper established multiscale FE model to study strain and stress distributions and the fluid flow based on the poroelastic theory and provided a foundation for further exploration of the micromechanism of the growth and differentiation of osteocyte. Furthermore, we were restricted to simplified finite element models, whereas the bone tissue has more complex geometric structure, boundary conditions, and physiological load conditions. Incorporation of a more realistic model of bone tissue would prove useful for a more accurate representation of the different responses of stress and strain distribution and fluid flow when under physiological load. Future refinements of the model will address the inherent limitations. Nevertheless, this study represented a significant step in developing a multiscale model of bone tissue that incorporated explicit representation of osteon lamellae and cement line parameters. The future work is to establish microscale model based on the multiscale method and further quantify the load and fluid flow signal transfer from macroscale and mesoscale to microscale.

Conclusions
In summary, a multiscale poroelastic finite element method for the hierarchical structure of bone tissue was developed under axial compressive cyclic load, and the distribution of stress and strain and fluid flow in bone were investigated in different scales. Multiscale method can reflect the real physiological environment of different layers more accurately than research each component of bone separately. The further model can be used to analyze the effect of bone scaffold, bone substitute or implants on the stress and fluid flow distributions of bone tissue to more accurately assess the potential beneficial and harmful effects, so as to accurately achieve better individual matching. This work provides a better understanding of fluid flow and mechanotransduction in bone remodeling.

Governing equations for poroelastic bone model
Poroelasticity was first developed by Biot and widely used in solid-liquid coupling poroelastic materials [31]. The following governing equations can be used to describe the poroelastic behavior of bone, and no body forces are taken into account [14,26]. Constitutive laws for the solid matrix material and the saturating fluid can be written as: where σ is the total stress tensor, M is the stiffness tensor of the drained porous matrix, ε is the total strain tensor, α is the Biot's tensor, p is the pore pressure, M is the Biot's modulus, ξ is the variations in fluid content, and tr() is the trace operator. The equilibrium equation: here, the total density ρ is related to the porosity φ , the density of solid phase ρ s , and the density of fluid phase ρ f , and the ρ can be defined by relation ρ = φρ f + (1 − φ)ρ s , and ü s is the second derivative of the displacement. Fluid mass conservation equation: Darcy's law: here, V is the velocity vector and k is the intrinsic permeability tensor which can be defined by the relation K = k/u , where K is the permeability tensor and u is the dynamic viscosity of interstitial fluid. Inertia items in Eq. (3) and (5) can be ignored, because the bone is always subjected to low frequency cyclic loading in daily life and it has little effect on computing results. The simplified governing Eq. (6) is obtained by plugging (1) into (3) and plugging (2) and (5) into (4):

Establishment of macroscale poroelastic mathematical model and FE model
In macroscale model, we neglected the bone marrow cavity and trabecular bone and established a hollow bone model in the cylindrical coordinate system (r, φ, z). As shown in Fig. 1A, a, b, and h represent inner radius, outer radius, and height, respectively. The surface of periosteum was set no-flow conditions [32]. Equation (7), reported by literature [33,34], models a modified physiological arterial pressure pulse (in mmHg) as the hydrostatic pressure for the surface of endosteum.
The cyclic loads were axial to represent longitudinal compression on the top and bottom surface [32,35]. The solution of fluid pressure can be obtained: where I n and K n represent the first kind and the second kind modified Bessel function of order n (n = 0, 1), respectively. ε z0 is amplitude of axial strain, E is drained Young's modulus (Pa), and ν is Poisson's ratio of the drained porous matrix. α is the Biot-Willis coefficient, and c represents the integral constant determined by the boundary conditions. The E r , ν r and E z , ν z represent radial and axial drained Young's modulus and Poisson's ratio, respectively. The constant C can be determined by the relation: where i = √ −1 , ω represents load frequency, and k is intrinsic permeability.The fluid velocity V can be given by Darcy's Law: In order to provide the basis for researching the stress and strain field and fluid flow in various functional units by poroelastic FE method, the validity of poroelastic FE model was validated by comparing numerical result with simulation results. The analytical solution was obtained by MATLAB software. The poroelastic FE model was established by COMSOL Multiphysics. To simulate the mechanical environment of bone tissue, the load-displacement w of z direction was axial to represent longitudinal compression on the top and bottom surface. The amplitude of harmonic displacement (w) and frequency (f) is 0.5 μm and 1 Hz, respectively. The maximum axial strain ( ε ) is 1000 µε in a cycle: There are two kinds of porosities associated with bone fluid: the vascular porosity (order 20 μm) and the lacunar-canalicular porosity (order 0.1 μm). The range of permeability (k vp and k lcp ) that describes the fluid flow in the vascular porosity ( φ v ) and the lacunar-canalicular porosity ( φ lc ) are 10 −13 m 2 -10 −17 m 2 and 10 −17 m 2 -10 −25 m 2 , respectively [28]. We set k vp = 10 −15 m 2 and porosity φ v = 0.04 to describe fluid flow in the vascular porosity and k lcp = 10 −19 m 2 and porosity φ lc = 0.05 to describe fluid flow in the lacunar-canalicular porosity. Croker et al. found that the cross-sectional diameter and thickness of cortical bone of human femur (22.1-32.8 mm) were thicker than animals by analysis and comparison of human, sheep, and kangaroos using statistical methods [36]. In order to simplify calculation and facilitate comparison of the upcoming animal experiment, the cross-sectional diameter and the height of segment bone are set 1 cm and 1 mm, respectively. The related parameters are shown in Table 1 [14,27,37].

Establishment of macro-mesoscopic FE model
It is difficult to study the relationship of fluid flow within different levels and functional units by experimental and theoretical methods; therefore, the FE simulation method becomes the first choice. As shown in Fig. 1B, the macro-mesoscopic model considered the endosteum, osteon, interstitial bone, and periosteum, but the Haversian canal and the marrow cavity in bone were neglected (Fig. 1C). As shown in Fig. 1d, due to symmetry, only 1/8 model was observed in the computations. The radius of osteon was set R o = 150 μm, and the Haversian canal radius is set r o = 50 μm. The periosteum is a strong connective tissue envelope covering the bone surface except joints. It was tightly bound to the outer wall of bone tissue, and the thickness was set 150 μm [38]. The endosteum was a thin connective tissue envelope covering the bone marrow cavity and bone trabecula, and the thickness was set 50 μm.

Boundary conditions and material parameters
The surface of periosteum and bone marrow were set to the same value as the macroscopic model. The FP of Haversian canal was always ignored in previous study [10,14,27]; however, the span of both k vp and k lcp was very large [28], and when the gap of permeability value between k vp and k lcp was not that big, the FP of Haversian canal couldn't be ignored. In macro-mesoscopic model, the FP of Haversian canal was derived from the calculation results of the corresponding locations of macroscale model. The loaddisplacement is the same as macroscale model.
Due to the symmetry of full model, the both sides and bottom surface of the segment 1/8 model which cut from the whole model applied constrained symmetrically to prevent rigid body motion. Because the endosteum, osteon, and interstitial bone and periosteum were tightly bound together, all geometric objects were tied together by using the function of Form Union in COMSOL Multiphysics software to make the condition of the mesh border lines continuous. (11)  The Young's modulus (E) of interstitial bone was 10% larger than osteon, and the Poisson ratio (ν) was 10% smaller than osteon [39]. The transverse isotropic elastic constants for cortical bone were used. McBride et al. measured that the axial and radial elastic modulus of periosteum were 18.8-32.5 MPa and 3.19-5.62 MPa [40], respectively, and the average values 25.65 MPa and 4.41 MPa were used in this model. The Poisson's ratio was 0.49 [19]. The permeability of the periosteum from the bone to muscle surface of periosteum and the muscle to bone surface of the periosteum was much different, and we only considered the permeability (2.7 × 10 −16 m 2 ) between the periosteum and the bone tissue [41]. For other poroelastic material parameters, we set the same value on macroscopic and macro-mesoscopic model. The material parameters used in the FE model are shown in Table 2 [14,19,37,40,41].
Different precision of meshes were studied, to demonstrate that the results were converged with respect to mesh refinement-any further refinement of the mesh would only marginally improve the precision of the results. A Free Tetrahedral mesh was used and included 48,375 elements.

Establishment of FE mesoscale model including macro-meso interface
The establishment of mesoscale model consisted of two steps. Firstly, the macroscopic and macro-mesoscopic models were analyzed in order to capture the general trends and to identify the critical part. Secondly, a fine mesoscale model containing the critical part was made, and the study was resolved. In order to allow for smooth transition between macroscale and mesoscale, the FE mesoscopic model including macro-meso interface was established (Fig. 1g). The change of material parameters was consistent with the mesoscale model.

Establishment of mesoscopic FE model
The osteon was a basic unit for organizing and constructing the cortical bone, which composes of multilayered quasi-cylindrical composites of compact tissue arranged in 7-10-µm-thick lamella around the Haversian canal (Fig. 1i) [20]. Each osteon lamella had different properties [26]. Each osteon was encircled in a 1-5-μm-thick interface structure form, which was called the cement line, to separate osteon from interstitial bone (Fig. 1j). In the analysis, 10 bone lamellae were established. The cement line was set 1 µm [42]. The global effects from the macroscale model were transferred to the mesoscale model via appropriate boundary conditions. By the analysis of the simulation results of macroscale model, we confirmed the critical part near the outer region of periosteum as the position of mesoscale model. As shown in Fig. 1h, the mesoscale model was a cube (340 µm) that cut from the macro-mesoscopic model near the periosteum, which contained one osteon.

Boundary conditions and material parameters of the mesoscale model
The mesoscale model was based on the coupling of Structural Materials Module and Darcy's law in COMSOL Multiphysics. Since the model contained the properties of structural mechanics and fluid flow behaviors that were modeled with a poroelastic material, many degrees of freedom were required in order to obtain the accurate