A multiscale 3D f i nite element analysis of f l uid/solute transport in mechanically loaded bone

Lixia Fan*,Shaopeng Pei*,X Lucas Luand Liyun Wang

A multiscale 3D f i nite element analysis of f l uid/solute transport in mechanically loaded bone

Lixia Fan1,2*,Shaopeng Pei1*,X Lucas Lu1and Liyun Wang1

The transport of f l uid,nutrients,and signaling molecules in the bone lacunar–canalicular system(LCS)is critical for osteocyte survival and function.We have applied the f l uorescence recovery after photobleaching (FRAP)approach to quantify load-induced f l uid and solute transport in the LCS in situ,but the measurements were limited to cortical regions 30–50 μm underneath the periosteum due to the constrains of laser penetration.With this work,we aimed to expand our understanding of load-induced f l uid and solute transport in both trabecular and cortical bone using a multiscaled image-based f i nite element analysis(FEA) approach.An intact murine tibia was f i rst re-constructed from microCT images into a three-dimensional(3D) linear elastic FEA model,and the matrix deformations at various locations were calculated under axial loading.A segment of the above 3D model was then imported to the biphasic poroelasticity analysis platform (FEBio)to predict load-induced f l uid pressure f i elds,and interstitial solute/f l uid f l ows through LCS in both cortical and trabecular regions.Further,secondary f l ow effects such as the shear stress and/or drag force acting on osteocytes,the presumed mechano-sensors in bone,were derived using the previously developed ultrastructural model of Brinkman f l ow in the canaliculi.The material properties assumed in the FEA models were validated against previously obtained strain and FRAP transport data measured on the cortical cortex. Our results demonstrated the feasibility of this computational approach in estimating the f l uid f l ux in the LCS and the cellular stimulation forces(shear and drag forces)for osteocytes in any cortical and trabecular bone locations,allowing further studies of how the activation of osteocytes correlates with in vivo functional bone formation.The study provides a promising platform to reveal potential cellular mechanisms underlying the anabolic power of exercises and physical activities in treating patients with skeletal de f i ciencies.

INTRODUCTION

It is well known that bone tissue is capable of adapting its mass and structure in response to mechanical cues.1Although the cellular and molecular mechanisms of how the mechanical environment affects bone tissue are still not well understood,various mechanical parameters including(but not limited to)matrix deformations(strains), interstitial f l uid pressure,and f l uid f l ow-induced shear and drag forces have been found to impact bone’s responses to mechanical loading at cellular and tissue levels.2Quanti f i cation of the strains associated with physiological mechanical stimuli in bone has been performed at both tissue and cellular levels.3–5Using strain gages,the tissuelevel strains were found to vary from～600 μ? during the light activity of walking to～2 000 μ? during vigorous activities such as running and jumping.3–4Strains at the cellular levels have been mapped recently using fi nite element analysis(FEA)or imaging correlation techniques, and inhomogeneous strains(0.8%–3%)were recorded near the lacunar pores.5Whether these matrix strains directly excite osteocytes,the presumed mechanical sensors in bone remain debatable.Previous experiments have shown that bone cells were more sensitive to loading-induced fl uid fl ow than matrix strains.6–7Recent studies stronglysuggest that osteocytes,due to their large number and strategic positioning in the bone matrix,have very important roles in bone adaptation and metabolism.2,8–10These multi-functioning cells form an extensive and wellconnected network that optimizes them for detecting external mechanical stimuli,for example,f l uid f l ow in the canaliculi driven by load-induced matrix deformations.In response,osteocytes release various soluble bioactive factors,which modulate the functions of other bone cells and trigger biological processes such as osteoclastic resorption during overuse and disuse as well as load-induced osteoblastic bone formation.11–13The microscopic lacunar–canalicular system(LCS)that houses the osteocytes within the mineralized matrix is one key feature that enables the osteocytes to perform these important functions.2,9–10Not only does the LCS provide the critical“l(fā)ife line”for nutrient supply and“network”for cell signaling14–16but it also possesses the structural components(for example,tethering f i bers within the f l uid annulus)to amplify the loading signals and convert the overall loading to cellular stimulation forces such as shear stress and/or drag force acting on the osteocytes.17–18

