Three-Scale Multiphysics Modeling of Transport Phenomena within Cortical Bone

Bone tissue can adapt its properties and geometry to its physical environment. This ability is a key point in the osteointegration of bone implants since it controls the tissue remodeling in the vicinity of the treated site. Since interstitial fluid and ionic transport taking place in the fluid compartments of bone plays a major role in the mechanotransduction of bone remodeling, this theoretical study presents a three-scale model of the multiphysical transport phenomena taking place within the vasculature porosity and the lacunocanalicular network of cortical bone. These two porosity levels exchange mass and ions through the permeable outer wall of the Haversian-Volkmann canals. Thus, coupled equations of electrochemohydraulic transport are derived from the nanoscale of the canaliculi toward the cortical tissue, considering the intermediate scale of the intraosteonal tissue. In particular, the Onsager reciprocity relations that govern the coupled transport are checked.


Introduction
Bone remodeling corresponds to a process by which bone tissue gradually alters its morphology and properties in order to adapt to its external environment [1,2].This ability is a key point of the success of different protocols in bone tissue engineering.It has thus been demonstrated that bone remodeling around implants is essential to obtain a good osteointegration of the implant [3].Similarly, in bone tissue engineering, when grafting a scaffold to fill a bone loss, there is increased bone remodeling within the treated site until the degradation of the entire scaffold [4].
That is why the study of the mechanisms governing bone remodeling process is a very important topic in bone biomechanical engineering.Interstitial fluid and ionic transport taking place in the fluid compartments of bone is thought to play a major role in the bone ability to transmit physical stimuli toward bone cells and thus trigger the remodeling process [5,6].
Cortical bone, which corresponds to the dense tissue located at the periphery of bones, is mostly made of cylindrical structures called osteons.These cylinders are crossed by a hierarchical porous network saturated by a fluid with both mechanical and chemical functions.Axial Havers and radial Volkmann canals (HVC), few tens of microns in diameter, contain vasculature and nerves and provide a fast way to transmit chemomechanical information between different osteons.Local messaging takes place inside the smaller lacunocanalicular porosity (LCP), a dense network of cavities (lacunae) and pseudocylindrical canals (canaliculi, few hundred nanometers in radius) crossing the osteonal matrix.The LCP hosts the cellular network of the bone mechanosensing cells, the osteocytes, and is connected with the HVC.
Since the experimental description of the fluid flow within the cortical tissue is still a very challenging topic [7,8], theoretical approaches are often carried out to understand the in vivo phenomena that govern bone behavior.Mathematical descriptions of bone fluid flow often only consider the osteonal scale to describe the transport within the LCP [9][10][11][12] or sometimes involve a biporous treatment of this medium [13][14][15], even if recent works tend to show that the nanopores inside the collagen-apatite matrix may also participate in the fluid movement inside bone volume [16,17].
In geoscience, fluid and ionic transport in monoporous media is a classical application of upscaling methods to derive effective transport coefficients [18][19][20][21][22]. Similarly, hydrochemical transport phenomena in the LCP were extensively studied by our group [23,24] highlighting their effects on bone physiology [25,26].This paper stems from those results and extends them to the osteonal scale, focusing on the fluid and ionic transport along the osteon.An asymptotic homogenization technique is used to upscale a coupled description of the transport phenomena toward the tissue scale [23,24].The outcome is an explicit description of the velocity and the ionic flux along the osteon which turns out to be ruled by generalized forms of the Darcy law and macroscopic convectiondiffusion equations.
To the best of our knowledge, this theoretical work is the first attempt to describe the multiphysical transport phenomena taking place within the HVC and the LCP at once.These two fluid compartments exchange mass and ions through the permeable outer wall of the HVC.The model developed in this paper accounts for these exchanges by identifying the normal fluxes of mass and ions through the outer walls of the HVC with the longitudinal fluxes inside the canaliculi.Furthermore, a rewriting of the macroscopic equations when convection and diffusion mechanisms are comparable resulted in checking the Onsager reciprocal relations linking the fluxes and driving forces that govern the coupled transport phenomena.

