J.-G Minonzio

N Bochud

Q Vallet

Y Bala

D Ramiandrisoa

H Follet

D Mitton

P Laugier

Bone cortical thickness and porosity assessment using ultrasound guided waves: An ex vivo validation study

Keywords: Elastic guided waves, Cortical bone, Cortical porosity, Cortical thickness, Quantitative ultrasound

Several studies showed the ability of the cortex of long bones such as the radius and tibia to guide mechanical waves. Such experimental evidence has given rise to the emergence of a category of quantitative ultrasound techniques, referred to as the axial transmission, specifically developed to measure the propagation of ultrasound guided waves in the cortical shell along the axis of long bones. An ultrasound axial transmission technique, with an automated approach to quantify cortical thickness and porosity is described. The guided modes propagating in the cortex are recorded with a 1-MHz custom made linear transducer array. Measurement of the dispersion curves is achieved using a two-dimensional spatio-temporal Fourier transform combined with singular value decomposition. Automatic parameters identification is obtained through the solution of an inverse problem in which the dispersion curves are predicted with a two-dimensional transverse isotropic free plate model.

Thirty-one radii and fifteen tibiae harvested from human cadavers underwent axial transmission measurements. Estimates of cortical thickness and porosity were obtained on 40 samples out of 46. The reproducibility, given by the root mean square error of the standard deviation of estimates, was 0.11 mm for thickness and 1.9% for porosity. To assess accuracy, site-matched micro-computed tomography images of the bone specimens imaged at 9 µm voxel size served as the gold standard. Agreement between micro-computed tomography and axial transmission for quantification of thickness and porosity at the radius and tibia

Introduction

Bone strength is determined by a multiplicity of bone characteristics. Bone mineral density (BMD), assessed in vivo by dual energy X-ray densitometry (DXA), remains the most widely used clinical biomarker to predict fracture risk. However, strength-related bone characteristics also include material and micro-architectural properties of trabecular and cortical bone. While the early clinical assessment of skeletal status has been focused mainly on trabecular bone, recent findings about the critical role of cortical bone in bone strength [START_REF] Holzer | Hip fractures and the contribution of cortical versus trabecular bone to femoral neck strength[END_REF][START_REF] Rmd Zebaze | Intracortical remodelling and porosity in the distal radius and post-mortem femurs of women: a cross-sectional study[END_REF] have suggested that fracture risk assessment should also include accurate evaluation of cortical bone [START_REF] Bala | Role of cortical bone in bone fragility[END_REF][START_REF] Cooper | Cortical bone porosity: What is it, why is it important, and how can we detect it?[END_REF].

Porosity and thickness are important cortical bone features in determining its mechanical properties. Cortical porosity has been shown to increase with age [START_REF] V Bousson | Distribution of intracortical porosity in human midfemoral cortex by age and gender[END_REF][START_REF] Cooper | Age-dependent change in the 3D structure of cortical porosity at the human femoral midshaft[END_REF], as well as with the progression of osteoporosis [START_REF] Kk Nishiyama | Postmenopausal women with osteopenia have higher cortical porosity and thinner cortices at the distal radius and tibia than women with normal aBMD: an in vivo HR-pQCT study[END_REF]. It has been found to be correlated with the incidence of vertebral [START_REF] Ostertag | Vertebral fractures are associated with increased cortical porosity in iliac crest bone biopsy of men with idiopathic osteoporosis[END_REF] and femoral [START_REF] Bell | Regional differences in cortical porosity in the fractured femoral neck[END_REF] fractures. In addition, an increasing porosity can result in the trabecularisation of the endosteal surface of cortical bone and thinning of the cortex, further compromising the strength of bone [START_REF] Rmd Zebaze | Intracortical remodelling and porosity in the distal radius and post-mortem femurs of women: a cross-sectional study[END_REF][START_REF] Shigdel | Bone turnover markers are associated with higher cortical porosity, thinner cortices, and larger size of the proximal femur and non-vertebral fractures[END_REF].

Routine in vivo evaluation of cortical porosity and thickness remains challenging with X-ray based imaging technologies. Their assessment is now attainable with the recent advent of the new generation of high-resolution peripheral quantitative computed tomography devices (HR-pQCT) [START_REF] Aj Burghardt | Reproducibility of direct quantitative measures of cortical bone microarchitecture of the distal radius and tibia by HR-pQCT[END_REF][START_REF] Zebaze | A new method of segmentation of compact-appearing, transitional and trabecular compartments and quantification of cortical porosity from high resolution peripheral quantitative computed tomographic images[END_REF][START_REF] Vilayphiou | Agerelated changes in bone strength from HR-pQCT derived microarchitectural parameters with an emphasis on the role of cortical porosity[END_REF]. However, the limited spatial resolution of images, the cost and size of devices and the use of ionizing radiation is a barrier to the use of the technique in clinical routine [START_REF] Bl Jorgenson | A comparison of methods for in vivo assessment of cortical porosity in the human appendicular skeleton[END_REF]. Furthermore, the noninvasive in vivo assessment of the mechanical properties (e.g., elastic anisotropy) of cortical bone remains unattainable to date.