We and others have built mathematical models to predict the magnitudes of load-induced f l uid f l ow and solute transport in the LCS19–26based on the groundbreaking paper by Weinbaum et al.8A mathematical model was developed to investigate the net solute(for example,nutrients and tracers)transport in the discrete LCS channels during cyclic loading of bone,and solute mixing within the extracellular space in lacunae was found to be responsible for the net solute transport.25To guide the experimental investigation of solute convections using the f l uorescence recovery after photobleaching(FRAP) approach that was developed initially for diffusion studies,16Zhou et al.26developed a one-dimensional, three-compartment model to simulate load-induced solute transport in the LCS during FRAP experiments. This modeling work enabled the use of experimental FRAP measurements at the lacunar level to predict fl uid fl ow in the canalicular level in a later study.27We further expanded the three-compartment model by considering the solute–matrix interaction as exogenous tracer probes moved through the osteocytic pericellular matrix,28allowing us for the fi rst time to quantify the average fi ber spacing in bones from young and aged mice with normal or altered pericellular matrix conditions.29However, due to limited laser penetration into mineralized bone tissue,FRAP measurements have been limited only to the cortical bone 30–50μm underneath the periosteum of the tibial cortex.16,27–29Currently,transport measurements on deeper cortex and trabecular bone areas are lacking. These locations are of biological signi fi cance because they typically undergo changes under mechanical loading or disuse conditions.30–31

FEA is a powerful numerical tool of analyzing stress–strain fi elds in objects of irregular shapes and has been extensively used in bone research.A multilevel FEA model of a human femur cortex was developed to predict local matrix deformations at the osteocyte level during normal gait.32To predict interstitial fl uid fl ow in bone,a two-dimensional FEA poroelasticity model was developed to investigate the load-induced fl uid fl ow in cortical bone.21In another study,solute/ fl uid fl ow fi elds were predicted in an FEA model of a rat tibia under four-point bending.33All of these models have focused on cortical bone sites.In this study,we aimed to develop a multiscale approach combining FEA and ultrastructural modeling. Our objective was to test the feasibility of the multiscale approach to predict functional outcomes resulted from mechanical loading.These outcomes include macroscopic(whole-bone level)and microscopic(LCS level) distribution fi elds of matrix strains,interstitial fl uid pressure,as well as stimulation( fl uid shearing and drag forces)at the cellular level.The rich pool of experimental data that we have obtained using FRAP on cortical bone16,27–28,34was used to validate the material properties assumed in the multiscale model.This greatly improved the fi delity of our model predictions on deeper cortex and trabecular bone regions,which are inaccessible with current FRAP techniques.

MATERIALS AND METHODS

Whole-bone model