Multiscale Structure of the Cortical Tissue
In Figure 1, the multiscale structure of cortical bone is presented.From the organ to the cell, the hierarchical description of cortical bone reads the following: (a) cortical tissue scale (10 −2 m): cortical bone is made of cylinders called osteons that are interconnected by macropores (axial Havers and radial Volkmann canals) containing the vasculature and the nerves; (b) osteon scale (10 −3 m): the osteon is made of concentric lamellae and contains a porous network of lacunae interconnected thanks to thin channels (canaliculi); (c) LCP scale (10 −5 m): located between the lamellae, each lacuna hosts a mechanosensitive cell called osteocyte that develops dendritic protrusions inside the canaliculi to form a stellar network; (d) canalicular pore scale (10 −8 m): the canalicular pore space, which is roughly annular between the collagenapatite matrix and the cell dendrite surface, is partially occupied by a pericellular matrix made of charged proteins (mainly glycosaminoglycans).Interstitial fluid transport at this scale should take into account diffusive, convective, pericellular friction and electrochemical effects.
In this paper, we consider the interstitial fluid transport in two levels of these nested porous networks, namely, the Havers-Volkmann channels and the lacunocanalicular network.Thus, in our three-scale homogenization process, the microscale corresponds to the canalicular pore scale (Figure 1(d); ℓ  ∼ 10 −8 m), the mesoscale corresponds to the intraosteonal scale of the LCP and the pores correspond to the Haversian-Volkmann channels (Figures 1(b) and 1(c); ℓ  ∼ 10 −5 m), and the macroscale corresponds to the cortical tissue scale (Figure 1(a); ℓ  ∼ 10 −2 m).As a consequence, the micro-to-meso and meso-to-macro characteristic lengths ratios are comparable: ℓ  /ℓ  ∼ ℓ  /ℓ  ∼  ≡ 10 −3 .Hereinafter, sub-and superscripts , , and  will denote quantities at the micro-, meso-, and macroscale, respectively.

From the Microscale to the Mesoscale
At the microscale, the representative volume element (RVE) is a bone fraction surrounding the canaliculus (see Figure 1(d)).Canaliculi are seen as two concentric cylinders.The interstitial fluid, which is assumed to be an incompressible Newtonian monovalent electrolyte (typically Na-Cl), occupies the annular space between the outer canalicular wall and the inner osteocyte process membrane.The cationic and anionic concentrations are noted  + and  − , respectively.Under bulk conditions, electroneutrality reads  + =  − =   , where   stands for the salinity.
Accompanying the pressure induced hydraulic flow, other electrochemical phenomena are generated by the electrolyte solution movement [27].Indeed, due to the pore surface charge, double layers develop in the pore volume to maintain electroneutrality.When advected by the fluid flow, the mobile charge population of the double layer generates the macroscopically observed streaming potentials.The spatial variability of these potentials generates electroosmotic flow [28].
Thus, the microscopic model combines an electrostatics balance to describe the double layers, Nernst-Planck equations for the ionic species transport that include possible electromigration effects, and a Stokes equation considering the effects of the Coulombic body force.Moreover, the tethering pericellular fibers can sensibly reduce the permeability of the canalicular space.This friction effect is represented by a Brinkman-like term involving a pericellular permeability [29].