To overcome these limitations, sophisticated quantitative ultrasound (QUS) approaches, refer to as axial transmission techniques, are currently being explored. Several recent studies have suggested that the cortex of long bones (e.g., radius or tibia) behaves as a guide to ultrasonic waves [START_REF] Phf Nicholson | Guided ultrasonic waves in long bones: modelling, experiment and in vivo application[END_REF][START_REF] Lefebvre | Development of a new ultrasonic technique for bone and biomaterials in vitro characterization[END_REF]. Such a waveguide has the characteristic of being multi-mode in the frequency range used (200 kHz-2 MHz), which means that several guided modes coexist.

Axial transmission measurement consists in recording the guided modes that propagate into the cortex in response to an ultrasonic excitation transmitted to its surface, and then studying their frequency response in the form of dispersion curves (represented as the variation of the wave number k = 2πf /c(f ) as a function of the frequency f , c(f ) being the phase velocity of the mode considered). Waveguide characteristics, such as thickness and porosity, can then be deduced from the dispersion curves by fitting a theoretical waveguide model to the experimental data. This stage, referred to as the resolution of an inverse problem, consists in adjusting the parameters of the model until the theoretical and experimental dispersion curves are matched.

Moilanen et al. [START_REF] Moilanen | Assessment of the cortical bone thickness using ultrasonic guided waves: Modelling and in vitro study[END_REF] reported reliable estimates of cortical thickness of human radius measured ex vivo using a free isotropic tube model (with fixed material properties). Similarly, Ta and colleagues retrieved the thickness of bovine tibia by manually fitting measured phase velocities of longitudinal modes to an isotropic hollow tube model filled with a viscous liquid [START_REF] Ta | Measurement of the dispersion and attenuation of cylindrical ultrasonic guided waves in long bone[END_REF][START_REF] Song | Analysis of superimposed ultrasonic guided waves in long bones by the joint approximate diagonalization of eigen-matrices algorithm[END_REF]. In the same vein, Lefebvre et al. [START_REF] Lefebvre | Development of a new ultrasonic technique for bone and biomaterials in vitro characterization[END_REF] identified the Young's modulus of ox bone specimens using an isotropic free plate model (with fixed thickness). In all these studies, only one characteristic of the waveguide (either the thickness or the modulus) was obtained. In addition, they all considered bone as isotropic, whereas it is expected that anisotropy has a significant influence for recovering reliable estimates of cortical bone properties [START_REF] Bossy | Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models[END_REF]. More recently, a pioneering proposal provided a combined estimation of both the cortical thickness and the elastic coefficients from ex vivo measurements [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: A pioneering study on in vitro cortical bone samples[END_REF]. The proposed inverse method, however, required strong prior knowledge to fit the experimental data to the theoretical guided modes. Moreover, only the cortical thickness could be validated as no reference values for the elastic coefficients were available. To face these limitations, our group recently introduced an automated inversion procedure to provide concurrent estimates of cortical thickness and porosity [START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF]. This procedure, applied in vivo at the distal radius of healthy subjects, yielded thickness estimates in good agreement with ground truth values derived from site-matched HR-pQCT measurements [START_REF] Vallet | In vivo characterization of cortical bone using guided waves measured by axial transmission[END_REF]. However, the spatial resolution of the HR-pQCT device used in this study was not sufficient to yield reference values for the porosity.

Here, we go one step further and we validate our approach against traditional microcomputed tomography (µ-CT). The focus of this article is on the ex vivo validation of cortical thickness and porosity estimated at the one-third distal radius and mid-shaft tibia using our custom made axial transmission device. The following two steps were necessary: (i) measurement of the guided modes, and (ii) application of an appropriate and operatorindependent signal processing to automatically extract the dispersion curves and solve the inverse problem. To the authors' best knowledge, this is the first face-to-face comparison of guided wave-based concurrent estimates of cortical thickness and porosity to reference values.

With this information, the goal was to gain better confidence in QUS-based cortical bone biomarkers before implementing them in future clinical trials.

Material and Methods

Samples

Thirty-one excised human radii without soft tissue and fifteen tibiae were measured.

The bone specimens were removed from fresh cadavers and kept frozen (-20 • C) before use.

Seventeen donors were females and fourteen were males (mean age: 78.9 ± 11.9 year; range: 50 -98 year). The bones were provided by the Département d'Anatomie Rockefeller (Lyon, France), through the French program on voluntary corpse donation to science (declaration number DC-2015-2357).

Axial transmission measurements

For each bone, axial transmission measurements were performed on a standardized region of interest (ROI). The center of the probe, aligned along the main bone axis, was placed on a marked position corresponding to the the posterolateral face of the specimens at the distal one-third of the radius (i.e., 70 mm away from the radial styloid, Fig. 1) and the anteromedial face of the specimens at the mid-shaft for the tibia. The length of the ROI, corresponding to the length of the receiver array, was 20 mm. The probe was placed directly in contact with the specimens and coupling between the probe and the bone was ensured with an ultrasonic gel. Before data acquisition, the specimens were gently warmed to room temperature (21 • C). The axial transmission device (Azalée, Paris, France) included an ultrasonic probe (Vermon, Tours, France), a driving electronics (Althaïs, Tours, France) and a human machine interface (BleuSolid, Paris, France). The custom made probe consisted of one array of 24 receivers surrounded by two arrays of 5 emitters each. The reader is referred to our previous publications for detailed characteristics of the probe [START_REF] Bossy | Bidirectional axial transmission can improve accuracy and precision of ultrasonic velocity measurement in cortical bone: A validation on test materials[END_REF][START_REF] Moreau | Accurate measurement of guided modes in a plate using a bidirectional approach[END_REF]. A wideband pulse with a central frequency of 1-MHz (-6 dB frequency bandwidth from 0.4 to 1.6 MHz) was used to excite each emitter. Sixteen radio-frequency signals were averaged by hardware and the resulting signal was sampled at the rate of 20 MHz and then converted by a 12-bit analog-to-digital converter. The human machine interface has been developed to display the spectrum of guided waves in quasi real time (at a frame rate of 4 Hz) and to guide the operator in finding during measurement the optimal position of the probe.