An intact mouse tibia from an adult C57BL/6J male mouse was imaged by a Scanco μCT35 scanner(Scanco USA, Inc.,Wayne,PA,USA)using a standard protocol(55 keV, 145 μA,200 ms integration time,3 600 projections,and 20 μm voxel size).The raw image slices(998 slices)were imported in the DICOM format into ScanIP(Simpleware, Chantilly,VA,USA),with which the entire tibia,including the cortical and trabecular bone,was thresholded and meshed with 5 112 690 tetra elements(Figure 1a).In Hypermesh(Altair/HyperWorks;http://www.altairhyper works.com/),f i xed displacement constraints were imposed at the elements of the proximal tibial plateau.Similar to our experimental setup,27a 3 N compressive load was applied to the distal end of the tibia(Figure 1a).Assuming bone elements to be an elastic material with 20 GPa Young’s modulus and 0.33 Poisson’s ratio,8the strain f i eld was obtained using OptiStruct,a FEA linear solver in the HyperWork software package.The average strain of a 1×3 mm area on the medial–anterior surface that was 20%–40%distal from the tibial proximal end was compared with the strain measurement of a similar area in ourprevious studies.27,35Good agreement between the comparisons would validate the material properties and boundary conditions assigned to the whole-bone FEA model.

Segment biphasic model

The model elements were de f i ned as an isotropic porous linear elastic material(20 GPa modulus,0.33 Poisson’s ratio) saturated with interstitial f l uid(viscosity=0.001 Pa·s).8As reviewed previously,9there are three levels of porosity in the bone tissue:the large vascular pore(order 10μm),the lacunar–canalicular pores(order 1–0.1 μm),and intercollagen hydroxyapatite pores(order 1–10 nm).These intertwined pores make direct measurements of bone permeability quite challenging.The reported permeability varied from 10-12to 10-23m2,depending on the levels of pores probed.9,38–39Because the murine cortex is relatively thin with pores being predominantly smaller LCS ones,the permeability at the tissue(element)level was chosen to be in the lower range of the reported values associated with the smaller pores.In the model,the permeability was varied parametrically over three orders of magnitude(2.8×10-20, 2.8×10-21,or 2.8×10-22m2).The permeability was assumed to be isotropic and identical at both trabecular and cortical bones.Solute diffusivity was assumed to be isotropic for small strain cases as in loaded tibia.27

Measurements of LCS porosity in the murine cortical bone One key parameter for the above poroelastic FEBio modeling is the volume fraction of the LCS pores.Large variations of LCS porosity(1%–22%)were reported in literature using either two-dimensional or 3D imaging techniques.39We thus measured the volume fractions of lacuna and canaliculi in the murine cortical bone using in situ confocal microscopy and a protocol modi f i ed from a previous study.40In brief,femoral mid-shafts were dissected and sectioned into 0.3 mm-thick segments,f i xed in 10%formaldehyde for 24 h,polished to 0.05–0.10 mmthickness,dehydrated in ascending grades of ethanol, stained in sodium f l uorescein solution for 4 h,and mounted on an imaging chamber with a bottom cover glass.Highresolution 3D stacks(122 slices)of a f i eld of 512×256 pixels (pixel 0.199 μm;z-step=0.2 μm),which contained several lacunae and numerous canaliculi,were captured using a 60×oil objective and a Zeiss LSM 510 confocal microscope (Carl Zeiss Microscopy,LLC,Peabody,MA,USA).The raw images were then imported to the Volocity software package(PerkinElmer,Tempe,AZ,USA),where an adaptive threshold and a size-dependent segmentation scheme were applied sequentially to separate the larger lacunae from the smaller canaliculi.The surface area of the lacunae and the total volume of the canaliculi was then obtained. Unlike previous measurements,we quanti f i ed the volume fraction of the LCS pores by subtracting the cell body volume.Our previous transmission electron microscopy studies showed that within the lacunae there was a 0.49 +0.15 μm gap between the lacunar matrix wall and the cell body,and that the pericellular f l uid space occupied 79%of the total canalicular cross-sectional area.15These measures were used to calculate the volume fractions of the lacunar pores and the canalicular pores,respectively.The lacunar pore volume was the product of the total lacunar surface area and the lacunar gap.The total canalicular pore volume was a fraction(79%)of the total measured canalicular volume.Our result(shown later in the Results section)demonstrated a total lacunar and canalicular porosity(15.4%)in mouse,where the majority was contributed by canaliculi(87.5%).This value was adopted for all the simulations shown below.Please note that this value may represent an overestimation of the LCS porosity due to potential artifacts such as the partial volume effect and axial distortion associated with light microscopy.41

Simulation of FRAP experiments:validation of the segment model

As part of validation of the biphasic transport model,we simulated FRAP experiments as performed in our previous experiments27because they provided the most relevant experimental data in literature.One element(20 μm3), which was～30 μm below the anterior–medial periosteum and similar to the dimensions of a single lacuna,was chosen to be the photobleached lacuna(the center of the highlighted green area,Figure 1c).The immediate postphotobleaching tracer concentration of the photobleached element and those within the 90olaser cone below and above the photobleached element were set to be 0.2 due to the effects of photobleaching,whereas a concentration of 1.0 was assigned to the rest of elements in the model.Due to the relatively small hydraulic conductance and solute permeability of the mineralized bone tissue,as a f i rst estimation,impermeable boundary conditions(no f l uid and solute f l ux)were assumed for the periosteal surface and the bone cross-sections at the two ends of the model(Figure 1c).26The marrow cavity was assumed to have a constant pressure and tracer concentration due to the presence of interstitial osmotic/hydraulic pressures and the tracer-rich vasculature,respectively.A free draining solute/f l uid f l ux was assumed on the endosteal surface(Figure 1c).25Transport of solutes with diffusivity varying from 27 to 100 μm2·s-1was simulated as below.

Two FRAP experiments under non-loaded and loaded conditions27were simulated using the transport model.For the non-loaded condition,the model was run for a total duration of 36 s with a time step of 0.2 s.The time course of solute concentration in the photobleached element was obtained from the simulation results.This time course was predicted to be an approximately exponential function of time.16,27A transport rate Kdiff,was de f i ned as the slope of the curve of ln[(C?C0)/(Cb?C0)]vs time,where C was the concentration at time t,C0(=1)the concentration before photobleaching,and Cb(=0.2)the concentration immediately after photobleaching,Kdiffmeasured the speed of the tracer recovery,which was the reciprocal of the time constant for the exponential recovery.16Based on the experimental result of Kdiff(=0.017 s-1)for sodium fl uorescein,27a best- fi t diffusion coef fi cient was determined and compared with that obtained using the mathematical model developed previously.16For the loaded condition, we fi rst obtained the dynamic displacements on the two transverse surfaces that comprise the distal and proximal boundaries of the transport model(indicated by the dashed lines in Figure 1a and the purple surfaces in Figure 1c).These data were obtained from the wholebone model simulation under the 3N peak load(Figure 1b) with a 0.5 Hz sinusoidal waveform followed by 2 s resting periods(total 4 s for a cycle).Up to eight cycles of loading were simulated.To capture the enhanced convective transport into photobleached lacuna through both loading and unloading phases of the loading cycle,26the local solute concentration was superimposed with the two phases in sequence,as mixing in larger lacuna helped the entrapment of the tracer locally.25A total of eight loading cycles(32 s)were simulated,from which the transport enhancement(Kload/Kdiff)was obtained for sodium f l uorescein and compared with the experimental measurements.27A good agreement between the modeling results and the experimental data would justify the use of material properties assumed in the FEBio model.

Outputs from the segment model

Once validation of the transport model was con f i rmed,the spatiotemporal distributions of pore pressure and f l uid f l uxesat the tissue level were obtained from the FEA model under cyclical mechanical loading.Several locations in the midtransverse plane were chosen to show the results at both cortical and trabecular sites.

Outputs from the ultrastructural canalicular fl ow model

To accomplish the goal of predicting the load-induced mechanical stimulation on osteocytes located on the entire bone,the outputs from the above segmental transport model,which re fl ected the tissue-level measures, were converted onto the cellular level.The tissue-level fl uid fl uxes were scaled to the canalicular level by a factor of 6.5 based on the average LCS porosity(15.4%),that is,the canalicular-level fl uid fl ux would be 6.5 times of that at the tissue level predicted by the FEBio model.Using the Brinkman fl uid fl ow model for a single canaliculus developed by Weinbaum et al.,8we were able to predict the forces acting on the osteocytes including the shearing force and fl uid drag force.17The essential component of the Weinbaum model was the presence of fi bers/tethers spanning the annular fl uid space between the membrane of the cell process and the mineralized wall(Figure 1e). Recent studies have measured the radii of the canaliculi (160 nm)and osteocyte process(76 nm)in adult mouse bone,18and a fi ber spacing(center to center 14.3 nm).29,42The formula of the canalicular fl ow pro fi le,the shear stress,the shearing force,and drag force per unit length were given in the appendix of reference.43As the tissue-level fl ow was found to be sensitive to tissue permeability,the sensitivity of the canalicular fl ow velocity was also tested by varying the permeability over three orders of magnitude(2.8×10-20,2.8×10-21,and 2.8×10-22m2).

RESULTS

Model validations

Validating the whole-bone FEA model.The intact tibia bone was deformed by a combined mode of compression and bending under the 3 N compressive load applied at the ends.A tensile strain of～450 μ? was predicted on the relatively f l at anteromedial tibial surface around the FRAP site(30%distal from the proximal end),and compressive strains were mostly found on the posterior cortex(Figure 2).In this linear FEA model,the strain values were proportional to the assumed Young’s modulus.Our predicted tensile strain of 450 μ? at the region of interest (FRAP site)matched well with the experimentally measured data of～400 μ?.27Thus,we concluded that the assumed material properties of 20 GPa Young’s modulus and 0.33 Poisson’s ratio were justi f i ed for the subsequent transport modeling.Measuring the LCS porosity(an input to the segment model).Extensive LCS pores were labeled with high-intensity green fl uorescence in the high-resolution 3D stacks of confocal images(Figure 3a),and the images were thresholded and segmented into lacunae and canaliculi categories(based on size criteria)for calculation of the porosity(Figure 2b and c).The volume fractions of the lacunae and canaliculi were found to be 1.9%and 13.5%,respectively.The total LCS porosity (15.4%)was then used as an input to the segment model. This porosity was also used to scale the tissue-level fl uid fl ow predicted by the segment model to that at the canalicular level in the ultrastructural model.

Validating the segment biphasic model.For the nonloaded condition,the solute concentration at the FRAP site was shown to increase with time and the rate of increase varied with diffusivity(Figure 4a).The transport rate,K,shown as the slope of ln[(C?C0)/(Cb?C0)]vs time curve(Figure 4b),was slightly higher initially and gradually reached a steady state(constant slope).This steady-state transport rate(Kdiff)was nearly linearly proportional to the solute diffusivity(Kdiff=4.11×10-4×D, Figure 4c).In particular,a diffusivity of 31.8 μm2·s-1in the 3D porous model was found to best match the experimentally observed Kdiff=0.017 s-1of sodium f l uorescein (376 Da).27

This tissue diffusivity of 31.8 μm2·s-1was then used to simulate the loaded condition(Figure 5).From the time courses of solute concentration(Figure 5a)and the logarithm of the recovery at the FRAP site under loaded and non-loaded conditions(Figure 5b),the transport enhancement(Kload/Kdiff)was 1.24 for the 3N loading 3N,which fell within 1 s.d.above or below from the previously obtained experimental mean value(Kload/Kdiff=1.31±0.24).27In addition,we ran FRAP simulations with higher loads(5 N and 7 N).As anticipated, we observed faster f l uorescence recovery and greater transport enhancement as load magnitude increased (Figure 5).Taken together,these results provide strong evidence that supports the use of FEBio segment transport model to predict the pore pressure and f l uid f l ux in mechanically loaded bone.

Results from the segment transport model

If not stated otherwise,the results presented in this section were obtained using the following parameters:a peak force of 3 N,material properties of 20 GPa Young’s modulus,0.33 Poisson’s ratio,15.4%porosity,and permeability of 2.8×10-20m2.

Pore f l uid pressure f i eld.The pore f l uid pressure distribution(Figure 6)was obtained at t=3 s when the loading peaked at 3 N(Figure 1b).In general,negative pressures were found in the regions under tension and positive pressures in the compressed regions,and the pressure magnitude in the cortex increases with the distance to the neutral plane(Figure 6a).Negative pressures were found in the trabecular area adjacent the cortex under compression(location G shown in blue shades),suggesting a more complex local loading pattern there.The temporal changes of the f l uid pressure were shown duringone loading cycle in Figure 6b,where location A (corresponding to the FRAP site)experienced a peak pressure of-8 MPa,and a high 17 MPa pressure was found in location D.Pressures in the trabecular locations(G and H)were relatively smaller in magnitude. Reducing permeability by one or two orders of magnitude (2.8×10-21and 2.8×10-22m2),peak pressure was found to vary slightly(＜10%)or modestly(＜25%)at locations A(+7%and+15%),B(+6%and+15%),C(-9%and-17%), D(+7%and+18%),E(+4%and+14%),F(-6%and-17%), G(+8%and+24%),and H(-6%and-18%),respectively.

Fluid fl ow fi eld.The local distribution of load-induced fl uid fl ux could be obtained from the segment model.The local fl uid fl ow varied cyclically as a function of time and a snapshot of the fl uid fl ow fi eld at t=2.6 s is shown in Figure 7a.The fl ows for the surface elements are shown with the vectors,with the length indicating fl ow magnitude and arrow indicating the fl ow direction;and the fl ow magnitude for other elements are indicated with pseudocolors(Figure 7a).Overall,higher f l ow rates were found near the endosteal surfaces.The temporal pro f i les of f l ow magnitude at the selected locations are shown in Figure 7b.Comparing with f l uid pressure that dropped to zero after t=4 s in most locations(Figure 6b),f l uid f l ow at location B and C persisted till t=5 s(Figure 7b).Among all those selected locations,location C near the endosteal surface experienced the largest f l ux.Fluid f l ow was also found in the trabecular site(as shown in locations G and H),although the f l ux was relatively smaller than that in the adjacent endosteal cortical bone(Figure 7).

Results from the ultrastructural canalicular f l ow model

Scaling based on the LCS porosity(15.4%),the peak canalicular f l ux was predicted to be 6.5-fold higher than that at the tissue level,varying from 0.02 to 1.84μm·s-1(Table 1)for cortical sites(A–F)and trabecular sites(G and H)(Figure 6).Site-dependent variations were also found for the shear stress acting on the osteocyte cell process and the two fl uid-related stimulating forces such as fl uid shearing force and drag force(Table 1).Consistently, the fl uid drag force was approximately one order greater than the shear force.

DISCUSSION

The anabolic effects of mechanical forces have long been appreciated by the musculoskeletal research community.1,4Clinicians routinely recommend exercise and physical activity to patients when treating and managing osteoporosis,osteoarthritis,and other medical conditions(www.cdc.gov).Quantifying the mechanical stimulation,however,remains a challenge,because of the complex structure and inhomogeneous material properties of the bone tissue32and the different mechanical loading parameters associated with various exercise regimens.9–10With the advances of imaging techniques,we are now able to reconstruct the 3D structures of bones with suf fi cient details from the whole-bone level(on the order of 0.1–1 m), individual osteon/trabecula level(on the order of 0.2 mm), and down to the cellular level(on the order of 0.01 mm).39One striking feature of the bone tissue is the multiscaled, interconnected pore system housing the cellular components in bone,spanning the central marrow cavity, vascular channels,and the LCS.9These pores are saturated with interstitial fl uid,providing a continuous pathway for nutrient supply,waste removal,and exchange of signaling molecules.Furthermore,when bone is subjected to mechanical strains during exercise,it behaves as a stiff sponge:pore pressure builds up and fl uid is forced to fl ow within the bone matrix through the pores.9The osteocytes residing in the bone thus experience matrix deformation, fl uid pore pressure, fl uid shear force,as well as fl uid drag force due to the pericellular tethering fi bers.9These mechanical parameters have been quanti f i ed using various engineering techniques.Matrix deformations have been measured with strain gages and imaging correlation methods.3,5FRAP-based imaging techniques in combination with mathematical models provide measures of f l uid and solute f l ows at the canalicular level under loading,27allowing us to predict the various forces(shear and drag)that can trigger cellular responses.28–29These data greatly enhance our understanding of the mechanosensing and mechanotransduction processes in bone that underlie the anabolic power of physical activities.However,these measures are currently limited to regions close to bone surfaces.The objective of this study was to expand the mapping of mechanical stimulations to the inner portions of the bone in a multiscale manner.

This study demonstrated the feasibility of performing such comprehensive mapping using an image-based FEA platform.This approach is consisted of three coupled models at the whole bone,bone segment,and single canaliculus levels.Comparing with previous models,19–26the current model is unique by incorporating the ultrastructural model to predict loading effects(in terms of f l uid shear and drag forces)that are highly relevant to the functioning of osteocytes.The model is physically sound,where thesegment model is coupled with the whole-bone level through load-induced displacement fi elds,and the intrinsic LCS pores couple the segment model with the single canaliculus model.We also took advantage the available experimental data(on cortical bone)to extensively validate the model.The optically measured strain data27,35were used to compare with the model outputs,ensuring that the material properties and the deformation coupling scheme were appropriate.The FRAP diffusion and transport enhancement measures16,27,34were used to compare with the biphasic model in FEBio,enhancing our confi dence of the choices of model parameters and boundary conditions.Using the validated models,we showed that loading-induced cellular stimulations such as pore pressure, fl uid fl ow-induced shear stress and/or drag force acting on osteocytes could be mapped at both cortical and trabecular sites(Figures 6 and 7;Table 1).These quantitative data would help us to better understand which loading-induced parameters,including but not limited to matrix deformation,pore pressure,f l uid shear, and f l uid drag,correlate well with the spatial distribution of in vivo bone formation.

Table 1.Cellular stimulation forces

As with any model,our model has its limitations and assumptions.(1)The FE model was generated from a single biological sample,limiting the adoption of statistical analysis.As the goal of this study was to demonstrate the feasibility of modeling bone f l uid f l ow,this simpli f i ed approach was chosen to remove any potential variabilities from the bone geometry and thus allow the study to better focus on validating model parameters and coupling schemes.With the procedure being streamlined and computational power ever increasing,the methodologies described herein can be applied for multiple samples to account for variations on bone anatomy.(2)The scaling factor between tissue-level f l uid f l ow and the canalicularlevel f l uid f l ux was assumed to depend on the LCS porosity (15.4%),which was measured using confocal imaging.As noted earlier,this value may be an overestimation due to the scattering of f l uorescence signals and axial stretching in the point spread function.41Indeed,smaller porosity values(1%–5%)have been reported using methods based on electron microscopy and x-ray computation topography.39Were the LCS assumed to be 1.5%,the model outputs(Table 1)would be expected to be nearly one order of magnitude higher.(3)The permeability we used in the model(2.8×10-20m2)was based on our experimental measurement in dog bone.38Large variations(in several orders of magnitude)of permeability have been reported in the literature.39We also tested the sensitivity of the model outputs to permeability.As the permeability was reduced for one or two orders of magnitude(2.8×10-21and 2.8×10-22m2),f l uid f l ow velocity was found to decrease compared with the values presented here(that is,22%and 0.25%for locations A–C, and 69%and 15%for locations D–F,respectively).We thus conclude that accurate permeability measurement is the key to predict fl uid fl ow velocity in the model.(4)We assumed sealed boundaries in our segment model for faster convergence in solute concentration simulations.This idealized condition was not fully compatible with in vivo situation where periosteum was found to be permeable for fl uid and small solute.44–45Leaky permeability should be considered for future modeling.(5)We assumed a biphasic material with isotropic linear elastic solid phase,which has constant isotropic permeability and constant isotropic solute diffusion coef fi cients.Previous studies33indicated that the anisotropy of bone has an important role in the occurrence and distribution of the fl uid fl ow in bone,which should be quanti fi ed in future experimental and modeling studies.Despite these limitations,our model was demonstrated to serve as a promising platform that would allow indepth studies of local loading environment,which may help identify the important mechanical factors that drive bone’s response to loading and disuse in normal and disease conditions.

Acknowledgements

The study was supported by grants from NIH(P30GM103333 and RO1AR054385 to LW),a China CSC fellowship(to LF),and DOD W81XWH-13-1-0148(to XLL).

1 Duncan RL,Turner CH.Mechanotransduction and the functional response of bone to mechanical strain.Calcif Tissue Int 1995;57:344–358.

2 Klein-Nulend J,Bacabac RG,Mullender MG.Mechanobiology of bone tissue.Pathol Biol(Paris)2005;53:576–580.

3 Lanyon LE,Hampson WG,Goodship AE et al.Bone deformation recorded in vivo from strain gauges attached to the human tibial shaft. Acta Orthop Scand 1975;46:256–268.

4 Burr DB,Milgrom C,Fyhrie D et al.In vivo measurement of human tibial strains during vigorous activity.Bone 1996;18:405–410.

5 Nicolella DP,Moravits DE,Gale AM et al.Osteocyte lacunae tissue strain in cortical bone.J Biomech 2006;39:1735–1743.

6 Owan I,Burr DB,Turner CH et al.Mechanotransduction in bone: osteoblasts are more responsive to f l uid forces than mechanical strain. Am J Physiol 1997;273:C810–C815.

7 You J,Yellowley CE,Donahue HJ et al.Substrate deformation levels associated with routine physical activity are less stimulatory to bone cells relative to loading-induced oscillatory f l uid f l ow.J Biomech Eng 2000;122:387–393.

8 Weinbaum S,Cowin SC,Zeng Y.A model for the excitation of osteocytes by mechanical loading-induced bone f l uid shear stresses.J Biomech 1994;27:339–360.

9 Fritton SP,Weinbaum S.Fluid and solute transport in bone:f l owinduced mechanotransduction.Annu Rev Fluid Mech 2009;41:347–374.

10 Bonewald LF.The amazing osteocyte.J Bone Miner Res 2011;26:229–238.

11 Robling AG,Niziolek PJ,Baldridge LA et al.Mechanical stimulation of bone in vivo reduces osteocyte expression of Sost/sclerostin.J Biol Chem 2008;283:5866–5875.

12 You L,Temiyasathit S,Lee P et al.Osteocytes as mechanosensors in the inhibition of bone resorption due to mechanical loading.Bone 2008;42: 172–179.

13 Dallas SL,Prideaux M,Bonewald LF.The osteocyte:an endocrine cell and more.Endocr Rev 2013;34:658–690.

碑志是中國古代最為常見常用的石刻文字。朱熹所撰書的石刻文字，主要有兩類：碑志文和摩崖題刻。這些石刻文字不僅豐富了石刻的斯文內(nèi)涵，也提升了石刻的文化境界。

14 Piekarski K,Munro M.Transport mechanism operating between blood supply and osteocytes in long bones.Nature 1977;269:80–82.

15 Lai X,Price C,Modla S et al.The dependences of osteocyte network on bone compartment,age,and disease.Bone Res 2015;3:15009.

16 Wang L,Wang Y,Han Y et al.In situ measurement of solute transport in the bone lacunar-canalicular system.Proc Natl Acad Sci USA 2005;102: 11911–11916.

17 You L,Cowin SC,Schaf f l er MB et al.A model for strain ampli f i cation in the actin cytoskeleton of osteocytes due to f l uid drag on pericellular matrix.J Biomech 2001;34:1375–1386.

18 You LD,Weinbaum S,Cowin SC et al.Ultrastructure of the osteocyte process and its pericellular matrix.Anat Rec A Discov Mol Cell Evol Biol 2004;278:505–513.

19 Goulet GC,Coombe D,Martinuzzi RJ et al.Poroelastic evaluation of fl uid movement through the lacunocanalicular system.Ann Biomed Eng 2009;37:1390–1402.

20 Goulet GC,Hamilton N,Cooper D et al.In fl uence of vascular porosity on fl uid fl ow and nutrient transport in loaded cortical bone.J Biomech 2008;41:2169–2175.

21 Gururaja S,Kim HJ,Swan CC et al.Modeling deformation-induced fl uid fl ow in cortical bone's canalicular-lacunar system.Ann Biomed Eng 2005; 33:7–25.

23 McCarthy ID,Yang L.A distributed model of exchange processes within the osteon.J Biomech 1992;25:441–450.

24 Srinivasan S,Gross TS.Canalicular f l uid f l ow induced by bending of a long bone.Med Eng Phys 2000;22:127–133.

25 Wang L,Cowin SC,Weinbaum S et al.Modeling tracer transport in an osteon under cyclic loading.Ann Biomed Eng 2000;28:1200–1209.

26 Zhou X,Novotny JE,Wang L.Modeling f l uorescence recovery after photobleaching in loaded boneotential applications in measuring f l uid and solute transport in the osteocytic lacunar-canalicular system.Ann Biomed Eng 2008;36:1961–1977.

27 Price C,Zhou X,Li W et al.Real-time measurement of solute transport within the lacunar-canalicular system of mechanically loaded bone: direct evidence for load-induced f l uid f l ow.J Bone Miner Res 2011;26: 277–285.

28 Wang B,Zhou X,Price C et al.Quantifying load-induced solute transport and solute-matrix interaction within the osteocyte lacunarcanalicular system.J Bone Miner Res 2013;28:1075–1086.

29 Wang B,Lai X,Price C et al.Perlecan-containing pericellular matrix regulates solute transport and mechanosensing within the osteocyte lacunar-canalicular system.J Bone Miner Res 2014;29:878–891.

30 Oxlund H,Ejersted C,Andreassen TT et al.Parathyroid hormone(1-34) and(1-84)stimulate cortical bone formation both from periosteum and endosteum.Calcif Tissue Int 1993;53:394–399.

31 Par f i tt AM,Mathews CH,Villanueva AR et al.Relationships between surface,volume,and thickness of iliac trabecular bone in aging and in osteoporosis.Implications for the microanatomic and cellular mechanisms of bone loss.J Clin Invest 1983;72:1396–1409.

32 Deligianni DD,Apostolopoulos CA.Multilevel f i nite element modeling for the prediction of local cellular deformation in bone.Biomech Model Mechanobiol 2008;7:151–159.

33 Steck R,Niederer P,Knothe Tate ML.A f i nite difference model of loadinduced f l uid displacements within bone under mechanical loading. Med Eng Phys 2000;22:117–125.

34 Li W,You L,Schaf f l er MB et al.The dependency of solute diffusion on molecular weight and shape in intact bone.Bone 2009;45:1017–1023.

35 Price C,Li W,Novotny JE et al.An in-situ f l uorescence-based optical extensometry system for imaging mechanically loaded bone.J Orthop Res 2010;28:805–811.

36 Maas SA,Ellis BJ,Ateshian GA et al.FEBio:f i nite elements for biomechanics.J Biomech Eng 2012;134:011005.

37 Mauck RL,Hung CT,Ateshian GA.Modeling of neutral solute transport in a dynamically loaded porous permeable gel:implications for articular cartilage biosynthesis and tissue engineering.J Biomech Eng 2003;125: 602–614.

38 Gardinier JD,Townend CW,Jen KP et al.In situ permeability measurement of the mammalian lacunar-canalicular system.Bone 2010;46: 1075–1081.

39 Cardoso L,Fritton SP,Gailani G et al.Advances in assessment of bone porosity,permeability and interstitial f l uid f l ow.J Biomech 2013;46: 253–265.

40 Ciani C,Doty SB,Fritton SP.An effective histological staining process to visualize bone interstitial f l uid space using confocal microscopy.Bone 2009;44:1015–1017.

41 Sharma D,Ciani C,Marin PA et al.Alterations in the osteocyte lacunarcanalicular microenvironment due to estrogen de f i ciency.Bone 2012;51: 488–497.

42 Wijeratne SS,Martinez JR,Grindel BJ et al.Single molecule force measurements of perlecan/HSPG2:A key component of the osteocyte pericellular matrix.Matrix Biol 2016;50:27–38.

43 Jing D,Baik AD,Lu XL et al.In situ intracellular calcium oscillations in osteocytes in intact mouse long bones under dynamic mechanical loading.FASEB J 2014;28:1582–1592.

44 Evans SF,Parent JB,Lasko CE et al.Periosteum,bone's"smart"bounding membrane,exhibits direction-dependent permeability.J Bone Miner Res 2013;28:608–617.

45 Lai X,Price C,Lu XL et al.Imaging and quantifying solute transport across periosteum:implications for muscle-bone crosstalk.Bone 2014;66: 82–89.

This work is licensed under a Creative Commons Attribution 4.0 International License.The images or other third party material in this article are included in the article’s Creative Commons license,unless indicated otherwise in the credit line;if the material is not included under the Creative Commons license,users will need to obtain permission from the license holder to reproduce the material.To view a copy of this license,visit http://creativecommons.org/licenses/by/4.0/

?The Author(s)2016

Research(2016)4,16032;

10.1038/boneres.2016.32;published online:27 September 2016

1Department of Mechanical Engineering,University of Delaware,Newark,DE,USA and2School of Mechanical Engineering,Nanjing University of Science and Technology,Nanjing,China

Correspondence:Liyun Wang(lywang@udel.edu)

*The two authors contributed equally to this work.

Received:16 May 2016;Revised:15 August 2016;Accepted:21 August 2016