Statement of the Transport Model at the Microscale.
This subsection briefly recalls the main results of the two-scale treatment of multiphysical transport within the lacunocanalicular system as presented in [24].
First, electrostatics is expressed thanks to the Poisson-Boltzmann equation: In this equation, two reduced electric potentials appear: the reduced streaming potential   and the reduced double layer potential   .Moreover the Debye length    roughly characterizes the thickness of the electrical double layer.Electric flux continuity at the pore surfaces is the associated boundary condition.
Then, the Coulombic force effects can be split to make electrochemical driving effects appear in addition to the hydraulic pressure effect.Fluid flow has thus threefold nature: the Poiseuille velocity k   in response to the gradient of the fluid pressure   , a chemoosmotic velocity k   in response to the gradient of the salinity   , and an electroosmotic velocity k   in response to the gradient of the streaming potential   .That is why, introducing the fluid dynamic viscosity   and the pericellular permeability   , the three velocities are governed by the following equations: ( The coefficients   • (with • = , , ) on the right-hand side of these equations read [26] where  is the gas constant and  the temperature.The total fluid velocity vector at the micropore scale is thereby the sum of the three velocities . Moreover, at the canalicular lateral surfaces, no-slip boundary conditions are assumed for each elementary velocity.Note that this assumption to consider no-slip condition for each coupled velocity is necessary to perform the homogenization process as explained in [29].
Finally, noting  ± the ionic diffusion parameters, the ionic fluxes at the microscale J  ± = − ± exp(∓  )(∇  ±   ∇  ) are due to diffusion and electromigration effects.Then, the Nernst-Planck convection-electrodiffusion equations read then Ionic exchanges are possible at the pore surface and quantified thanks to exchanges coefficients   ± :

Description of the Transport Phenomena at the Mesoscale of the Lamellar Tissue.
The first upscaling procedure to derive the transport phenomena at the mesoscale from the microscopic description is based on the asymptotic periodic homogenization process as proposed in [23,30,31].The two-level homogenization procedure for the same Poisson-Boltzmann equation ( 1) and Brinkman equations (2) can be found in [23,29], whereas the upscaling of the Nernst-Planck equations is given in [24].In particular, it was shown in these references that the fluid pressure   , the salinity   , and the streaming potential   =   do not vary at the microscale whereas the double layer potential   is purely microscopic.
At the mesoscale, the average fluid velocity along the canalicular network V  (the superscript  referring to the mesoscale) can be described by a modified Darcy law: Thus, the elementary velocities at the mesoscale are where K   , K   , and K   represent the effective permeability tensors quantifying the Poiseuille, osmotic, and electroosmotic effects at the mesoscale.The expressions of these tensors are detailed in our previous works (see, e.g., [23]).Moreover, ∇  represents the space derivative operator with respect to the coordinates of the mesoscopic level whereas ⟨⋆⟩  represents the averaging operator on the microscopic RVE.
Finally, to derive the effective ionic fluxes at the mesoscale ⟨J  ± ⟩  , we have to consider two cases.On the one hand, if the diffusive and convective effects are of the same order of magnitude, the canalicular Péclet number   ≡ (1) and we obtain The complete expression of the effective diffusion tensor at the mesoscale D  ± is detailed in [24].This quantity combines textural and electrical effects.These ionic fluxes along the canalicular system can then be separated into convective and electrodiffusive parts J  ± and J  ± : where On the other hand, if diffusion is the main transport process at the canalicular scale, the Péclet number is small and the convective term J  ± vanishes so that

From the Lamellar Tissue to the Cortical Tissue
The purpose of this section is now to describe the flow phenomena and chemical transport at the macroscale of the cortical tissue.At this scale, the porous network corresponds to the Havers-Volkmann canals (named Haversian pores hereafter), whereas the effective solid phase of the porous medium corresponds to the collagen-apatite tissue and the microporosity of the lacunocanalicular system.Hydroelectrochemical phenomena occurring across the canalicular network are nevertheless integrated in the description of transport phenomena across Haversian porosity.

Statement of the Transport Model at the Haversian Pore
Scale.Since capillaries and nerves are present in the central part of Haversian pores [32], the interstitial fluid flows in the pseudo-annular space between the outer wall of the Haversian canal and the walls of the capillaries and nerves.Again, the pore geometry can be roughly represented by two concentric cylinders representing the Haversian wall (outer cylinder) and capillary and nerves (inner cylinder).
The walls of the Haversian pores are permeable [33], allowing exchanges of matter and ions with the blood vessels (inner wall, indexed vas) and with the lacunocanalicular network (outer wall, indexed LCP).To use again periodic homogenization, we assume that a 3D periodic representative cell can be defined.Note that the possibility to obtain a reliable elementary volume at the scale of cortical tissue is a tricky point in bone biomechanics due to the high spatial heterogeneity of this tissue.
Similarly to the multiphysical approach proposed at the canalicular scale, the description of the transport phenomena at the scale of the Haversian porosity should take into account both hydraulic and electrochemical effects.
The electric flux continuity at the walls completes the electrostatics problem.

Fluid Movement with Coulombic
Force.Since the friction effect of the pericellular matrix is no more present at the mesoscale of the Haversian porosity, the coupled fluid flow is now represented by the following Stokes equations: where the coefficients   • are the same as before when transposed at the mesoscale: Note that these equations correspond to (2) considering an infinite fiber permeability (  → ∞).Again, the velocity in the Haversian porosity k  is the sum of three contributions: Noting n and t the outer normal and longitudinal tangent unit vectors, the boundary conditions at the walls of the Haversian porosity are given hereafter.
(i) Tangent no-slip condition: (ii) Exchanges with the vasculature or the LCP: where • = ,, and ⋆ = vas, LCP depending on the location of the wall.The filtration velocities across the blood vessels membranes V vas

•
have to be given.The filtration velocities at the outer wall of the pore V LCP • are equal to the velocities ⟨k  • ⟩  determined from the previous upscaling procedure in (7).
Note finally that pressure continuity between the Haversian and canalicular pores should be considered [34] and that the incompressibility of the fluid should also be assumed for each elementary velocity: 4.1.3.Ionic Transport.As at the microscale, the ionic flux vector at the mesoscale of the vascular porosity At the inner pore walls, the normal fluxes J  ± ⋅n are equal to the ionic fluxes through the blood vessels walls J vas ± ⋅ n (which has to be provided): At the outer pore walls, the normal fluxes J  ± ⋅ n are equal to the fluxes derived from the canalicular network analysis ⟨J  ± ⟩  ⋅ n (see (9) or (10) according to the value of the canalicular Péclet number):

Homogenization to Reach the Macroscale.
The homogenization procedure consists in the nondimensional writing of the problems, the asymptotic expansion, the collection of the slow variables, and the proposition of the closure problems (see, e.g., [23]).Here, this meso-to macroscale upscaling stage is often similar to the one performed to get to the mesoscopic description from the microscale.Thus, we will briefly present the similarities and be more specific about the differences.

Electrostatics. The homogenization process of the
Poisson-Boltzmann equation ( 11) is the same as the one proposed to perform the passage between the micro-and mesoscale.

Treatment of the Fluid Flow.
Again, the homogenization procedure of the fluid transport equation and the mass conservation are similar as previously done considering   → ∞.Only the boundary conditions should be treated in a different way.Indeed, the effects of the hydraulic exchanges between canalicular and Haversian pores may vary depending on the scale ratios ℓ  /ℓ  and ℓ  /ℓ  [34,35].Here, if the no-slip condition of (15) remains the same, according to the arguments of [36], the canalicular permeability has to be scaled by  2 , giving the following scaling rules of the nondimensional exchange conditions of ( 16): with ⋆ = vas, LCP depending on the considered pore surface.

Treatment of the Ionic
Transport.First, we consider that the convective and diffusive driving effects can be compared; that is to say that   ≡  (1).Analogously to the results of the previous paragraph, if the upscaling of the Nernst-Planck equations ( 18) is very close to the one performed previously (see ( 8)), only the specificity of the nondimensional writing of the boundary conditions ( 19) and ( 20) should be here discussed.
On the one hand, remembering the decomposition of ⟨J  ± ⟩  in ( 9), the nondimensional writing of the boundary equations ( 20) at the outer wall of the Haversian pores reads where the nondimensional convective flux numbers   ± and diffusive flux numbers   ±  resulting from the LCP transport are introduced.To propose scaling rules of these parameters, we have to estimate their order of magnitude.A reference value of the convective flux could be obtained from the lacunocanalicular porosity   , the reference fluid velocity in the LCP V  ref , the reference length of the mesoscale   ref , and the reference ionic diffusion coefficients  ref as follows: Since we have V  ref ≡ 10 −7 m⋅s −1 [7],   ≡ 10 −3 [37],   ref ≡ 10 −5 (see Figure 1(c)), and  ref ≡ 10 −9 m 2 ⋅s −1 (sodium bulk diffusivity in water), we obtain   ± ≡  2 .Furthermore, a reference value of the diffusive flux reads immediately Then, at the outer wall of the Haversian pores, the exchanges conditions of (22) read On the other hand, at the vasculature wall, if we consider that the ionic exchange term with vasculature is mainly due to convective terms, we have no dimensional conditions involving  2 :

Slow Variables.
Through the homogenization process, we recover that the fluid pressure   , the salinity   , and the streaming potential   =   are macroscopic fields.

Macroscopic Model at the Scale of Cortical Tissue.
In this subsection, the macroscopic equations of the model are given.We note ⟨⋆⟩  the volume average on the representative cell used to pass from the meso-to the macroscale and ⟨⋆⟩   the solid-fluid interface average (integral over the surface per unit of volume) of a given flux.
At the macroscale, the Poisson-Boltzmann problem is no more visible.However, it still has to be solved inside the representative volume since the double layer potential consequences do exist in the transport phenomena.Moreover, the Darcy fluid velocity k  is expressed through a coupled law as obtained before: where ∇  is the macroscopic space derivation operator and K  • are the effective coupled permeability tensors (with • = , , ).Their expressions can be derived from cell problems similar to those presented in [23].
Finally, the ionic transport at the macroscopic scale involves submacroscopic exchange terms where the effective velocity k * = k  −   , which is the difference between the fluid velocity in the mesopores k  and a chemical accumulation velocity term   .This chemical accumulation velocity represents the effects due to the downer scale.Furthermore, the complete expressions of the effective diffusion tensors at the macroscale D M ± are similar to those proposed in [24].

Discussion on the Onsager Reciprocal Relations
The Onsager reciprocal relations allow highlighting the interactions between transport phenomena occurring simultaneously in a thermodynamic system.In this section, the question of the validity of the Onsager reciprocal relationships is examined for the above problems.This analysis is based on a study of transport phenomena in coupled clayey media [38].

Discussion on the Fluid Flow in the Haversian Porosity.
The coupled Darcy law (27) is the average expression of a generalized Hagen-Poiseuille law expressing the fluid velocity at the Haversian pore scale k  : where   • are conductivity parameters quantifying the coupled contributions in the Haversian porosity.These transport parameters at the Haversian pore scale are obtained from the cellular problem appearing in the meso-to-macro homogenization process.Indeed, the classical closure statements for coupled flow are (see [23]) Introducing these equations in the Stokes equations ( 12), we obtain the cell problems: where the   • are those defined by (13).These cell problems are completed by introducing the closure statements in the previous boundary conditions (21).The coupled permeability parameters at the macroscale are then determined by an averaging procedure: 29) can thus be rewritten remembering that the definition of the reduced streaming potential involves the Faraday constant , the gas constant , and the temperature , so that   =   /, and noticing that ∇  = (  /)∇ ln   : The averaged version of this equation provides thus an alternative expression of the Darcy law (27):

