Xiran Cai
email: xirancai@stanford.edu

Laura Peralta

Pierre-Jean Gouttenoire

Christine Chappard

Renald Brenner

Cécile Olivier

Françoise Peyrin

Didier Cassereau

Pascal Laugier

Quentin Grimal

Cortical bone characterization using ultrasonic guided waves measured in axial transmission View project

Keywords: Bone microstructure, Elasticity, Homogenization, Resonant ultrasound spectroscopy

With ageing and various diseases, the vascular pore volume fraction (porosity) in cortical bone increases, and the morphology of the pore network is altered. Cortical bone elasticity is known to decrease with increasing porosity, but the effect of the microstructure is largely unknown, while it has been thoroughly studied for trabecular bone. Also, popular micromechanical models have disregarded several microarchitectural features, idealizing pores as cylinders aligned with the axis of the diaphysis. The aim of this paper is to quantify the relative effects on cortical bone anisotropic elasticity of porosity and other descriptors of the pore network microarchitecture associated with pore number, size and shape. The five stiffness constants of bone assumed to be a transversely isotropic material were measured with resonant ultrasound spectroscopy in 55 specimens from the femoral diaphysis of 29 donors.

The pore network, imaged with synchrotron radiation X-ray micro-computed tomography, was used to derive the pore descriptors and to build a homogenization model using the Fast Fourier Transform (FFT) method. The model was calibrated using experimental elasticity.

A detailed analysis of the computed effective elasticity revealed in particular that porosity explains most of the variations of the five stiffness constants and that the effects of other microarchitectural features are small compared to usual experimental errors. We also evidence that modeling the pore network as an ensemble of cylinders yields biased elasticity values compared to predictions based on the real micro-architecture. The FFT homogenization method is shown to be particularly efficient to model cortical bone.

Introduction

Cortical bone is a porous material which forms the outer shell of bones. Pores are observed at different length scales. The main component of the porosity is the so-called vascular porosity network [START_REF] Cooper | Cortical bone porosity: What is it, why is it important, and how can we detect it?[END_REF], constituted of quasi-cylindrical channels (Haversian canals, approximately 50 µm in diameter) oriented along the long axis of the bones and interconnected to each other by a network of smaller transverse channels (Volkmann's canals).

Larger pores formed by the removal of old bone by osteoclasts (from 50 µm to 250 µm [START_REF] Britz | The relation of femoral osteon geometry to age, sex, height and weight[END_REF]) are also present in bone and contribute to the vascular pore volume fraction depending on the balance of the remodeling activity. The vascular pores host blood vessels and nerves.

Smaller pores of characteristic dimensions less than 10 µm (as those hosting bone cells) also contribute to the overall porosity. The volume fraction of the vascular porosity, hereafter referred to as porosity, typically varies between a few percent and 20 %. Porosity and the morphology of the pore network can be regarded as a signature of the physiological processes regulating bone mass which are altered in the elderly and in pathologies such as osteoporosis [START_REF] Bala | Role of cortical bone in bone fragility[END_REF] or osteogenesis imperfecta [START_REF] Imbert | Microstructure and compressive mechanical properties of cortical bone in children with osteogenesis imperfecta treated with bisphosphonates compared with healthy children[END_REF]. An increased cortical porosity can be a consequence of an increased number of canals, an increased canal size, or both [START_REF] Cooper | Age-dependent change in the 3d structure of cortical porosity at the human femoral midshaft[END_REF]. At low porosity values, Haversian canals and small resorption cavities are more or less regularly distributed on a hexagonal lattice [START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -ii. measurements of elastic properties and microhardness[END_REF][START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -i. theoretical considerations for hexagonal symmetry[END_REF]. As the porosity increases, more irregular and large resorption cavities are found, resulting in a less organized network with a disparity of pore dimensions.

Understanding the determinants of cortical bone elastic properties is important because the bone cortex bears a substantial amount of the mechanical loads applied to the skeleton.

Millimeter-scale (mesoscale) elasticity [START_REF] Grimal | A determination of the minimum sizes of representative volume elements for the prediction of cortical bone elastic properties[END_REF] of cortical bone is known to strongly depend on porosity [START_REF] Currey | The effects of strain rate, reconstruction and mineral content on some mechanical properties of bovine bone[END_REF][START_REF] Dong | The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity[END_REF][START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF]. It was suggested that not only the porosity but also the changes in pore network micro-architecture may impact elasticity during growth and ageing [START_REF] Bala | Pore network microarchitecture influences human cortical bone elasticity during growth and aging[END_REF]. As a result, strain levels in a bone subjected to mechanical loading are partially driven by the pore volume fraction. These strains in return mediate bone homoeostatic activity and bone resistance to fracture [START_REF] Klein-Nulend | Mechanobiology of bone tissue[END_REF][START_REF] Robling | Mechanical signaling for bone modeling and remodeling[END_REF]. Also, recently, ultrasound methods [START_REF] Vallet | Predicting bone strength with ultrasonic guided waves[END_REF][START_REF] Minonzio | Bone cortical thickness and porosity assessment using ultrasonic guided waves: an ex vivo validation study[END_REF] have been proposed to assess in vivo cortical bone porosity, a clinical marker of bone health [START_REF] Bala | Role of cortical bone in bone fragility[END_REF]: in short, ultrasound probes mesoscale elasticity from which porosity can be back-calculated relying on a continuum mechanics homogenization model. Such a model requires extensive validation. The above considerations are strong motivations to investigate the variation of mesoscale elasticity associated with modifications of the pore network, which is the focus of this work.

Investigations of cortical bone elastic properties at the mesoscale are conveniently conducted resorting to a two-phase material model, i.e., a solid matrix (collagen fibers reinforced by mineral particles of hydroxyapatite) pervaded by the vascular porosity, and disregarding any fluid flow [START_REF] Crolet | Compact bone: numerical simulation of mechanical characteristics[END_REF][START_REF] Hellmich | Can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independent phase properties and their interactions?[END_REF][START_REF] Grimal | A two-parameter model of the effective elastic tensor for cortical bone[END_REF][START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF]. In this model, small pores are considered to contribute to the effective properties of the matrix. This two-phase model is particularly well suited to investigate the contribution of the vascular pores to the mechanical properties. As the healthy vascular pore network of a young individual resembles a regular assembly of cylindrical pores [START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -ii. measurements of elastic properties and microhardness[END_REF], several authors have used an idealized representation of the network, namely infinite cylindrical pores distributed randomly or on a hexagonal lattice, and associated analytical homogenization methods such as Mori-Tanaka method [START_REF] Hellmich | Can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independent phase properties and their interactions?[END_REF][START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] or homogenization of periodic media (asymptotic homogenization, AH) [START_REF] Crolet | Compact bone: numerical simulation of mechanical characteristics[END_REF][START_REF] Parnell | The influence of mesoscale porosity on cortical bone anisotropy. investigations via asymptotic homogenization[END_REF]. However, with ageing the structure of the vascular pore network departs from a regular assembly of cylindrical pores, in particular the average pore size increases with age as well as the range of pore size and canal interconnections [START_REF] Cooper | Age-dependent change in the 3d structure of cortical porosity at the human femoral midshaft[END_REF][START_REF] Bousson | Distribution of intracortical porosity in human midfemoral cortex by age and gender[END_REF]. Furthermore, the idealized representation of infinite cylindrical pores disregards the effect of Volkman's canal which are transverse to the bone axis. It is unknown to which extent idealized microstructures are suitable to accurately account for the variations of elastic properties of real cortical bones. This calls for models accounting for the exact structure of the pore network, which can be conveniently captured by X-ray micro-computed tomography (µCT). Few authors have attempted to do such modelling.

Baumann et al. [START_REF] Baumann | The relative influence of apatite crystal orientations and intracortical porosity on the elastic anisotropy of human cortical bone[END_REF] have used a µFE model to investigate the effect of mineral particles orientation and porosity on the elastic anisotropy by analyzing three longitudinal stiffness constants, disregarding shear stiffness constants. Donaldson et al. [START_REF] Donaldson | Relating age and micro-architecture with apparent-level elastic constants: a micro-finite element study of female cortical bone from the anterior femoral midshaft[END_REF] used a µFE model to investigate the correlations between parameters of the bone microstructure (porosity, pore surface to tissue volume, mean pore diameter and degree of geometric anisotropy) specifically focusing on heterogeneity between endosteal and periosteal sites, and evolution with aging. They used small volumes of interest (VOI) leading to a large porosity range (about 0 -50%). They found that porosity could explain almost all the variations of the effective elastic constants and that the contribution of other micro-architecture parameters was small. However, the use of VOIs smaller than the representative volume element [START_REF] Grimal | A determination of the minimum sizes of representative volume elements for the prediction of cortical bone elastic properties[END_REF] in the latter study may obscure the dependence of elasticity on some micro-architecture parameters. Sansalone et al. [START_REF] Sansalone | Anatomical distribution of the degree of mineralization of bone tissue in human femoral neck: Impact on biomechanical properties[END_REF] used synchrotron radiation micro-computed tomography images coupled to a multiscale micromechanical model to account for the variation of local porosity and mineral density. They described the heterogeneity of mesoscale elasticity in the femoral neck. However, they modeled the pores as cylinders, disregarding the details of the micro-architecture.

It is interesting to draw a parallel between the current development of knowledge regarding cortical bone and earlier researches regarding trabecular bone. Trabecular bone is a spongy tissue found in large quantity in the inner part of bones, formed of a network of thin rods and plates of mineralized collagen (trabeculae). Its porosity is typically larger than 80 % in human. The mechanical properties of trabecular bone were first accounted for with idealized geometrical models [START_REF] Gibson | The mechanical behaviour of cancellous bone[END_REF][START_REF] Beaupre | Finite element analysis of a three-dimensional open-celled model for trabecular bone[END_REF][START_REF] Hollister | Application of homogenization theory to the study of trabecular bone mechanics[END_REF]. The morphological complexity of the trabecular network has been accessible ex vivo with non-invasive techniques such as µCT for several decades [START_REF] Feldkamp | The direct examination of three-dimensional bone architecture in vitro by computed tomography[END_REF]. The analysis of the digital images provides parameters such as volume fraction, trabecular thickness, inter-trabecular spacing, trabecular number [START_REF] Hildebrand | Direct threedimensional morphometric analysis of human cancellous bone: microstructural data from spine, femur, iliac crest, and calcaneus[END_REF], or architectural anisotropy using the mean intercept length method and a fabric tensor [START_REF] Harrigan | Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor[END_REF]. Popular analytical methods to predict mechanical properties of trabecular bone consist in implementing formula using the fabric tensor [START_REF] Cowin | The relationship between the elasticity tensor and the fabric tensor[END_REF][START_REF] Zysset | An alternative model for anisotropic elasticity based on fabric tensors[END_REF][START_REF] Zysset | A review of morphology-elasticity relationships in human trabecular bone: theories and experiments[END_REF], or multiscale models in the framework of continuum micromechanics [START_REF] Fritsch | universal' microstructural patterns in cortical and trabecular, extracellular and extravascular bone materials: Micromechanics-based prediction of anisotropic elasticity[END_REF][START_REF] Oftadeh | Hierarchical analysis and multi-scale modelling of rat cortical and trabecular bone[END_REF]. A transition occurred in the 1990's when it became possible to account for the full complexity of the trabecular network architecture with the advent of the so-called micro-finite element modelling (µFEM) [START_REF] Van Rietbergen | A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models[END_REF][START_REF] Van Rietbergen | Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture[END_REF] enabling to solve the homogenization problems using a three-dimensional (3D) numerical model of the real micro-architecture obtained with µCT. Micro-finite element models, which in general assume that bone tissue is isotropic, can predict with a good accuracy the elastic properties measured along different directions of a bone sample [START_REF] Kabel | The role of an effective isotropic tissue modulus in the elastic properties of cancellous bone[END_REF]. µFE modelization is considered as a reference method to estimate elastic properties of trabecular bone samples and of skeletal sites [START_REF] Van Rietbergen | A survey of micro-finite element analysis for clinical assessment of bone strength: The first decade[END_REF]. In comparison, the research on the role of microstructure in cortical bone stiffness is in a relative infancy because only idealized micro-architectures have been considered.

Twenty years after the successful introduction of µFEM to compute trabecular bone effective properties, µFEM or other numerical techniques such as Fast Fourier Transform (FFT) methods [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] can be conveniently used to compute the homogenized mechanical properties of cortical bone. We believe this opportunity will allow rapid progress in the understanding of the effects of cortical bone micro-architecture on different mechanical properties.

The aim of this paper is to quantify the relative contribution of the vascular porosity and its microstructure (i.e., pore number, pore size and pore shape) to bone elasticity.

Considering inter-individual variations of bone matrix properties caused by various factors such as the physical activity, age, nutrition and health conditions among donors, the question can hardly be addressed experimentally. Indeed, bone matrix material properties should be controlled in order to evidence the sole contribution of the micro-architecture to mesoscale elasticity. Consequently we resort to a model of bone which is as a two-phase material with a matrix having fixed 'universal' material properties but various possible micro-architectures corresponding to real bone microstructures imaged with synchrotron radiation micro-computed tomography imaging (SR-µCT) [START_REF] Nuzzo | Quantification of the degree of mineralization of bone in three dimensions using synchrotron radiation microtomography[END_REF] enabling an accurate assessment of micro-architecture descriptors [START_REF] Chappard | 3d characterization of pores in the cortical bone of human femur in the elderly at different locations as determined by synchrotron micro-computed tomography images[END_REF][START_REF] Ostertag | Multiscale and multimodality computed tomography for cortical bone analysis[END_REF]. The full set of anisotropic effective elastic constants predicted by the model are calibrated with experimental data obtained with reso-nant ultrasound spectroscopy (RUS) [START_REF] Bernard | Accurate measurement of cortical bone elasticity tensor with resonant ultrasound spectroscopy[END_REF], a state-of-the-art method to measure the stiffness tensor of cortical bone specimens.

An ancillary objective is to introduce the FFT homogenization technique [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] to model cortical bone stiffness. The method, which has recently been used to study trabecular bone elasticity [START_REF] Colabella | Calculation of cancellous bone elastic properties with the polarization-based fft iterative scheme[END_REF], is commonly used to study the mechanics of engineering materials. FFT-based numerical methods are efficient in term of memory storage and computation cost to solve the heterogeneous local problems using directly images of bone microstructure of relatively large specimens. It has been applied to evaluate the effective properties of composites and porous materials, including elastoplastic responses [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF], permeability [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF], conductivity [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] and piezoelectricity [START_REF] Brenner | Numerical computation of the response of piezoelectric composites using fourier transform[END_REF][START_REF] Brenner | Computational approach for composite materials with coupled constitutive laws[END_REF].

Materials and methods

Specimens

Bone specimens were harvested from the left femur of 29 human cadavers. The femurs were provided by the Départment Universitaire d'Anatomie Rockefeller (Lyon, France) through the French program on voluntary corpse donation to science. The tissue donors or their legal guardians provided informed written consent to give their tissue for investigations, in accord with legal clauses stated in the French Code of Public Health. Among the donors, 16 were females and 13 were males (50 -95 years old, 77.8 ± 11.4, mean ± SD).

The fresh material was frozen and stored at -20 • C.

The samples were slowly thawed and then, for each femur, approximately a 10 mm-thick cross-section was cut perpendicular to the bone axis from the mid-diaphysis. Using a watercooled low speed diamond wire saw (Model 3241, Well, Lyon, France), two rectangular parallelepiped shaped specimens were prepared in the lateral and medial anatomical quadrants of each cross-section (Figure 1). The nominal specimen size was 3×4×5 mm 3 in radial (axis 1), circumferential (axis 2) and axial direction (axis 3), respectively, defined by the anatomical shape of the femoral diaphysis. This specimen size was chosen as a compromise between the limited cortical shell thickness and the need to reduce stiffness measurement errors caused by small sample dimensions. We took care not to include samples presenting large resorption cavities close to the endosteum. Also, this size is believed to be sufficient to include statistically representative micro-architectural information in the µCT images.

The mass density of each specimen which is needed in stiffness determination was derived from the average values of four mass (Sartorius CPA224s, precision: 0.1 mg) and dimensions measurements (Mitutoyo Coolant Proof Caliper 500-606, precision: 0.01 mm). All specimens were kept hydrated during sample preparation. Specimens from one subject with a porosity larger than 30% were not considered as representative of cortical bone and were discarded. One specimen was broken during sample preparation. Finally, 55 specimens from 28 subjects were included in the analysis.

Stiffness measurements with resonant ultrasound spectroscopy

The determination of the stiffness constants of bone specimens by RUS follows the procedure previously described [START_REF] Bernard | Accurate measurement of cortical bone elasticity tensor with resonant ultrasound spectroscopy[END_REF][START_REF] Bernard | Bayesian normal modes identification and estimation of elastic coefficients in resonant ultrasound spectroscopy[END_REF] and consists of the following steps:

(1) A bone specimen was placed on two opposite corners between two ultrasonic transducers (V154RM, Panametrics, Waltham, MA), one for emission and one for reception, to achieve a free boundary condition for vibration (see the experimental setup in Figure 1-b).

(2) The frequency response of the vibration in a specified bandwidth was amplified by a broadband charge amplifier (HQA-15 M-10 T, Femto Messtechnik GmbH, Berlin, Germany)

and then recorded by a vector network analyzer (Bode 100, Omicron Electronics GmbH, Klaus, Austria), from which the 20 to 30 first resonant frequencies were extracted.

(3) Assuming a transversely isotropic symmetry [START_REF] Yoon | Ultrasonic wave propagation in human cortical bone -ii. measurements of elastic properties and microhardness[END_REF][START_REF] Orías | Anatomic variation in the elastic anisotropy of cortical bone tissue in the human femur[END_REF], the stiffness constants C ij (ij = 11, 33, 13, 44, 66) (Voigt notation), were automatically calculated by optimizing the misfit function between the experimental and model predicted resonant frequencies (inverse problem), which is formulated in a Bayesian framework [START_REF] Bernard | Bayesian normal modes identification and estimation of elastic coefficients in resonant ultrasound spectroscopy[END_REF].

The prior information on the distribution of the stiffness constants, required for the Bayesian analysis, was taken from a previous study [START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF] The quality of the experimental data and of the determined C ij are reflected in the misfit error (RMSE) σ f between the experimental and predicted resonant frequencies. We found for all the specimens σ f < 0.7% (mean value = 0.4% ) which is small in RUS measurements of cortical bone [START_REF] Bernard | Elasticity-density and viscoelasticity-density relationships at the tibia mid-diaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF]. This indicates that the assumption of transverse isotropy is very reasonable. Moreover, the differences between the effective stiffness constants in the radial and circumferential directions evaluated with homogenization using the SR-µCT images of bone microstructure were small (see Section 3.2). This means that the orthotropic behavior created by the micro-architecture is small and it supports the assumption of transverse isotropy.

Bone microstructure

After RUS measurements, bone specimens were defatted for 12 hours in a chemical bath of diethylether and methanol (1:1) and rinsed in distilled water before SR-µCT scanning in order to comply with the local regulation by the European Synchrotron Radiation Facility (ESRF, Grenoble, France). Details of the defatting protocol can be found in Cai et al. [START_REF] Cai | Cortical bone elasticity measured by resonant ultrasound spectroscopy is not altered by defatting and synchrotron x-ray imaging[END_REF].

Then, the specimens were scanned using SR-µCT 3D imaging, which was performed on the beamline ID19 at ESRF. SR-µCT is the gold standard to assess cortical bone microstructure due to the high flux, the monochromaticity and the parallel beam geometry [START_REF] Ostertag | Multiscale and multimodality computed tomography for cortical bone analysis[END_REF].

The SR-µCT setup is based on a 3D parallel beam geometry acquisition [START_REF] Salomé | A synchrotron radiation microtomography system for the analysis of trabecular bone samples[END_REF][START_REF] Weitkamp | Status and evolution of the esrf beamline id19[END_REF]. The beam energy was tuned to 26 keV by using a (Si111) double crystal monochromator. A full set of 2D radiographic images were recorded using a CDD detector (Gadox scintillator, optic lenses, 2048 × 2048 Frelon Camera) by rotating the specimen in 1999 steps within a 360 • range of rotation. The detector system was fixed to get a pixel size of 6.5 µm in the recorded images in which a region of interest of 1400×940 pixels was selected to fit the specimen.

For each specimen, the SR-µCT images were reconstructed to obtain the 3D volume of the specimens.

The orientation of the 3D volume of each specimen was corrected by slightly rotating the ensemble of voxel using Fiji [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF] so that the image reference frame coincides with the orientation of the specimen's faces. As the shape quality of our specimens was quite good (perpendicularity and parallelism errors were about 1 • , see [START_REF] Cai | Quantification of stiffness measurement errors in resonant ultrasound spectroscopy of human cortical bone[END_REF]), the potential misalignments due to orientation correction should be small. In each specimen, a rectangular parallelepiped volume of interest (VOI, sized ∼ 2.8 × 3.9 × 4.8 mm 3 ) was selected manually on the principle of retaining the maximum volume for morphometric analysis, as well as for the evaluation of the effective stiffness using FFT homogenization (see Section 2.6). Note that VOIs are much larger than the recommended representative volume element size [START_REF] Grimal | A determination of the minimum sizes of representative volume elements for the prediction of cortical bone elastic properties[END_REF]. Each VOI was chosen to be comparable to the specimen's volume in order to allow direct comparison between measured and calculated stiffness. Following Bala et al. [START_REF] Bala | Pore network microarchitecture influences human cortical bone elasticity during growth and aging[END_REF], the VOIs for the morphometric analysis were binarized treating the void volumes as the solid and the bone phase as a background (Figure 1 (d)). Then, the following microstrutural variables (Table 1)

were calculated using the software CTAnalyser (V 1.16.1, Skyscan NV, Kontich, Belgium): pore volume fraction (φ), pore surface to pore volume ratio (PoS/PoV), the average diameter of the pores (PoDm), the average separation between pores (PoSp), pore number (PoN), connectivity density (ConnD), pore pattern factor (PoPf), structure model index (SMI) and degree of anisotropy (DA). In the process of the DA calculation, the fabric tensor of the porous network measured by mean intercept length (MIL) analysis [START_REF] Harrigan | Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor[END_REF][START_REF] Odgaard | Three-dimensional methods for quantification of cancellous bone architecture[END_REF] was also obtained, from which the main orientation of the Haversian canals can be evaluated. The angle (α) of the misalignment of the Haversian canals relative to the longitudinal direction (material basis) was also extracted. A summary of the microstrutural variables and their definitions are listed in Table 1.

To study the effect of pore shape on the elastic symmetry in the transverse plane, we defined a transverse plane geometrical factor R d calculated as follows. Pores were regarded as ellipsoids and the average ratio between the two short (circumferential and radial) axis lengths of the pores were calculated, in analogy to the MIL analysis for fabric tensor. Specifically, in each two-dimensional image of the CT dataset (plane defined by axes 1 and 2), the axis length of the pores in one direction was estimated as the total number of pixels contained inside the pores divided by the number of pores along each line of the image.

Then, the ratio between the average axis length estimated on each line in the circumferential (axis 2) and radial (axis 1) directions was calculated. Finally, R d was calculated as the average ratio over the entire volume. Note that R d does not quantify the aspect ratio of Haversian canals in the transverse plane but rather a combination of the aspect ratios of Haversian and Volkmann's canals in that plane. The aforementioned steps to calculate R d was implemented in a customized Matlab script. ) were evaluated with the FFT method in the framework of linear elasticity using the VOI as the unit cell for homogenization. After binarization, one VOI contains bone matrix voxels and pore voxels. They were respectively allocated elastic properties C m , a transverse isotropic stiffness tensor to be defined as described in Section 2.5, and C p , corresponding to water, with bulk modulus = 2.2 GPa and a null shear modulus. Note that inter-and intra-specimen mineralization heterogeneity was disregarded after binarization, hence the same matrix elasticity was attributed to every matrix voxel for every specimen.

The local mechanical problem in the VOI (V ) consists of the equilibrium equation, generalized Hooke's law as constitutive equation, compatibility of the displacement field and boundary conditions. The phases (bone matrix and the material filling the pores) are assumed to be perfectly bonded. The local problem closed by periodic boundary conditions [START_REF] Suquet | Homogenization Techniques for Composite Media[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF] can be solved by the FFT-based numerical approach proposed by Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. More details of the method can be found in Appendix A.

Developments of the FFT-based method have been proposed by several authors, in particular to improve its convergence in the case of a high mechanical contrast on the local properties [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF][START_REF] Brisard | Fft-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. In the present study, we have used the augmented Lagrangian scheme originally proposed by Michel et al. [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF] and later reinterpreted by Moulinec and Silva [START_REF] Moulinec | Comparison of three accelerated fft-based schemes for computing the mechanical response of composite materials[END_REF] as a special case of the polarization-based scheme of Monchiet and Bonnet [START_REF] Monchiet | A polarization-based fft iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. The latter is the one which has been chosen in the recent study on trabecular bone elasticity [START_REF] Colabella | Calculation of cancellous bone elastic properties with the polarization-based fft iterative scheme[END_REF].

Once the local stress and strain fields are calculated, the effective stiffness tensor C F F T is obtained with the relation between average stress and strain tensors

σ = C F F T : ε . ( 1
)
The computation of the twenty-one independent coefficients of C F F T is classically done by considering six independent macroscopic strain loadings.

Preliminary tests (Appendix B.2) showed that increasing the size of the voxels to 35 µm leads to an acceptable error on the solution. The total number of voxel in each VOI was about 1.2 million. The 3D images of VOIs are provided in the online supplementary material.

One loading in the FFT homogenization method took around 50 iterations in our specimens and the effective stiffness tensor could be computed in about 1 minute on a workstation (CPU, Intel Xeon E5-2695 v3, 8 threads per loading).

Our implementation of the FFT method was validated by comparison with the analytic solution of the homogenized properties of hexagonal periodic pavement as reported by Bravo-Castillero et al. [START_REF] Bravo-Castillero | Unified analytical formulae for the effective properties of periodic fibrous composites[END_REF] (see Appendix B.1). After the calculation of the 21 terms in C F F T , we set to 0 the terms which are zero for orthotropy. We assume negligible deviation of the material orthotropy axes compared to the sample axes. The orthotropy error due to this assumption was evaluated following van Rietbergen et al. [START_REF] Van Rietbergen | Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture[END_REF] (see equations 6 and 7 in [START_REF] Van Rietbergen | Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture[END_REF]).

Model with idealized microstructure

For comparison purposes, we also calculated sample-specific effective properties using the pore volume fraction φ as a parameter and assuming that the pore network is an ensemble of infinite cylindrical pores periodically distributed on a hexagonal lattice (Figure 1 (e)).

The material in pores and the bone matrix are allocated the stiffness C p as defined above and C m as defined in Section. 2.6. The effective properties C cyl ij were calculated using the AH method [START_REF] Grimal | A two-parameter model of the effective elastic tensor for cortical bone[END_REF][START_REF] Parnell | The influence of mesoscale porosity on cortical bone anisotropy. investigations via asymptotic homogenization[END_REF].

Calibration of bone matrix stiffness

For all the specimens' models, we used a unique bone matrix stiffness tensor C m chosen such that the effective tensors C F F T ij for the different specimens of the set match experimental data C EXP ij . The elastic symmetry of C m was assumed to be transverse isotropy with the axis of transverse isotropy coincide with the longitudinal axis (axis 3). This calibrated bone matrix stiffness C m was accordingly obtained by minimizing the objective function defined as

F (C m ) = 55 k=1 ij C EXP ij;k -C F F T ij;k (C m ; C p ) C EXP ij;k 2 , (2)
where C EXP ij;k and C F F T ij;k refer to the experimental and effective stiffness constants of the kth specimen.

Data analysis

All analysis on modeled stiffness C F F T ij;k were conducted after the calibration of the numerical model with experimental data.

A possible anisotropy in the transverse plane (defined by axes 1 and 2) due to the shape of the pores was investigated using R d and the relative difference between C F F T To investigate the relationships between each stiffness constant C F F T ij and the microstructural variables, linear correlation analysis and stepwise multiple linear regression analysis were conducted. We assessed the extent to which variations of bone stiffness can be accounted for by variations of porosity only. The stiffness constants C F F T ij and C EXP ij were fitted to an exponential function of φ because the variation of porosity were found to explain a major amount of the variations of bone stiffness [START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF]. The exponential fit was found to be superior to the linear fit (see Figure 2),

C f it ij = a ij e bij φ , (3)
where a ij and b ij are the coefficients to be optimized. The root-mean-square-error (RMSE) between the fit and data was calculated as,

RM SE = 1 55 55 k=1 C f it ij;k -C ij C ij 2 × 100%, (4)
where

C ij stands for either C F F T ij or C EXP ij
. RMSE values were used to reveal the part of the variations of stiffness that cannot be explained by the variations of porosity.

The model with the realistic micro-architecture (C F F T ij ) and that with the idealized micro-architecture (C cyl ij ) were finally compared to reveal a potential bias in predicted stiffness when details of the micro-architecture are disregarded. The significance level was set to p < 0.05 for all statistical analyses. As illustrated by the RMSE with the experimental fits (equation ( 3)), the fluctuations of C EXP ij (RMSE between 3.1 -6.4%) at a given porosity value are larger than those of C F F T ij (RMSE between 0.6 -1.7%). This is expected because of the experimental uncertainty to measure stiffness and because of possible differences of bone matrix stiffness between specimens. Note that bigger discrepancies between C F F T 13 and C EXP 13 was observed especially when φ > 12% which was not fully understood and requires further investigation. A reason that could partially explain this observation could be that RUS measurement errors are larger for off-diagonal stiffness constants compared to shear and longitudinal stiffness constants [START_REF] Cai | Quantification of stiffness measurement errors in resonant ultrasound spectroscopy of human cortical bone[END_REF][START_REF] Zadler | Resonant ultrasound spectroscopy: theory and application[END_REF].

Results

Calibration of the model

Pore shape and the anisotropy in the transverse plane

The range of R d was between 0.79 and 1.21 (Mean±SD, 1.02±0.08) and slight anisotropy between C F F T was observed as the δ 21 was between -5.0% and 6.5% (Mean±SD, 0.5% ± 1.6%) and δ 45 was between -5.5% and 7.3% (Mean±SD, 0.8% ± 1.7%). Pore shape parameter R d significantly correlates with δ 21 (r = 0.77, p < 10 -11 ) and δ 45 (r = 0.73, p < 10 -9 ) (Figure 3) suggesting that pore shape contributes to the orthotropic elastic behaviour of cortical bone. However, this anisotropy effect is small, as the relative differences δ 21 or δ 45 of only 5 out of 55 specimens exceed 3% (between 3.4% and 7.3%). In order to compare transverse isotropic stiffness (experimental data C EXP ij and results from the cylinder model C cyl ij ) to effective stiffness constants C F F T ij obtained in an orthotropic framework, we retained the mean value of (C F F T

11

, C F F T

22

) and the mean value of (C F F T

44

, C F F T

55

).

Descriptive statistics

The orthotropic errors were between 0.6 -15.1% (Mean ± SD = 4.7 ± 3.5%). The values of C F F T ij and the microstructural variables are summarized in Table 2, in which one can observe that the coefficients of variation (CV) of the microstructural variables (7.8 -51.1%)

were wider compared to that of C F F T ij (6.9 -10.8%). The mean value of the misalignment angle (α) was small (6 • ). Nevertheless, the relatively large (16.4 • ) misalignment angle of one specimen may be caused by the misalignment between cutting direction and bone axis during 4.

sample preparation. The orthotropic error of this specimen was 8.4%. The ratio between the maximum and second biggest eigenvalues (corresponding to two transverse directions) was 1.2±0.1 (Mean±SD) much smaller than the values of DA (4.2 ± 0.7, Table 2), indicating that the microstrucure has a quasi-isotropic pattern in the radial-circumferential plane.

Correlations between bone stiffness and microstructural variables

C F F T 11 (GPa) C F F T 33 (GPa) C F F T 13 (GPa) C F F T 44 (GPa) C F F
C F F T ij
remained significant (-0.57 < adj-r < 0.65), but lower than the non-adjusted ones.

The results of stepwise multiple linear regression analysis further showed that the adj-r 2 values with φ in the model alone were between 98.9 -99.7% for all the C F F T ij . Adding the other microstructural variables in the regression models, the adj-r 2 values were improved by less than 1.0%.

Significant but weak correlations between DA and the anisotropy ratios of the stiffness constants were observed (Figure 4) (r = 0.48 for C F F T

33

/C F F T 11 and r = 0.36 for

C F F T 44 /C F F T

66

). (in GPa) and microstructural variables and their partial correlation coefficients (adj-r) after controlling for the contribution of φ. -0.45

C F F T 11 C F F T 33 C F F T 13 C F F T 44 C F F T 66 φ (%) r -0.
2 3 4 5 6 1.2 1.3 1.4 1.5 1.6 1.7
Figure 4: Relationship between DA and the anisotropy ratios of the stiffness constants. r is the Pearson's correlation coefficient.

Modeling of stiffness accounting only for variations of porosity

The RMSEs (equation ( 4)

) between C EXP ij , C F F T ij
and their corresponding exponential models (equation ( 3)) were between 3.1-6.4% and 0.6-1.7%, respectively. As expected, the in Table 5.

Discussion

In this work, we have investigated the relationships between the micro-architecture of the cortical bone vascular pore network and the effective stiffness at the mesoscale. We used the elastic data (transverse isotropic stiffness tensor measured with RUS) and microstructures (obtained with SR-µCT) of a collection of 55 human cortical bone specimens from the mid-diaphysis of the femur. Specimen-specific numerical models were created using each specimen's microstructure and the bone matrix stiffness, which is identical for all the specimen's models, was determined in a calibration procedure to minimize the discrepancies between the experimental and numerical data sets.

Doing so, we have built a numerical model that accounts reasonably well for the observed variations of experimental stiffness (Figure 2) in which the elastic properties of the bone matrix are controlled. This last point is crucial to analyze the minute influence of the details of pore micro-architecture on stiffness.

The small orthotropy error in C F F T indicates that orthotropy assumption is reasonable within the sample. The mineral density and mechanical properties of the osteonal tissue are heterogeneous in particular because of the presence of recently remodeled tissues [START_REF] Raum | Site-matched assessment of structural and tissue properties of cortical bone using scan-29 ning acoustic microscopy and synchrotron radiation µct[END_REF]. The effect of modeling the heterogeneity of elastic properties within trabeculae, associated to variable mineral contents, has been investigated. It was found that modeling the tissue as homogeneous lead to a slightly overestimated trabecular bone effective stiffness compared to the stiffness obtained with a model incorporating tissue heterogeneity [START_REF] Gross | Mineral heterogeneity has a minor influence on the apparent elastic properties of human cancellous bone: a srµct-based finite element study[END_REF][START_REF] Kaynia | Effect of intraspecimen spatial variation in tissue mineral density on the apparent stiffness of trabecular bone[END_REF]. However, the effect of including matrix heterogeneity in the calculation of effective properties of cortical bone has not been studied as far as we know.

The overall symmetry of the effective stiffness tensor results from the elastic anisotropy of the constituents as well as their shape and spatial distribution. Real pores in cortical bone usually do not have perfect circular cross-sections, especially for the resorption cavities. In this work, we quantified this deviation from a circular cross-section by the average ratio R d between pore axis length in the circumferential and radial directions. We must remind the readers that R d does not define the shape of the Haversian canals, but it corresponds to both the information about the axis length and orientation of Haversian and Volkmann's canals.

The correlations found between R d and δ 21 and δ 45 (Figure 3) support the hypothesis that non-circular pore shapes may cause the anisotropy of the elastic properties in the transverse plane observed here, leading to mesoscopic orthotropy. However, the magnitude of the differences δ 21 and δ 45 were averaged due to the randomized distribution of the pores evidenced by the observation that δ 21 or δ 45 of only 5 specimens exceed 3%. On the other hand, this observation highlights the fact that the FFT homogenization is able to reveal such a small anisotropy in the transverse plane caused by bone microstructural features. Another source of orthotropy could be related to porosity gradient. It is well known that a gradient of porosity exists across cortical bone thickness both in males and females. Specifically, there are more pores and pore size are bigger in the endosteal region than in the periosteal region [START_REF] Bousson | Distribution of intracortical porosity in human midfemoral cortex by age and gender[END_REF][START_REF] Thomas | Regional variation of intracortical porosity in the midshaft of the human femur: age and sex differences[END_REF]. Our specimens were retrieved in a typical periosteal zone with a low porosity and mainly composed of the Haversian system from the modeling process. In a similar study in femoral diaphyses imaged at ID 19 ESRF with a voxel size of 7.5 µm, we have found very close results of porosity about 7.3 ± 6.7% [START_REF] Chappard | 3d characterization of pores in the cortical bone of human femur in the elderly at different locations as determined by synchrotron micro-computed tomography images[END_REF]. However, although we observed a moderate porosity gradient in some of our specimens (data not shown), we did not find significant correlations between δ 21 , δ 45 and the porosity gradient. As the specimens were prepared with caution by avoiding the endosteal region, we assume the gradient was too small to affect the elastic symmetry in our study.

The quantitative relationships between the microstructural characteristics of the porous network and the elastic properties of cortical bone have not been thoroughly investigated, so far. The only study that compares to ours is that of Bala et al. [START_REF] Bala | Pore network microarchitecture influences human cortical bone elasticity during growth and aging[END_REF]. However, a different skeletal site (fibular) was investigated (for children and adults) and different methods were used to derive bone microstructure and stiffness. Also, the potentially confounding factor (matrix stiffness) was controlled in our work, whereas this was not the case in their purely experimental study. Nevertheless, consistent relationships between bone microstructural variables and stiffness constants were observed, e.g., negative correlations between φ, PoN, PoDm and bone stiffness, positive correlations between PoS/PoV, PoSp, PoPf and bone stiffness. Both works showed that significant correlations still exist between several microstructural variables and stiffness constants (e.g., PoSp) after controlling for the contribution of φ. In addition, we observed negative correlations between DA estimated from the fabric tensor and all the effective stiffness constants and the correlations between DA and the elastic anisotropy ratios were positive. Multiple regression analyses further confirmed that almost all the variations of the effective stiffness constants could be explained by porosity which is consistent with the observations of Donaldson et al. [START_REF] Donaldson | Relating age and micro-architecture with apparent-level elastic constants: a micro-finite element study of female cortical bone from the anterior femoral midshaft[END_REF].

As mentioned above, fixing bone matrix stiffness in the FFT homogenization, this study has demonstrated the contribution of details of bone micro-architecture to mesoscale stiffness. However, most of the correlations between microstructural variables and stiffness constants vanished or were weak after removing the effect of porosity (Table 3). This can be explained by the interdependence between microstructural parameters and porosity, which means that the micro-architecture of all the specimens closely follows a specific pattern depending on porosity. It is noteworthy that the RMSE for the exponential fit of C F F T ij (0.6 -1.7%) (Figure 2) is smaller than the estimated experimental errors (in the range 1.7 -3.1% for the different elastic constants) which suggests that the small effects on stiffness caused by the details of the micro-architecture, independent of porosity, cannot be observed in the experimental data due to the limited accuracy of experimental techniques.

The above observation indicates that porosity variations explain a major part of stiffness variations, which has already been observed in several datasets [START_REF] Currey | The effects of strain rate, reconstruction and mineral content on some mechanical properties of bovine bone[END_REF][START_REF] Dong | The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity[END_REF][START_REF] Granke | Change in porosity is the major determinant of the variation of cortical bone elasticity at the millimeter scale in aged women[END_REF]. A consequence is that it may not be necessary to record details of the microstructure (other than porosity), to predict bone stiffness. Such prediction can be made on the basis of empirical models, such as the exponential fit discussed above or a model physically rooted in continuum mechanics and based on an assumption of pores architecture.

The model with vascular pores idealized as infinite cylinders aligned with the bone axis is a popular model in bone biomechanics [START_REF] Crolet | Compact bone: numerical simulation of mechanical characteristics[END_REF][START_REF] Hellmich | Can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissue-independent phase properties and their interactions?[END_REF][START_REF] Grimal | A two-parameter model of the effective elastic tensor for cortical bone[END_REF][START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF]. For the first time, as far as we know, we present data to critically test this model. The idealized stiffness C cyl ij was consistently larger than C F F T ij for all coefficients (1.0 -14.7%) in the range of the investigated porosity (Figure 5). The difference becomes measurable if it is greater than experimental error and the significance of such a difference depends on the acceptable error for one application.

Previous calculations in SiC/SiC composite materials [START_REF] Chateau | Micromechanical modeling of the elastic behavior of unidirectional cvi sic/sic composites[END_REF] suggest that the decrease of C F F T 11 relative to C cyl 11 may be explained by non-cylindrical shape of the pores and random spatial distribution in real microstructures. Indeed, the authors have shown that the non-cylindrical shape results in a strong decrease of the Young's modulus in the transverse direction whereas no effect of both random distribution and irregular shape of the pores was observed on the Young's modulus in the axial direction (see also [START_REF] Aghdam | Micromechanics based analysis of randomly distributed fiber reinforced composites using simplified unit cell model[END_REF]). The differences between AH and FFT for axial moduli thus have to be linked to other microstructural discrepancies than the spatial distribution and the cross-section shape of cylindrical pores. An important microstructural feature of cortical bone, as compared to a parallel cylindrical microstructure, is the presence of transverse connections (Volkmann's canals) between quasi-cylindrical pores (Figure 6).

It can be thus be conjectured that they are responsible for the decrease in the effective axial elasticity (C 33 and C 44 ). For the transverse shear stiffness constants C 66 , we noticed that the results from AH are very close to the upper Hashin-Rosen bound [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF] which has been depicted in Parnell et al. [START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF]. Based on the aforementioned observations, we think it is important to point out that the idealized model assuming infinite long cylindrical pores yields, to some extent, biased effective stiffness. An ancillary purpose of this work is to show that FFT homogenization is an attractive alternative to µFEM which is popular in the bone micromechanics community to compute local mechanical fields and effective properties. The meshing free feature and the low numerical cost (problems with several millions of degrees of freedom can be solved in a few minutes) pave the way for the investigations involving large series of specimens and specimens with large size.

Limitations in this work are that (i) the specimens came from only a group in the elderly with unknown healthy status and (ii) only one site (mid-diaphysis) of the femoral bone was investigated. Therefore, the observations and conclusion should be restricted to the specimens used in this work. In a future work, the findings in this work should be tested using specimens from different ages (young to adults), unhealthy bones and various skeletal sites. It may be also worthy to relate bone microstructure and anisotropic elasticity using fabric tensor based models, an approach widely used in trabecular bone [START_REF] Zysset | A review of morphology-elasticity relationships in human trabecular bone: theories and experiments[END_REF]. A comparison with fabric tensor based models which provide analytical solutions may provide different insights on the morphology-elasticity relationship for cortical bone.

The present study has been limited to moderate loading intensities for which the local mechanical response remains elastic. However, the FFT numerical method can as well be used to study the nonlinear behavior after yield. Besides, the relationship between bone microstructure and local mechanical behavior such as local strain which regulates bone remodelling process may be interesting to investigate.

Conclusion

Providing the evidence from the multiscale experimental data in this work, as well as the results from the two homogenization approaches, we found that, among the microstructural variables, porosity explains most of the variations of bone elasticity, all tissue properties being equal otherwise. The effects caused by other microstructural features are smaller than the experimental errors and therefore are masked by experimental uncertainties. Hence, we propose that measuring detailed microstructure for cortical bone stiffness prediction is unnecessary, i.e., knowing porosity is enough. However, a proper biomechanical model, yet to be established, should be chosen as idealizing the pores as infinite cylinders may yield biased stiffness values compared to stiffness predicted from the real microstructure. The conclusion should hold true at least for human femoral mid-diaphysis in the elderly.

FFTFigure 1 :

 1 Figure 1: Sample preparation procedure and the experiments, (a) a cross-section of femoral bone from which two specimens were extracted at the lateral and medial quadrants. (b) RUS measurements for bone stiffness. (c) A representation of the 3D microstructure of a bone specimen imaged by SR-µCT adapted from [11]. (d) Pore network obtained treating bone phase as background. (e) Bone microstructure idealized as infinite cylindrical pores periodically distributed in bone matrix on a hexagonal lattice (see Section 2.5).

  α • misalignment angle, the angle between the eigenvector of minimum eigenvalue and the vector of the longitudinal direction (material basis) R d a.u. transverse plane geometrical factor, average aspect ratio of pores in the transverse plane 2.4. Numerical model and method of solution Each VOI was considered as a representative volume element of bone material. The homogenized effective stiffness constants (C F F T ij

11 and C F F T 22 , 22 - 55 ×

 11222255 and between C F F T 44 and C F F T 55 , denoted as δ 21 = C F F T 100%, respectively.

  After minimizing the objective function between C EXP ij and C F F T ij (equation (2)), bone matrix stiffness C m ij was found to be C m 11 = 24.4 GPa, C m 33 = 33.5 GPa, C m 13 = 14.7 GPa, C m 44 = 6.9 GPa, C m 66 = 5.4 GPa. The constants in C EXP ij agree well with C F F T ij except for C 13 when φ > 12% (Figure 2).

11 and C F F T 22 ,

 1122 and between C F F T 44 and C F F T 55

Figure 2 :

 2 Figure 2: Calibration and validation of the model. Stiffness tensors C EXP ij

Pearson's correlationFigure 3 :C F F T 22 and C F F T 11 , and between C F F T 44 and C F F T 55 .

 322114455 Figure 3: Relationship between the transverse plane geometrical factor R d and the stiffness computed from the microstructure of the specimens. The vertical axis gives the relative difference (δ 21 and δ 45 ) between C F F T 22 and C F F T 11 , and between C F F T 44 and C F F T 55 . r is the Pearson's correlation coefficient.

  RMSEs for C F F T ij were smaller than for C EXP ij and even smaller than experimental errors (1.7 -3.1%). The coefficients of the exponential fit (equation (3)) for C F F T ij and C EXP ij are summarized in Table 4. The comparison between C cyl ij and C F F T ij plotted against porosity (Figure 5) shows that all the constants in C cyl ij were consistently larger than C F F T ij which considers real bone microstructure. The discrepancy increases as porosity increases. The discrepancy between each constant in C cyl ij and C F F T ij in the range of the investigated porosity are summarized

[ 39 ]

 39 . The effective stiffness constants C F F T ij estimated by FFT homogenization compared well with the experimental ones C EXP ij (Figure2). As expected, the values of C F F T ij exhibited a smaller dispersion (Figure2) compared to experimental values. This is related to the facts that, (i) the homogenization model did not take into account the variations of bone matrix stiffness between individuals and (ii) the experimental data includes random measurement errors. Moreover, the discrepancies between the RMSE for C F F T ij and C EXP ij cannot be entirely explained by experimental uncertainties. Overall, the above observation suggests that bone matrix stiffness variations have a measurable contribution to the variations of bone mesoscopic stiffness. In our model, the matrix was assumed to be homogeneous

Figure 5 :

 5 Figure 5: Comparison between C F F T ij and C cyl ij to reveal the effect of the simplified bone microstruture (Figure 1-e) on the effective bone stiffness. Stiffness data was normalized by dividing the value of the corresponding matrix stiffness coefficient (C m ij ).

Figure 6 :

 6 Figure 6: Quasi-cylindrical pores (Haversian canals) connected by the transversely oriented Volkmann's canals.

Figure B. 1 :

 1 Figure B.1: A hexagonal cell of 15% porosity and the evolution of the relative differences of stiffness constants using different values of η.

Figure B. 2 :

 2 Figure B.2: The evolution of the relative differences of stiffness constants (C 11 and C 44 ) of four specimens (5.0 -21.2% porosity) varying pixel size of the images.

  achieved on human femoral cortical bone as well. In the stiffness tensor, C 12 = C 11 -2C 66 and (1 -2) is the isotropy plane; C 11 and C 33 which correspond to pure longitudinal waves are denoted as longitudinal stiffness constants, C 44 and C 66 which correspond to pure shear waves are denoted as shear stiffness constants[START_REF] Auld | Acoustic fields and waves in solids[END_REF], and C 12 and C 13 are the off-diagonal stiffness constants. The experimental errors, e.g., irregularity of specimen geometry and uncertainties of the extracted resonant frequencies, following this protocol, typically cause an error of approximately 1.7% for the shear stiffness constants and approximately 3.1% for the longitudinal and off-diagonal stiffness constants[START_REF] Cai | Quantification of stiffness measurement errors in resonant ultrasound spectroscopy of human cortical bone[END_REF]. Bone stiffness constants of the specimens measured by RUS are denoted

	as C EXP ij	in the following text.

Table 1 :

 1 Microstructural variables assessed by CTAnalyser and their definitions.

	Variable	Unit	Definition
	φ	%	pore volume fraction
	PoS/PoV mm -1	pore surface to pore volume ratio
	PoN	mm -1	pore number per mm
	PoDm	µm	average diameter of the pores
	PoSp	µm	average separation between pores
	PoPf	mm -1	pore pattern factor, lower PoPf indicates
			higher concavity, i.e., better-connected pore
			network
	ConnD	mm -3	connectivity density, a measure of the degree
			to which a pore is multiply connected
	SMI	a.u.	structure model index, the relative prevalence
			of rods and plates in a 3D pore network
	DA	a.u.	degree of anisotropy, the ratio between the
			biggest and smallest eigenvalue of the fab-
			ric tensor measured by mean intercept length
			analysis [32, 62]

Table 2 :

 2 A summary (Mean±SD, range) of the values of C F F T

ij and microstructural variables.

Table 3 :

 3 Pearson's correlation coefficients r between C F F T ij

Table 4 :

 4 Coefficients of the exponential fit (see Figure 2 also) for C F F T 11 , b 11 ] [a 33 , b 33 ] [a 13 , b 13 ] [a 44 , b 44 ] [a 66 , b 66 ] EXP [23.28, -0.022] [32.64, -0.016] [13.33, -0.015] [6.98, -0.025] [5.55, -0.033]

	ij	and C EXP ij	.

Table 5 :

 5 The discrepancy between C cyl ij and C F F T ij in the range of the investigated porosity.

		C 11	C 33	C 13	C 44	C 66
	Min	1.3%	1.8%	1.5% 1.0%	1.2%
	Max 9.9% 12.1% 8.0% 5.9% 14.7%

Acknowledgment

This work has received financial support from the Agency National Research under the ANR-13-BS09-0006 MULTIPS project.

Appendix A. FFT homogenization

The unit-cell problem reads: find the displacement field u and the stress field σ such

with periodic boundary conditions. The displacement field u can thus be split into an affine part E•x and a correction term u * (x)

where E is the overall strain, u * is a periodic field characterizing the fluctuation of the displacement due to the presence of the heterogeneities. Consequently, the average of its gradient vanishes and the overall strain E coincides with the average strain in the unit-cell ε .

A usual practice in micromechanics is the introduction of a reference homogeneous medium C (0) to transform an elastic heterogeneous problem into an homogeneous elastic problem with fictitious body forces. The solution strain field may thus be written as

with Γ (0) the strain Green operator corresponding to the homogeneous medium with elasticity C (0) . This is an implicit integral equation for the strain field ε which is obtained as a series expansion reading

Based on this expansion, Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] have proposed a fixed-point algorithm to solve the local problem (Equation (A.1)).

Appendix B. Convergence study

Appendix B.1. Convergence criterion in FFT homogenization

A convergence criterion on the compatibility of the strain field and the fulfillment of the pointwise constitutive relation has been adopted [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF]. The iterative resolution stops when the relative error is less than a chosen value η. method [START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF][START_REF] Parnell | The influence of mesoscale porosity on cortical bone anisotropy. investigations via asymptotic homogenization[END_REF]. The stiffness of the bone and pore phases from Section 2.6 were used. The relative differences of the stiffness constants were calculated as

As shown in Figure B.1, when η = 0.01, the relative differences on all the stiffness constants are less than 0.3%.

Appendix B.2. Choice of the pixel size

The convergence study was carried out on four specimens of 5%, 9.9%, 15.4% and 21.2% porosity, respectively. The original pixel size (ps) of the 3D images was 6.5 µm. The pixel size of the images of each specimen was increased to ps = 15, 25, until 115 µm, respectively, with 10 µm interval. The stiffness constants of the specimens at each ps were calculated by the FFT homogenization and compared with the results when ps = 6.5 µm to estimate the relative differences (equation (B.1)). When ps = 35 µm, the relative differences of all the stiffness constants from all the specimens (from low to high porosity) were less than 1.4%.


K Haddadi
email: kamel.haddadi@univ-lille.fr

E Okada #

K Daffé #

F Mubarak $

D Théron #

G Dambrine

Multiport Vector Network Analyzer Configured in RF Interferometric Mode for Reference Impedance Renormalization

Keywords: multiport vector network analyzer, microwave interferometry, extreme impedance, vector calibration

A novel active microwave interferometric technique is implemented on a multiport vector network analyzer for renormalizing the reference impedance 50 Ohms into any desired complex impedance. The resulting measured reflection coefficient around the new reference impedance is around zero, resulting in high measurement sensitivity. The method proposed avoids any external component commonly found in interferometric set-ups. In addition, a zeroing process including vector calibration is developed for broad frequency range and requires only a software procedure to be implemented in the system framework.

I. INTRODUCTION

Vector network analyzers (VNA) have been introduced in the sixties to measure the electrical properties of microwave devices and circuits [START_REF]An advanced new network analyzer for sweep measuring amplitude and phase from 0,1 to 12.4 GHz[END_REF]- [START_REF] Rytting | ARFTG 50 year network analyzer history[END_REF]. Since then, there has been considerable work to extend their measurement capabilities such as frequency coverage from low RF to hundreds of GHz, multi-port devices characterization, balanced, differential, harmonic, modulation and pulsed measurements. In addition to these advances in the field of vector network analysis, there is an urgent need to further extend the measurement capabilities to meet the specific needs of characterization tools for impedances clearly higher than the 50 Ω reference impedance of the VNA [START_REF] Happy | Measurement techniques for RF nanoelectronic devices : new equipment to overcome the problems of impedance and scale mismatch[END_REF]. In particular, the electrical characterization of high impedance nanodevices in the microwave regime is still challenging [START_REF] Rice | Broadband electrical characterization of multiwalled carbon nanotubes and contacts[END_REF]- [START_REF] Imtiaz | Near-Field Scanning Microwave Microscopy: An Emerging Research Tool for Nanoscale Metrology[END_REF].

Among the solutions, the interferometric principle commonly used in optics through Michelson or Mach Zender configurations has gained the interest from the microwave research community to measure small complex impedance contrasts. Basically, a monochromatic signal is split with one part traveling a fixed path and the other a variable one. The two resulting signals are reunited, amplified if necessary, and produce the well-known interference. In this effort, this paper proposes for the first time a complete interferometric solution in a broad frequency range that requires only a multiport VNA. In other words, the interferometric process is built inside the VNA itself and requires only software resource for its configuration. The method can be implemented easily in any multiport VNA. The principle of operation and methodology are described in Section II. In Section III, the interferometric configuration including calibration is proposed and validated experimentally.

II. BRIEF STATE OF THE ART

To situate the contribution of the proposed work, main related achievements are described in the following and practical limitations are highlighted. In 1949, a Michelson-type microwave interferometer is introduced for the first time at the free-space wavelength of 3.2 cm for free-space dielectric characterization [START_REF] Lengyel | A Michelson-Type Interferometer for Microwave Measurements[END_REF]. In 1960, a microwave Fabry Perot interferometer at wavelengths around 6 mm is described for frequency or complex permittivity measurements [START_REF] Culshaw | High Resolution Millimeter Wave Fabry-Perot Interferometer[END_REF]. In 1970, broadband frequency operation in the range 8.12-12 GHz is achieved by inserting a variable phase-shifter in the interferometer architecture [START_REF] Gehre | A Microwave Interferometer for Wide-Band Measurements of Complex Transmission Coefficients[END_REF]. In the early 1980s, microwave interferometry is applied to material thickness [START_REF] Ashkenazy | A simple microwave interferometer for measuring effective thickness of dielectric slabs[END_REF] and short distance [START_REF] Thansandote | Microwave interferometer for measurements of small displacements[END_REF] precise measurements. In 2007, microwave interferometry is applied for the first time to the measurement of high impedance devices around 1.8 GHz [START_REF] Randus | A simple method for extreme impedances measurement[END_REF]. In 2011, a mechanical tunable interferometric matching network in the range 1-18 GHz is described [START_REF] Haddadi | Near-field scanning microscopy in liquid media based on microwave interferometry[END_REF]. In 2015, an I/Q-mixer-based interferometric technique is developed for speed operation and noise reduction [START_REF] Vlachogiannakis | An I/Q-mixer-steering interferometric technique for high-sensitivity measurement of extreme impedances[END_REF]. In 2017, an active interferometric method using the internal sources of a multiport VNA and an external coupler is introduced [START_REF] Votsi | An active interferometric method for extreme impedance on-wafer device measurements[END_REF]. These studies have proven that RF interferometric techniques are candidates to tackle the problem of impedance mismatch in a variety of applications. Nevertheless, although efforts have been made to automate the instrumentation and simplify the measurement protocols, conduct such experiments remains a heavy task.

III. INTERFEROMETRIC MODE FOR VNA

Fig. 1 depicts the architecture of the multiport VNA. We connect a reflective dipole at Port 4 with reflection coefficient named source ; this dipole is assumed to be unknown. We propose to choose in this work a 'Short' standard. This dipole only fixes the ratio b4/a4 which is assumed to remain constant when we modify the dipole at the 'Measurement Port'. The wave quantities a1, b1, a4 and b4 correspond to the microwave signals measured at Ports 1 and 4 respectively. The microwave interferometry mode is achieved by modification of external strapping VNA connections as presented in Fig. 1. In this configuration, the Port 1 is the 'Measurement Port' and Port 4 is the 'Reference Source Port'. The internal coupler of Port 3 is used to combine the reflected waves flowing from the DUT (Port 1) and from the reference load (Port 4). And finally, the receiver at Port 2 (noted 'B' in Fig. 1) is used to measure the wave flowing from Port 3. The interferometric procedure consists in the following steps: STEP 1: Connect a one-port reference device at Port 1 with complex reflection coefficient Γref = b1/a1. This reference device will determine the origin of Smith Chart. STEP 2: Turn OFF/ON alternatively sources 1 and 4 to measure the quantities S21 and S24 respectively. We define the complex quantity Ψ (1). We note that the VNA is not calibrated.

Ψ = -

(1) STEP 3: Tune the magnitude and the phase of a4 with respect to a1 by the complex quantity Ψ. The total resulting signal b measured at port B with both RF sources 1 and 4 ON corresponds to 𝑏 = 𝑆 𝑎 + 𝑆 ψ𝑎 = 0

In this configuration, the total emerging wave on port 2 is zero resulting in high sensitivity when measuring a complex reflection close to Γref. It is important to mention that the zeroing condition is directly obtained by only measuring 𝑆 and 𝑆 . If we measure these quantities over a frequency range, we are able to perform the zeroing procedure for the whole frequency points. Compared to microwave interferometers, any tuning of variable attenuator and/or phase-shifter is required to perform this zeroing stage.

IV. EXPERIMENTAL IMPLEMENTATION The photography of the test-set is represented in Fig. 2. The measurement bench involves a four port ®PNA (10MHz-26.5 GHz) from ®Keysight. A 2.92 mm coaxial Short standard is connected at Port 4 via a 2.92 mm flexible coaxial cable. The instrument state of the VNA is following: Fstart = 1.45GHz; Fstop = 1.55GHz; the number of frequency points is set to 401; IF bandwidth = 100Hz; the initial power of the two internal microwave sources is fixed to -5dBm. The device under test (DUT) is built up with mechanical attenuator, phase shifter and short to cover a wide range of impedances. The DUT is connected to Port 1. We chose an impedance of 330  to perform the zeroing at 1.5GHz; the magnitude and phase of the quantity Ψ are 12.31dB and -196.06° respectively. We generate eight resistive impedances between 25 to 1000; these impedances are preliminary measured using the VNA in conventional configuration and the corresponding mechanical positions of the attenuator and the phase shifter are recorded. In Fig. 3 (a), we represent the magnitude in dB versus frequency of the reflection coefficient () of several impedances. Fig. 3 (b) depicts the magnitude in dB of B/R4 versus frequency for the same impedance when the VNA is in interferometer configuration. Finally, Fig. 3 (c) represents the magnitude of S and B/R4 versus the resistance values. The values of B/R4 or B/R1 can be also calibrated using one-port procedure.

To summarize, these results prove that it is possible to realize an auto-zeroing based microwave interferometer using only a multiport VNA without external devices. Indeed, data presented in Fig. 3 demonstrate renormalization of the reference impedance to 330 Ω. In the example proposed, we have chosen a real impedance (phase-shifter is used to keep the phase-shift to of the reflection coefficient phase-shift to 0°) but the reference impedance can be chosen to any desired impedance (real, imaginary or complex).

Fig. 1 .

 1 Fig. 1. Synoptic of VNA in interferometer configuration. (from ®Keysight Technology).

Fig. 2 .

 2 Fig. 2. Photography of VNA in interferometer configuration.

Fig. 3 .

 3 Magnitude of calibrated S11 in dB versus frequency; (b) Magnitude of B/R4 in dB versus frequency; (c) Magnitude of calibrated S11 (linear) and magnitude of B/R4 (linear) versus the real part of impedances.V. CONCLUSION A multiport vector network analyzer configured in interferometric mode is proposed to renormalize the standard reference impedance 50 Ω into any desired impedance. The technique does not require any additional microwave component. The main feature of the present work is that it is possible by including a new option in the firmware to propose a new broadband interferometer based one-port calibration procedure including the auto-zeroing technique. Indeed, vector network analyzers use narrowband tuned receivers to perform measurement frequency point by frequency point. Therefore, the interferometric process can be implemented for each frequency point to provide broad frequency coverage. From the hardware point of view, main limitations in the field of microwave interferometry, beyond the scope of this paper, are the phase noise of the sources and receivers noise. Such implementation will be investigated in further study.

ACKNOWLEDGMENT

This work is performed under EMPIR Planarcal project (http://www.planarcal.ptb.de) . This work has received funding from the European Union Horizon H2020 Programme (H2020-NMBP07-2017) under grant agreement n°761036 MMAMA. This work used the facilities within the EQPX ExCELSiOR (www.excelsior-ncc.eu) funded by the National Research Agency (ANR).