Inference of cortical bone properties

In what follows, we introduce the waveguide model that is used to represent the cortex and the objective function involved in the inverse problem as a metric to compare measured and predicted dispersion curves.

A. Forward problem

The cortex of long bones is modeled as a two-dimensional (2-D) transverse isotropic homogeneous free plate waveguide. Earlier observations guided the choice of this model: (i) in the frequency bandwidth of interest (around 1-MHz), the wavelengths of guided waves, ranging from 1 to 10 mm, are larger than the typical size of the heterogeneous structures (osteons, pores). Cortical bone can thus be considered as a homogeneous propagation medium [START_REF] Baron | Effect of porosity on effective diagonal stiffness coefficients (cii) and elastic anisotropy of cortical bone at 1 MHz: a finite-difference time domain study[END_REF]; (ii) cortical thickness is generally not uniform and may be subject to local variation along the axis of the bone. Given the rather limited length of the ROI (20 mm), we assume a uniform cortical thickness in the measurement region [START_REF] Moreau | Measuring the wavenumber of guided modes in waveguides with linearly varying thickness[END_REF]; and (iii) given our probe configuration (the elements width is much larger than the pitch) and its central frequency (1-MHz), bone curvature is negligible and guided waves mostly propagate in the axial direction. The tubular bone shape can thus be locally approximated by a plate [START_REF] Minonzio | A free plate model can predict guided modes propagating in tubular bone-mimicking phantoms[END_REF].

For such a 2-D waveguide model, guided waves propagation can be described with four independent elastic coefficients c ij (i.e., c 11 , c 13 , c 33 and c 55 ), the mass density ρ, and the thickness Ct.Th, where (1 -2) is the isotropy plane (Fig. 1). Alternatively, the dispersion equations can be formulated as function of bulk wave velocities and elastic anisotropy ratios, without the need to know the mass density which is implicitly included in the velocities [START_REF] Rhee | The group velocity variation of lamb wave in fiber reinforced composite plate[END_REF][START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: A pioneering study on in vitro cortical bone samples[END_REF]:

F A,S f, k; c 13 c 11 , c 33 c 11 , V ⊥ L , V T , Ct.Th = 0 (1)
where V ⊥ L = c 11 /ρ and V T = c 55 /ρ stand for the compression and shear bulk wave velocities in the transverse plan, respectively, c 13 /c 11 and

c 33 /c 11 = V L /V ⊥ L 2
are the elastic anisotropic ratios and V L is the axial compression bulk wave velocity. Subscripts A and S stand for the antisymmetric (A) or symmetric modes (S), respectively. Note that because cortical bone is considered as a homogenous medium, the c ij are the effective stiffness coefficients.

In addition, it has been shown that the effective stiffness coefficients cij can be reasonably predicted from a single parameter [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], cortical porosity (Ct.Po), which represents the average fraction of pore volumes within the cortical bone volume. Parnell et al. [START_REF] Parnell | The influence of mesoscale porosity on cortical bone anisotropy. investigations via asymptotic homogenization[END_REF][START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF] proposed an asymptotic homogenization model in which cortical bone is represented as a two-phase composite material made of a homogeneous transverse isotropic matrix pervaded by periodically distributed water-filled cylindrical pores. This model allows the computation of the effective stiffness coefficients c ij , given the elastic properties (stiffness tensor C m ) and mass density (ρ m ) of the matrix, the properties of the fluid within the pores (C f , ρ f ), and the volume fraction of pores Ct.Po. This representation leads to transversely isotropic elasticity at the mesoscale.

In this study, the effective stiffness is computed by applying this asymptotic homogenization approach, assuming the bone matrix being spatially homogeneous and uniform among individuals. The elastic tensors C m and C f and the mass densities ρ m and ρ f were assumed to be known with fixed values, set according to [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] L and V T ) of Eq. ( 1) can be computed from these values. Note that the model parametrization in terms of porosity Ct.Po does not modify the dispersion equations Eq. ( 1), but can simply be regarded as a sampling of the complete elasticity domain spanned by the four elastic coefficients c ij [START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF]. Therefore, the dispersion equations can now be written as:

F A,S (f, k; Ct.Th, Ct.Po) = 0 . (2)
The

B. Inverse problem

Estimation of model parameters Ct.Th and Ct.Po is obtained through an approach that is an extension of the signal processing applied to extract the experimental guided mode wave numbers from the maxima of the so-called Norm function. The signal processing steps have been extensively described in our previous works [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF][START_REF] Minonzio | Impact of attenuation on guided mode wavenumber measurement in axial transmission on bone mimicking plates[END_REF]. In short, (i) the radio-frequency signals are Fourier transformed with respect to time and stored in a response matrix; (ii) a singular value decomposition (SVD) is applied to the response matrix at each frequency; (iii) signal-to-noise ratio enhancement is achieved by removing the singular vectors associated with the lowest singular values; and (iv) the projection of a testing vector onto the singular vector basis yields the Norm function [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF], defined as:

Norm(f, k) = ||e test (k)|| 2 U (f ) , (3)
where e test denotes the testing vector being a normalized attenuated plane wave [START_REF] Minonzio | Impact of attenuation on guided mode wavenumber measurement in axial transmission on bone mimicking plates[END_REF] and U (f ) stands for the singular vector basis. Similarly, the objective function Proj(Ct.Th, Ct.Po) is defined as the projection of a tested model in the singular vector basis:

Proj(Ct.Th, Ct.Po) = 1 f max -f min fmax f min 1 M M m=1 e test (k m (f, Ct.Th, Ct.Po)) 2 U (f ) df (4)
where f min and f max correspond to the frequency bandwidth limits and M denotes the number of theoretical guided modes.

The Norm function (Eq. ( 3)), obtained through a SVD-enhanced 2-D spatio-temporal Fourier transform, can be interpreted as a guided wave spectrum [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF]. The testing vector spans all the waves measurable by the device, corresponding to all frequencies f and wave numbers k within the device bandwidth. Thus, each pixel (f , k) of the Norm function reflects in a 0-1 scale the presence rate of the tested plane wave in the measured signals. In the case of the objective function (Eq. ( 4)), instead of spanning all measurable waves, the testing vectors are limited to the guided modes of the model. This approach takes advantage of the sparsity of the (f , k) domain, i.e., for a considered model, only a finite number M of guided mode wave numbers k m (f, Ct.Th, Ct.Po) are present at each frequency [START_REF] Xu | Sparse svd method for high-resolution extraction of the dispersion curves of ultrasonic guided waves[END_REF][START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF][START_REF] Zhao | Reconstruction of lamb wave dispersion curves by sparse representation with continuity constraints[END_REF]. Each pixel of the objective function reflects in a 0-1 scale the presence rate of the tested model in the measured signals. Thus, the optimal model parameters Ct.Th and Ct.Po can be found by maximizing the objective function.

Measurement protocol

The measurement protocol consisted of one acquisition resulting from ten measurement sequences without moving the probe. The inversion procedure was applied on each of the ten measurement sequences of an acquisition. For each measurement sequence, the objective function was calculated with a database of waveguide models obtained by sweeping the thickness from 0.8 to 4.5 mm in steps of 0.1 mm and the porosity from 1 to 25% in steps of 1%. These ranges were defined according to literature values for thickness [START_REF] Vallet | In vivo characterization of cortical bone using guided waves measured by axial transmission[END_REF][START_REF] Karjalainen | Ultrasonic assessment of cortical bone thickness in vitro and in vivo[END_REF][START_REF] Bossy | An in vitro study of the ultrasonic axial transmission technique at the radius: 1-MHz velocity measurements are sensitive to both mineralization and intracortical porosity[END_REF] and porosity [START_REF] Vilayphiou | Agerelated changes in bone strength from HR-pQCT derived microarchitectural parameters with an emphasis on the role of cortical 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][START_REF] Bossy | An in vitro study of the ultrasonic axial transmission technique at the radius: 1-MHz velocity measurements are sensitive to both mineralization and intracortical porosity[END_REF][START_REF] Raum | Bone microstructure and elastic tissue properties are reflected in QUS axial transmission measurements[END_REF][START_REF] Ctm Eneh | Porosity predicted from ultrasound backscatter using multivariate analysis can improve accuracy of cortical bone thickness assessment[END_REF]. The grid steps were chosen considering a tradeoff between precision and computation time. Precision was improved to 0.01 mm and 0.1% using a 2-D parabolic fit of the objective function around its maximum value.

The highest value of the objective function was retained when three conditions were fulfilled: (i) the maximum value is not located on a domain boundary, (ii) the maximum value is larger than an heuristic threshold equal to 0.3 (whose aim is to remove poor optimal solutions that are associated with a low value in the 0-1 scale, similarly to the threshold applied on the Norm function [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF]), and (iii) the maximum value is at least 10% higher than the secondary maxima, in order to avoid ambiguity due to multiple local solutions. When one of these three conditions was not met, the measurement was considered as failed. For each acquisition, the ten values of thickness and porosity solution of the inverse problem were averaged.

Measurement reproducibility was assessed by repeating the acquisition four times with intermediate probe repositioning. Finally, the mean of the individual values derived from all successful acquisitions, denoted hereafter as Ct.Th US and Ct.Po US , was compared to the reference values derived from µ-CT.

Reference measurements

Site-matched reference values of cortical thickness and porosity were obtained by imaging the measured region using µ-CT. To this end, a 15 mm-thick cross-section of the posterolateral hemi-diaphysis was cut from the bone region located under the receiver array of the ultrasonic probe (Fig. 3). The specimen was imaged using a desktop micro-computed tomography system (µ-CT, Bruker Skyscan 1176, Kontich, Belgium). Scanning was performed with the specimen immersed in distilled water in a 25 mm inner diameter plastic tube and held by an adhesive paste. The bone longitudinal axis (corresponding to the Haversian canals principal orientation) was aligned to the rotation axis of the sample holder. Scans were achieved with a field of view of 4000 × 4000 pixels, source voltage of 80 kV, current of 309 µA, rotation step of 0.3 • over a 360 • rotation and a 1 mm-thick aluminum filter for beam hardening reduction. An isotropic voxel size of 9.08 µm was used with an exposure time of 1 second, 2 frames averaging, leading to a total scan time of 110 minutes for each specimen. The images were then reconstructed using a filtered back-projection algorithm (NRecon software, V 1.6.9, Skyscan NV, Kontich, Belgium).

For each bone specimen, a stack of 1200 sections (i.e., 11 mm) was reconstructed and analyzed. From the stack of contiguous reconstructed images, the volumes of interest were selected for morphometric analysis using a semi-automated method (CTAnalyser Software V 1.14.4, Skyscan NV, Kontich, Belgium). To this goal, the external borders were contoured by the operator and then automatically interpolated along the longitudinal axis defining the volume of interest (Fig. 3). Bone tissue and pore volumes were segmented using an adaptive threshold (Median-C, radius: 14 pixels, constant: -20, lower gray threshold: 70, upper gray threshold: 230). Reference porosity Ct.Po µCT was quantified as the ratio of the void volume over the total volume of the bone specimen. Cortical thickness was assessed after global thresholding to encompass bone and pore (i.e., gray levels between 0 and 255 on 8 bit-images) by determining the average of the local thickness at each voxel representing solid. Local thickness for a point in solid was defined as the diameter of a sphere that fulfilled two conditions [START_REF] Ulrich | The ability of three-dimensional structural indices to reflect mechanical aspects of trabecular bone[END_REF]: (i) the sphere encloses the point and (ii) the sphere is entirely bounded within the solid surfaces. Finally, the resulting reference thickness Ct.Th µCT was calculated by averaging the values over the whole volume of interest.

Statistical analysis

Reproducibility was estimated by calculating the root-mean-square average of the standard deviation of repeated measurements for each of the measured specimens [START_REF] Glüer | Accurate assessment of precision errors: how to measure the reproducibility of bone densitometry techniques[END_REF]. Correlation analysis was used to compare ultrasound estimates of cortical thickness and porosity with reference values. Normal distribution was tested using Shapiro-Wilk and Kolmogorov-Smirnov tests [START_REF] Ghasemi | Normality tests for statistical analysis: a guide for nonstatisticians[END_REF]. A Bland and Altman analysis was performed to quantify the bias between axial transmission and µ-CT. The difference between the means was tested with one-sample t-test. The level of statistical significance was determined at a p-value below 0.05. Statistical analyses have been performed using the Statistics and Machine Learning Toolbox provided by Matlab (MathWorks, Natick, MA USA).

Results

The inversion was successful for forty specimens (27 radii and 13 tibiae) out of the fortysix that were measured. Six failure cases (4 radii and 2 tibiae), typically associated with poor measurements and the presence of multiple local maxima in the objective function, will be commented in the discussion. Reference values Ct.Th µCT and Ct.Po µCT at the radius fell in the range 1.0 mm to 3.7 mm and 2.3% to 13.3%, respectively. At the tibia, Ct.Th µCT varied from 1.2 mm to 4.3 mm. Ct.Po µCT varied from 3.1% to 16.3%. The reproducibility, given by the root mean square error of the standard deviation of estimates, was 0.11 mm for thickness and 1.9% for porosity. The locations of the ten peak values are shown as black dots and the corresponding means and standard deviations of cortical thickness and porosity estimates are indicated in the title. The standard deviations are less than 0.02 mm and 0.7%. The differences between the estimated values for the different acquisitions are 0.08 mm for the thickness and 2.8% for the porosity. The dynamics of the objective function can be appreciated thanks to the contours ranging from the maximum value minus 0.04 to the maximum value minus 0.01. The area delimited by the highest (thick) line can be interpreted as the measurement resolution, i.e., the ability of the measurement system to differentiate two close waveguide models. For these cases, resolution for the two parameters is estimated to be about ±0.2 mm and ±2%. The optimal matching between experimental (dots) and predicted wave numbers (continuous and dashed lines) is shown in the right column for the first measurement of each sequence. The dispersive part of the modes are clearly visible for phase velocities higher than 4.5 mm.µs -1 .

As the distribution of variables was normal, Pearson's correlations were reported. Agreement between ultrasound-based estimates and reference values range from good (porosity) to excellent (thickness). A strong correlation was found between Ct.Th µCT and Ct.Th US (R 2 = 0.89, p < 10 -5 , RMSE = 0.3 mm, Fig. 5a). The Bland and Altman plot revealed virtually no systematic difference (0.06 mm, p=0.24) between both methods (Fig. 5c). For the cortical porosity, a significant correlation (R 2 = 0.63, p < 10 -5 , RMSE = 1.8%, Fig. 5b) and virtually no systematic difference (0.1%, p=0.86, Fig. 5d) were found between both methods.

The 95% limits of agreement (average difference ± 1.96 standard deviation of the difference), which tell us how far apart measurements by the ultrasound and reference methods are more likely to be for most individuals, are found to be about ±0.6 mm for Ct.Th and ±3.8% for Ct.Po. An excellent agreement is observed for Ct.Th and a relatively overall good agreement is obtained for Ct.Po, although the method performs less well for porosities below 10%. For example, for samples that all have their Ct.Po µCT around 5%, the estimated values Ct.Po US range from 2.0 to 8.5%.

Discussion

We report on an ultrasound axial transmission technique combined with an automated signal processing and data inversion allowing the quantification of cortical thickness and porosity. The determination of the cortical biomarkers was successful for the vast majority of the investigated bone specimens; however it failed for 6 out of the 46 inspected specimens.

To assess the accuracy of bone characteristics retrieval, estimated values were compared to reference values derived from site-matched µ-CT images. Identification of waveguide characteristics by means of guided waves was proven to be feasible, with moderate to excellent agreement between estimates and reference values.

The method developed in the present article consists of two main steps: (i) measurement of the guided modes using our custom made prototype device, and (ii) a signal processing step for determination of the dispersion curves and identification of the waveguide characteristics by inverse problem solving. It is important to note that the whole process of Ct.Th and Ct.Po retrieval (except the step of data acquisition) is automatized.

Some assumptions underlying our model-based approach -namely, that (i) the cortical bone layer behaves like a free plate waveguide of uniform thickness, (ii) the material properties are homogeneous, and (iii) the bone matrix is universal-have been explicitly incorporated into the waveguide model for ultrasonic wave propagation in bone used in the current study.

It should be noted, however, that there is not yet an established empirical consensus on how these assumptions match the cortex of long bones.

The excellent accuracy of Ct.Th retrieval (R 2 =0.89, RMSE=0.3 mm) using a simple homogeneous free plate model allows being confident in the robustness of the proposed method when the computational model used for inverse problem solving does not match the exact properties of the measurement environment, that is to say that it does not account for the bone shape, bone edges irregularity or presence of heterogeneities. The very good accuracy of Ct.Th retrieval is in contrast with the moderate success of Ct.Po retrieval. Despite the general good agreement for Ct.Po indicated by an R 2 of 0.63 and a RMSE of 1.8%, results tell us that the interval corresponding to the limits of agreement (about 7.5%) is large, which, compared to the dynamics of the variability of cortical porosity (about 15%), is important and indicates that the two methods (ultrasound and µ-CT) are not interchangeable.

The development of a robust and accurate acoustic method for material properties determination (here the focus is placed on porosity identification) using guided waves gives rise to challenging difficulties in complex waveguides such as the cortex of long bone. Our approach leads to good results in homogeneous and regular waveguides that can be reasonably approximated by a free plate [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: A pioneering study on in vitro cortical bone samples[END_REF][START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF]. In such a case, several Lamb modes could be clearly identified in the measured dispersion curves (i.e., A 0 , S 0 , A 1 , S 1 , S 2 ). In particular, dispersive branches of high order modes A 1 and S 2 are visible, i.e., phase velocity higher than 4.5 mm.µs -1 (Fig. 2). For these bones, the two-parameter inverse problem delivered reliable and unambiguous solutions and the measurement sequences were stable and reproducible (Fig. 4 and Figs. 6a and6b). Nonetheless, when some of the complexity and heterogeneity of the waveguide cannot be modeled accurately enough by the homogeneous plate model, the identification of waveguide characteristics can result in an ill-posed inverse problem with multiple local solutions or even no solution at all, which could, in turn, explain the moderate success in retrieving accurate Ct.Po estimates and the six failure cases. The analysis of dispersion curves recordings shows that, in some cases, the dispersion curves are incomplete or even entirely lacking for certain modes. In particular, it has been observed that if the low-order antisymmetric mode, i.e., A 0 , is misidentified or lacking, the precision and accuracy of Ct.Po retrieval decreases. In Fig. 6c a larger spread of the objective function can be observed in the porosity dimension, suggesting that a poor measure of the A 0 mode may result in a porosity assessment with a larger uncertainty (e.g., equal to 3.5%). The slope of this mode is known to be mostly sensitive to a variation of porosity for a frequency-thickness product larger than 2 MHz.mm (Fig. 2) [START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF]. In the light of this observation, a particular attention should be given to this Lamb mode during the data acquisition. Notwithstanding, it is noteworthy that this hurdle mostly occurs ex vivo due to the intrinsic complexity of bone geometry under the probe, whereas it was shown that the presence of 4 to 10 mm-thick soft tissue in vivo did neither hamper nor modify the measurement of this guided mode [START_REF] Vallet | In vivo characterization of cortical bone using guided waves measured by axial transmission[END_REF][START_REF] Bochud | Predicting bone strength with ultrasonic guided waves[END_REF].

We speculate that strongly deteriorated bones make the problem particularly difficult to solve. In essence, the notion of waveguide presupposes well-defined waveguide boundaries on which the propagating waves are reverberated, thus giving rise to the presence of guided modes. The disruption of the endosteal bone edge that can be observed in case of strongly deteriorated bones likely weakens the reverberation and subsequently the generation of guided modes. Another issue is the presence of large resorption cavities that has been described in strongly eroded bones of elderly cadavers [START_REF] Rmd Zebaze | Intracortical remodelling and porosity in the distal radius and post-mortem femurs of women: a cross-sectional study[END_REF]. Strong wave scattering may arise from the interaction of ultrasound waves with these large resorption cavities and the superposition of these scattered signals with the guided waveforms may corrupt the analysis and make difficult the retrieval of dispersion curves (Fig. 6d). In other specimens, a strong gradient of porosity from the periosteal to the endosteal cortical bone edges can be observed, and this heterogeneous pores distribution may be a source of misidentification of the cortical porosity. Finally, a last source of complication may be the cortical thickness variation along the propagation direction, which causes reduction of the length of high order modes, as discussed in [START_REF] Moreau | Measuring the wavenumber of guided modes in waveguides with linearly varying thickness[END_REF]. The general complexity of cortical bone appears to result in some deviation from the strict homogeneous plate model and the approximation becomes less satisfying.

A solution to address this complexity would be to conduct 3-D numerical simulations of the wave propagation based on real bone structures derived from 3-D high resolution X-ray computed tomography acquisition [START_REF] Bossy | Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models[END_REF]. Furthermore, somewhat more complicated waveguide models could be incorporated in the forward modeling if appropriate [START_REF] Baron | Using the gradient of human cortical bone properties to determine age-related bone changes via ultrasonic guided waves[END_REF][START_REF] Vavva | Velocity dispersion of guided waves propagating in a free gradient elastic plate: Application to cortical bone[END_REF], at the cost of a more challenging multi-parametric inverse problem analysis [START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF][START_REF] Bochud | Predicting bone strength with ultrasonic guided waves[END_REF]. In spite of the current limitations, the cortical bone estimated parameters appear to be reasonable and might in future studies be shown to be of diagnostic value. Although clinical applications using axial transmission remain challenging for assessing multiple bone biomarkers, as the impact of the soft tissue layer on top of bone is not well established yet, we recently showed that, with our probe configuration and within the inspected frequency range, our experimental data mainly exhibit sensitivity to the influence of cortical bone, in which case a free plate model provides an appropriate inverse model and allows identifying unbiased waveguide estimates [START_REF] Bochud | Predicting bone strength with ultrasonic guided waves[END_REF]. Future clinical studies are warranted to assess whether these cortical bone biomarkers could be used to identify patients at risk of fracture or to judge whether a patient has responded well to treatment.

Another potential complicating factor was that the model of elasticity used to obtain the porosity was based on the assumption of a homogeneous transverse isotropic matrix [START_REF] Parnell | Analytical methods to determine the effective mesoscopic and macroscopic elastic properties of cortical bone[END_REF].

Universal matrix elasticity has been assumed with constant values of density ρ m and stiffness C m , meaning that the properties of the matrix are spatially homogeneous and uniform among individuals. Because the model described here does not account for matrix heterogeneity and diversity, errors are introduced that may influence the parameters estimation, especially that of porosity (indeed, it was shown previously that Ct.Th could be retrieved with reasonable accuracy using fixed effective elastic properties [START_REF] Moilanen | Assessment of the cortical bone thickness using ultrasonic guided waves: Modelling and in vitro study[END_REF][START_REF] Vallet | In vivo characterization of cortical bone using guided waves measured by axial transmission[END_REF]). Currently, it is difficult to determine how much parameter estimates are affected and how the results would be modified if the true matrix properties were used as input into the model. First, porosity estimates depend on the choice of the matrix characteristics used in the model. In this article, we used available experimental matrix stiffness values that have been measured at the femur of aged female donors [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]. Although it is generally accepted that matrix properties do not vary to a large extent, small local differences in matrix properties between cortical bone at the femur and at other skeletal sites may still exist. However, despite this choice of matrix elasticity, no systematic difference could be observed between estimated and reference values (recall Figs. 5c,5d). Alternatively, matrix characteristics could have been chosen according to values available in the literature, either calculated from multi-scale models [START_REF] Morin | A multiscale poromicromechanical approach to wave propagation and attenuation in bone[END_REF] or measured at other skeletal sites (e.g., tibia) [START_REF] Bernard | Elasticitydensity and viscoelasticity-density relationships at the tibia mid-diaphysis assessed from resonant ultrasound spectroscopy measurements[END_REF]. Further studies are warranted to address this issue.

Second, because universal matrix elasticity does not account for inter-individual variability, small deviations that exist between individual matrix characteristics and those of the model may thus explain some of the disagreement between estimated values of porosity and their reference counterpart. For instance, we calculated that a ±5% relative variation of the matrix stiffness coefficients is sufficient to explain a ±3% absolute variation of the porosity estimate.

A solution to overcome the difficulty of choosing matrix elasticity would be to conduct the inversion procedure in terms of stiffness (Eq. ( 1)) as proposed in [START_REF] Foiret | Combined estimation of thickness and velocities using ultrasound guided waves: A pioneering study on in vitro cortical bone samples[END_REF][START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF].

As a further limitation, it should be noted that the measurement was not fully operatorindependent. While data processing was automated, the quality of the data depends on the quality of the scanning exerted by the operator. Our human machine interface, which displays the guided wave measurements in quasi real time, helps the operator to optimize alignment between the probe and the main bone axis, but it did not allow controlling automatically the alignment. Consequently, a misalignment between the probe and the bone axis may result in measurement errors that can partly explain the disagreement between estimated and reference values. Indeed, the most critical part of a measuring sequence is to find the optimal alignment, which is expected to provide the best ultrasonic response. It was observed during the experiments conducted on several bones that the probe alignment with the bone axis is a significant issue in the measurement. One source of improvement would be, for instance, to take a larger number of measurements by continuously scanning the probe around the bone axis. Then, the measurements would be selected using a criterion based on the first arriving signal (FAS) velocity, following earlier findings [START_REF] Bossy | Bidirectional axial transmission can improve accuracy and precision of ultrasonic velocity measurement in cortical bone: A validation on test materials[END_REF][START_REF] Talmant | In vivo performance evaluation of bi-directional ultrasonic axial transmission for cortical bone assessment[END_REF], who suggested that the highest stable FAS velocities correspond to a correct probe alignment.

Conclusion

A model-based inverse problem, in which cortical bone was modeled as a 2-D transverse isotropic free plate with homogenized elastic properties, is presented for the identification of cortical bone thickness and porosity. Cortical bone properties are estimated from the measured dispersion curves of guided modes recorded using an axial transmission technique.

The performance of the inverse problem was demonstrated using reference values of cortical thickness and porosity derived from µ-CT, and results obtained on 40 bone specimens show agreement ranging from good (for porosity) to excellent (for thickness). Our approach failed on 6 specimens. Despite an overall good agreement for porosity, the method performs less well for porosities lower than 10%. The heterogeneity and general complexity of cortical bone structure, which is not fully accounted for by our model, is suspected to weaken the model approximation. As a perspective for our future work we plan to investigate the case of irregular and inhomogeneous bone specimens in order to enhance the estimation of porosity and to test the clinical ability of these cortical bone biomarkers to discriminate between diseased and healthy patients.

Acknowledgment

This work received financial support from Azalée, the Fondation pour la Recherche Médicale (FRM DBS201311228444). Authors would like to thank Rémy Gauthier and Leila Benboubaker for their help during the specimens' preparation and measurement.

Figure 1 :

 1 Figure 1: Measurement configuration with the multi-channel axial transmission probe aligned along the main bone axis of an ex vivo radius specimen. The center of the probe is placed at the distal one-third of the radius and at the mid-shaft for the tibia.

: c m 11 = c m 22 = 26 . 8

 22268 GPa, c m 33 = 35.1 GPa, c m 44 = c m 55 = 7.3 GPa, c m 66 = 5.8 GPa, c m 13 = c m 23 = 15.3 GPa, thus c m 12 = c m 11 -2c m 66 = 15.2 GPa, and ρ m = 1.91 g.cm -3 . In this model, one particular value of porosity Ct.Po determines the four effective elastic coefficients (c 11 , c 13 , c 33 and c 55 ) and the effective mass density ρ = ρ m • (1 -Ct.Po) + ρ f • Ct.Po. The two stiffness ratios (i.e., c 13 /c 11 and c 33 /c 11 ) and the two bulk wave velocities (i.e., V ⊥

  solutions k A,S (f, Ct.Th, Ct.Po) of the previous equations provide the theoretical guided mode wave numbers. The influence of thickness Ct.Th and porosity Ct.Po on the Lamb modes is depicted in Fig.2, for two typical thickness values (1 mm and 2 mm), and two values of porosity (5 and 15%). The number of modes increases with thickness, whereas an increase in porosity is associated with an increase in wave number k (or a decrease in phase velocity c(f ) = 2πf /k)[START_REF] Bochud | Genetic algorithmsbased inversion of multimode guided waves for cortical bone characterization[END_REF].

Figure 2 :

 2 Figure 2: Lamb modes for a 2-D transverse isotropic free plate model with homogenized elastic properties. 1 mm-thick (a) and 2 mm-thick (b) plates with porosities of 5% (black) and 15% (grey). Antisymmetric (A) and symmetric (S) Lamb modes are represented as continuous and dashed lines, respectively. The dotted straight line corresponds to a phase velocity equal to 4.5 mm.µs -1 .

Figure 3 :

 3 Figure 3: An image of a postero-lateral cross-section of a radial hemi-diaphysis obtained with µ-CT. The blue area indicates the ROI used for determining reference cortical thickness and porosity.

Figure 4

 4 Figure 4 illustrates a typical inverse problem solution for three successful acquisitions with intermediate repositioning on one radius specimen. The objective functions (Eq. (4)) of the first measurement of the ten measurement sequences are shown in the left column.

Figure 4 :

 4 Figure 4: Left column: illustration of a typical objective function (Eq. (4)) for three acquisitions with intermediate repositioning on the radius specimen represented on Fig. 3. Right column: comparison between experimental dispersion curves (dots) and the optimal predicted modes (continuous and dashed lines). The values indicated in the title correspond to the mean and standard deviation of the ten measurements shown with black dots.

Figure 5 :

 5 Figure 5: Cortical thickness estimates (Ct.Th US ) plotted against reference values (Ct.Th µCT ) (a) and Bland and Altman plot (c). Cortical porosity estimates (Ct.Po US ) plotted against reference values (Ct.Po µCT ) (b) and Bland and Altman plot (d). Radius and tibia specimens are shown with light gray circles and dark gray circles, respectively. Vertical error bars correspond to the standard deviation from the inferred set of optimal model parameters obtained from the number of successful acquisitions of ten measurement sequences. Dashed lines on (a) and (b) are associated with RMSE and on (c) and (d) with the 95% limits of agreement, whose values are indicated into brackets in the title.

Figure 6 :

 6 Figure 6: Illustration of inverse problem solutions for typical radius specimens: µ-CT transverse crosssection (left), objective function Proj(Ct.Th,Ct.Po) associated with the best acquisition (center) and optimal matching between the experimental dispersion curves (dots) and the guided modes (continuous and dashed lines) (right). The positions of the values obtained by ultrasound and µ-CT are indicated in the (Ct.Th, Ct.Po) plane, with black circles and gray triangles (center), respectively. Contours range from maximum value minus 0.04 to maximum value minus 0.01.