Macroscopic Ionic Flux and Electric
Current.After having reformulated the fluid flow at the macroscale, we then intend to derive the macroscopic chemical flux and electric current.First, the sum of the cationic and anionic transport (28) reads This equation can be reformulated to obtain the conservation of the chemical flux: where the chemical flux J M is defined by Similarly, the difference of the cationic and anionic transport equations (28) gives the conservation of the electric current: where the electric flux I M is

Matrix Form of the Multiphysics Transport within Cortical
Tissue.If we assume that the chemical accumulation effects can be neglected, the effective velocity appearing in ( 28) is equal to the local fluid velocity in the Haversian porosity k * = k  and the use of ( 32) in (36) and (38) gives Finally, from the macroscopic Darcy law (33) and the expressions of the chemical flux (39) and of the electric current (40) we can derive the matrix form of the multiphysics transport phenomena at the cortical scale: The diagonal terms of the matrix represent the direct parameters whereas the off-diagonal terms quantify the coupling phenomena.Their expressions are then The Onsager reciprocal relations would require this matrix to be symmetric.This point is checked in Appendix by showing that L  = L  , L  = L  , and L  = L  .

Concluding Remarks
The use of modeling approaches to mimic the in vivo behavior of bone tissue is a common alternative to the difficulty to perform convenient in vivo experiments.When adopting a theoretical approach, it is necessary to find a way to validate the model results.In bone poromechanics, this validation is often performed through a method correlating the macroscopic streaming potentials induced by the stressgenerated fluid flow within the lacunocanalicular pores.This idea was first proposed by [39] and was later used by several other groups [10,23,40].However, all these studies are based on strong assumptions.First, they involve a two-scale approach connecting phenomena occurring in the canaliculi to phenomena at the intermediate scale of osteons or osteonal tissue whereas experimental in vivo measurements are performed at the organ scale (see, e.g., [41]).Furthermore, these studies consider homeostasis to neglect chemical gradient effects and put forward the Onsager reciprocity relations to obtain the streaming potential distribution from the upper scale fluid velocity.When steady state is reached, there is no net charge transfer, implying that the electric current is zero.This is equivalent to setting the Ohmic and convective currents equal and opposite.Thus, the streaming potential can be linearly connected to the macroscopic pressure field.
The present study makes it possible to reach the macroscopic scale of bone tissue and shows the validity of the Onsager reciprocity relations.
A multiscale model of the transport phenomena occurring at different scales of cortical tissue was developed adopting two successive upscaling procedures.It takes into account the existing hydroelectrochemical phenomena occurring at the different scales of pores (Haversian porosity and canaliculi).
Furthermore, a rewriting of the macroscopic equations when convection and diffusion mechanisms are comparable resulted in checking Onsager reciprocal relations linking the coupled quantities that govern the coupled transport phenomena.This approach thus reinforces the seminal idea of [39] to use macroscopic streaming potentials measurements to validate poromechanical models of bone that are extensively used in bone biomechanical engineering.

Recovery of the Onsager Reciprocity Property
If u is a quantity, the divergence theorem gives a link between the volume and surface average values.Moreover, due to the no-slip boundary condition which infers that   • = 0 on the pore surface, a strong simplification is obtained: