Cortical pore size distribution and viscoelastic tibia properties discriminate fragility fractures independent of bone mineral density

Osteoporosis is a disorder of bone remodeling leading to reduced bone mass, structural deterioration, and increased bone fragility. The established diagnosis is based on the measurement of areal bone mineral density by dual energy x-ray absorptiometry (DXA), which poorly captures individual bone loss and structural decay. Enlarged cortical pores in the tibia have been proposed to indicate structural deterioration and reduced bone strength in the hip. Here, we report for the first time the in-vivo assessment of the cortical pore size distribution together with viscoelastic bone properties by means of a novel ultrasonic cortical backscatter (CortBS) technology to discriminate between fractured and nonfractured postmenopausal women (N = 55). The discrimination performance was benchmarked with DXA and high-resolution peripheral computed tomography (HR-pQCT). The results suggest a superior discrimination performance of CortBS (area under the receiver operating characteristic curve: 0.69  AUC  0.75) compared to DXA (0.53  AUC  0.55) and a similar performance compared to HR-pQCT (0.68  AUC  0.73). A model-based measurement of Ct.Th and Ct.Po in radius and tibia bones has been achieved for the first time with the bidirectional axial transmission technology by means of multimode waveguide dispersion analysis 35,45 . The method considers variations of porosity as a major source of variations of cortical bone elasticity, sound velocity and fracture toughness in postmenopausal women 46-48 . Results of a first validation study in postmenopausal women confirmed a comparable fracture discrimination performance of the BDAT variables as aBMD for both vertebral and peripheral fractures 41 . However, axial transmission measurements do not provide direct image-guidance and are restricted to patients with low BMI.


Introduction
Osteoporosis (OP) is an age-associated disorder of bone remodeling leading to reduced bone mass, structural deterioration, and increased bone fragility 1 . According to a recent systematic review of burden and management of fragility fractures in the largest EU countries 2 , fragility fractures caused by OP are the 4 th leading cause of chronic disease morbidity. Although OP is generally thought of as a ''woman's disease'' 3 , men account for a third of OP-related hip fractures in Europe 4 . An estimated 1.0 million quality-adjusted life years (QALYs) were lost in 2017 due to fragility fractures. The current disability-adjusted life years (DALYs) of individuals age >50 years were higher than those for stroke 2 . Established guidelines for the diagnosis of OP recommend the assessment of fracture risk factors and the T-score, which is derived from the measurement of areal bone mineral density (aBMD) by means of dual energy x-ray absorptiometry (DXA) at major fracture sites, i.e., spine and proximal femur 5 . OP is defined for aBMD values 2.5 SD or more below the mean assessed in young adults (i.e., T-score  -2.5). Low bone mass (osteopenia, -1 ≥ T-score > -2.5) is currently not considered as a disease 5 , although a specific osteoprotective therapy is recommended if additional clinical risk factors are present 6 . However, bone strength is determined by a plethora of factors including size, shape, architecture, and composition 7 . Today, there is increasing evidence that the majority of individuals who have sustained an osteoporosis-related fracture or who are at high risk of fracture are not diagnosed as osteoporotic according to the BMD level 8,9 (Fig. 1). 9 . In this longitudinal observational study N=149.524 postmenopausal women (age range: 50 -104 yrs) were evaluated. The colors indicate the established diagnosis thresholds. 82% of women with fractures had a T-score higher than -2.5, i.e., they are not considered to have osteoporosis according to the BMD level. Adapted from 10  Bone tissue undergoes permanent remodeling. Under normal conditions, osteoclasts create resorption canals in the cortical bone tissue matrix, which are refilled by osteoblasts leaving a Haversian canal with a diameter of approximately 30 µm 11 . Bone loss occurs in both women and men as part of the natural aging process 12 . Unbalanced intracortical remodeling typically starts in the endosteal subcompartment and leaves partially refilled or even non-refilled bone multicellular units (BMUs). Clustering of BMUs enhances their chances to merge, creating "giant" pores with diameters larger than 385 µm 13 , leading to the so-called trabecularized cortex 14 , and ultimately to a thinning of the cortical bone shell. However, more than 60 % of intracortical pores are smaller than 100 µm in diameter 14 . In the femoral neck of elderly people (age between 60 and 90 years), decreases of cortical thickness (Ct.Th) and porosity (Ct.Po) by ~4 % and ~32 %, respectively, per decade with no significant changes in pore density (Ct.Po.Dn) have been reported 15 . Decreased cortical thickness and the prevalence of large BMU's have been shown to reduce bone strength 16 and structural deterioration 17 . In contrast, cortical thinning is partially compensated by the apposition of new tissue matrix at the periosteal interface leading to an increased bone diameter and a further increase of the pore diameter gradient in the radial bone direction. Therefore, the cortical pore diameter distribution Ct.Po.Dm.D (hereinafter simply called pore size distribution) represents a relevant target for identifying patients with high fracture risk. However, given their small size, the in-vivo imaging of cortical pores remains challenging, even with the most advanced high-resolution peripheral quantitative computed tomography (HR-pQCT). The imaging resolution of the first and second-generation HR-pQCT systems with voxel sizes of 82 and 61 μm allows direct visualization and segmentation of large pores (i.e., Ct.Po.Dm >100 μm) only, leaving the major fraction of smaller cortical pores unresolved. Iori et al. 18 have proposed a calibration rule for the estimation of Ct.Po locally from volumetric BMD (vBMD) distribution parameters. This method is more accurate (absolute error: 3.4 %) than established vBMD or thresholdbased approaches, as it approximates the contribution of unresolved pores (i.e., Ct.Po.Dm < 80 µm). However, the HR-pQCT technology is not widely distributed and is used mostly in clinical research so far. With DXA-based aBMD, i.e., the diagnostic gold-standard, cortical bone loss resulting from structural decay are poorly captured 9,19,20 .

Fig. 1 Distribution of fracture rate and absolute number of women with fractures vs. bone mineral density (expressed as Tscore values) from the National Osteoporosis Risk Assessment (NORA) study
Quantitative ultrasound (QUS) methods are non-ionizing alternatives for the diagnosis of osteoporosis and the prediction of fracture risk 21 . Many QUS approaches target trabecular sites, e.g. at the heel, and predicting BMD via empirical associations with the measured speed of sound (SOS), broadband ultrasound attenuation (BUA) and other parameters derived from the measurement 22 . More recent QUS technologies, e.g., bidirectional axial transmission measure cortical sites, e.g., distal radius and tibia and aim at the quantitative assessment of structural cortical bone properties, e.g., Ct.Th and Ct.Po 23 . However, none of the existing diagnostic technologies can assess quantitative information about the cortical pore micro-morphology. Particularly, the transition from a normal, age-and gender specific pore size distribution to a pathologically altered one caused by large BMUs could not be assessed in vivo so far.
We have recently developed theoretical cortical bone backscatter model (CortBS) and an ultrasonic multi-angle 3-D acquisition and data processing scheme to assess viscoelastic and microstructural properties in cortical bone 24 . The method measures the frequency-dependent attenuation and backscatter coefficients (f) and BSC(f) at the tibia and retrieves the cortical pore size distribution Ct.Po.Dm.D by fitting a theoretical backscatter coefficient to the measured BSC(f). In an ex-vivo study on bones from 19 human donors, pore size parameters describing the prevalence of large pores could be assessed. The combination of cortical thickness and CortBS parameters provided similar or better prediction accuracies of proximal femur stiffness and strength than aBMD.
In this cross-sectional study the CortBS method was applied for the first time in humans. Postmenopausal women with and without history of fragility fractures were included. We hypothesized that the viscoelastic and microstructural CortBS biomarkers can be assessed in-vivo and that they are associated with the occurrence of fragility fractures. The in-vivo short-term precision of the CortBS biomarkers was assessed and the fracture discrimination performance was compared with those of DXA and HR-pQCT.

Results
Study population. The patient characteristics are summarized in Table 1. Twentynine subjects had at least one fragility fracture. Age and anthropometric data were not different between Fx and nFx groups ( Table 1). More subjects with fragility fractures received antiresorptive treatment than subjects without fractures. The higher number of subjects treated with aromatase inhibitor was not significant (p = 0.06).

DXA:
A valid vertebral spine T-score could not be assessed in 8 subjects due to severe degenerative changes in the lumbar spine. DXA parameters were associated with subject's height (R² = 0.40), weight (R² = 0.34), and almost independent of BMI (R² = 0.23) and age (R² = 0.14) (see Table A.1 in the Supplementary Materials). None of the aBMD values and T-scores was significantly different between Fx and nFx groups ( Table 2), but the difference of the lowest total T-scores between fractured and nonfractured groups almost reached the significance level (p = 0.06). HR-pQCT: Data from one subject could not be evaluated due to an apparent motion artefact. From the remaining 54 subjects, 81 structure and material properties were extracted. HR-pQCT parameters obtained from the scanner software were associated with subject's weight (R² = 0.46), height (R² = 0.42), age (R² = 0.31) and almost independent of BMI (R² = 0.20). Except for BMI, the associations of cortical parameters derived from the custom analysis with anthropometric data and age were generally lower (see supplementary Table A.1). None of the parameters derived from the scanner software was significantly different between fractured and non-fractured groups (Table 3). In contrast, most parameters describing the local distributions of porosity and pore diameter in the antero-medial region of interest were significantly different between both groups. The most prominent differences were observed for skewness (p = 0.004) and kurtosis (p = 0.004) of the pore size distribution Ct.Po.Dm.D evaluated in the full cross sections. CortBS: An ultrasound compound image of the antero-medial region of the tibia of bone together with the normalized difference spectrum and representative backscatter and attenuation coefficients, and pore size distributions for subjects with and without fragility fractures are shown in Fig. 2. The shortterm precision of the individual parameter estimations was in the range between 1.7 and 13.9 % (Table  4). For 5 patients, the quality factor was below 77 % and therefore, data were not analyzed. CortBS parameters were associated with subject's age (R² = 0.45), height (R² = 0.36), and marginally with weight (R² = 0.25) and BMI (R² = 0.22) (see supplementary Table A.1). Parameter ranges and differences between fractured and non-fractured groups are summarized in Table 3. Slope Ct.f and intercept Ct.f values of the attenuation coefficient were significantly different between fractured and non-fractured groups. The change of the width of the pore size distribution Ct.Po.Dm.DFWHM almost reached the significance level (p = 0.06). * p < 0.05

Fragility Fracture Discrimination:
The results of discrimination performance analyses are summarized in Fig. 4 and Table 5. The DXA based T-score values reached the lowest discrimination performance for all types of fragility fractures (0.51  AUC  0.55). Incorporating subject's weight and height could slightly improve the discrimination performance for vertebral and other fractures. Among all HR-pQCT parameters, those describing the shape distributions of porosity and pore diameter were the most predictive ones. Distinct parameter combinations provided good discrimination models for vertebral, non-vertebral, and all fragility fractures (0.65  AUC  0.73). Age and anthropometric information could not further improve the discrimination models. CortBS parameters provided very good discrimination models for all types of fragility fractures (0.65  AUC  0.72). While for non-vertebral fractures only attenuation parameters (Ct.0 and Ct.f) were selected, vertebral and all fractures were discriminated by a combination of attenuation and pore size distribution parameters. The discrimination performance was further improved by weight and height information (0.69  AUC  0.75).  Ultrasound-based BMD prediction: Figure 5 shows the prediction of aBMD from CortBS parameters using multivariate PLS models. Although significant, the correlations were weak to moderate (0.25  R²  0.34).

Discussion
The diagnosis of osteoporosis based on the assessment of structural deterioration of the porous micromorphology has been prevented by the limitations of currently available diagnostic technologies 7 . The established diagnosis of OP based on aBMD captures the relative bone loss of an individual compared to a mean value of young reference population, but not the individual bone loss caused by impaired bone remodeling. Thereby, people with non-pathologically decreased T-scores, but at risk for fragility fractures are currently undiagnosed until the fracture occurs 8 . Alternative diagnostic modalities provide x-ray radiation free BMD and T-score surrogates 25-28 , but do not overcome the principal lack of sensitivity of BMD to the causal microstructural and tissue deteriorations caused by OP. To date, HR-pQCT is the most sophisticated in-vivo imaging modality for the microstructural analysis of bone. It can resolve pathologically enlarged cortical pores, but the spatial resolution is not sufficient to resolve the normal cortical pore micromorphology 18 . In this work, we have applied for the first time in humans a novel ultrasound technology that that provides non-invasively and without ionizing radiation quantitative information about the heterogeneous pore size distribution together with viscoelastic bone properties in cortical bone at the tibia midshaft. The fracture discrimination performance of the novel CortBS technology was benchmarked against conventional DXA-based diagnosis and state-of-the-art x-ray computed tomography (HR-pQCT).
The anteromedial tibia is a favorable site for the assessment of systemic viscoelastic and structural bone tissue deteriorations leading to fragility fractures. The standard locations for DXA measurements (i.e., L1-L4 lumbar spine, and hip and forearm) were initially selected because morbidity from fractures at these locations is high 7 . However, metal implants, previous fractures, scoliosis, osteosclerosis, and aortic calcifications render aBMD estimations at these sites inaccurate or even impossible 29 . The most standardized HR-pQCT measurement site is the distal tibia 30 . In an ex-vivo study, hip failure load has been reported to be associated with low vBMD and microstructural alterations measured at this site 31 . However, towards the epiphyses cortical bone becomes thinner and is increasingly replaced by a trabecular core. Therefore, parameters assessed in this region are susceptible to both positioning errors and inter-subject anatomical differences. In this study, we have used the tibia midshaft region for the following reasons. First, the midshaft contains predominantly cortical bone. Second, the cortical thickness in the antero-medial measurement midshaft region is relatively invariant with respect to the long-axis position and approximately 2 to 3 times larger 32 than reference values reported for the proximal tibia 33 . Thereby, the tibia midshaft provides a much larger and anatomically more invariant tissue volume for cortical bone microstructural characterization than the distal shaft. Recent ex-vivo studies provided evidence that structural deterioration at this measurement site is associated with reduced proximal femur strength 32 and that the biomarkers assessed by CortBS combined with Ct.Th provide superior predictions of proximal femur stiffness and strength compared to aBMD 24 . Ultrasound can be transmitted most easily to and along bone at the facies medialis of the tibia midshaft, where the periosteum is covered by a thin layer of soft tissue only. Already in 1995, Foldes et al. 34 suggested speed of sound measured by axial transmission at the tibia midshaft as an independent predictor of fracture risk in women with non-osteoporotic bone mineral density 34 . Since then, various novel bone QUS techniques have targeted this site for the measurement of Ct.Th, Ct.Po 35,36 , and speed of sound 37 . This study confirmed that i) pore architecture and viscoelastic properties of the cortical tibia midshaft decay with age 12 and ii) are associated with bone fragility. In line with the well-known microstructural deteriorations induced by OP, both HR-pQCT and CortBS revealed predominantly features describing the asymmetry of the cortical pore size distribution rather parameters describing the mean pore-tissue volume fraction as factors associated with fragility fractures. Moreover, frequency-dependent ultrasound attenuation, which is determined by both structural and viscoelastic tissue properties 24 was found to be significantly altered in subjects with fragility fractures.

Discrimination performance:
The results of this pilot study suggest a superior discrimination performance of the ultrasonic cortical backscatter measurement (0.69  AUC  0.75) compared to DXA (0.53  AUC  0.55) and a similar performance compared to HR-pQCT (0.68  AUC  0.73). The two attenuation parameters Ct.0 and Ct.f were the strongest predictors for all types of fragility fractures. Together with subject's height and weight, cortical bone attenuation provided the best discrimination performance for non-vertebral fractures (AUC = 0.69). Subject's height is a known risk factor for nonvertebral fractures 38 , which has been partly linked thinner and more porous cortices in taller women, as measured at the distal tibia by first-generation HR-pQCT 39 . While Ct.Th at the tibia mid-shaft was not a predictive variable in our study, the higher porosity together with other variations of the porosity and pore size distributions were observed both in HR-pQCT and CortBS measurements.
For vertebral fractures, width and 90-% quantile values of the pore diameter distribution were additional predictor variables, while subjects weight remained the only anthropometric factor (AUC = 0.75). This finding is in agreement a previous report suggesting risk factors, e.g. physical weakness, poor health, and weight loss as risk factors for vertebral, but not for non-vertebral fractures 40 .
Our AUC values were lower for DXA and comparable for QUS parameters than those reported in another study, in which cortical thickness and porosity were estimated from axial transmission ultrasound 41 . While in that study on 201 postmenopausal women Ct.Th in was found to be discriminant for hip fractures only (AUC = 0.72), Ct.Po was discriminant for all fractures (AUC = 0.71), for vertebral (AUC = 0.84) and wrist fractures (AUC = 0.71).
Several bone QUS technologies have been used in the past to measure cortical or cancellous bone sites, and at least some of them have demonstrated the potential to predict fracture risk with an equivalent efficiency compared to X-ray densitometry techniques 21,42 . Although ultrasound wave propagation is governed by the structural and material properties of the propagation medium, none of the currently available clinical devices provide any direct measurement of stiffness, strength, or tissue quality. Instead, they provide bone density, stiffness or quality surrogate markers derived from empirical correlations of acoustic properties (e.g., speed of sound SOS and broadband ultrasound attenuation BUA 21 , travel time delays 26,43 , or the shape of the backscatter spectrum 28  CortBS reflects viscoelastic and microstructural deteriorations of cortical bone, which are causally linked to natural aging process and the development of osteoporosis 12 . The crucial role of the porous microarchitecture, particularly the prevalence of large pores as biomarker for reduced bone strength 32 was also confirmed in the HR-pQCT analysis, which revealed the asymmetry of the porosity distribution, but not the total porosity as a fracture discriminating tissue property. In contrast to that ex-vivo study, which included bone from both male and female donors, Ct.Th was not found to be a fracture discriminating biomarker in our study. Limitations. This pilot study has several limitations. First, the cohort size was small and restricted to post-menopausal women with T-scores below -1. The included subjects had diverse fracture and medication histories as well as various co-morbidities. However, the selected cohort resembles the population that is i) most vulnerable for fragility fractures and ii) mostly undertreated based on the BMD diagnosis. Despite these limitations, a strong discrimination performance was achieved, which needs to be confirmed in larger studies covering a larger age range, both gender, larger BMI ranges, and T-scores above -1. Second, the cross-sectional study design did not allow to assess fracture risk. Future prospective studies should therefore evaluate the potential of CortBS parameters to identify people at risk and to assess the individual fracture risk. Second, no real-time assessment of the CortBS measurement quality was possible in this study, which led to the exclusion of data from 5 subjects during the post-hoc data analysis. For clinical applications, the data quality assessment needs to be incorporated into the measurement, providing real-time feedback to the operator and the possibility to repeat the measurement, until an appropriate data quality is achieved.

Conclusion
CortBS is the first quantitative bone imaging modality that can quantify viscoelastic and microstructural tissue deteriorations in cortical bone, which occur during normal aging and the development of osteoporosis. CortBS discriminates fragility fractures in postmenopausal women better than, and independent of aBMD. It could be used as a portable, low-cost, non-ionizing, and widely applicable screening tool to identify people at risk, particularly in the population with low bone mass. A widespread application of the method is anticipated to enable an early identification of people at increased risk, a timely initiation of preventive therapies, and subsequently to a reduction of the prevalence of fragility fractures in people with metabolic bone diseases.

Methods
Research participants. For this cross-sectional study, 55 female subjects (age ≥ 55 years) who have been referred to the Center for Muscle and Bone Research for a clinically indicated DXA bone density measurement were recruited. Height, weight, age, medical history regarding diseases affecting bone health, chronical diseases, fracture status with differentiation regarding adequate/inadequate trauma, medications negatively affecting bone health, as well as osteoprotective and osteoanabolic medications were assessed. To reflect the distribution of fracture rate with respect to BMD in postmenopausal women (Fig. 1), the patient recruitment was stratified into three groups according to the results of the DXA measurement (lowest T-Score of lumbar spine and proximal femur) and fracture status, i.e., OP: osteoporosis (T-Score  -2.5); OPE-Fx: osteopenia (T-Score between -1 and -2.5) and prevalent fragility fracture); OPE-nFx: osteopenia (T-Score between -1 and -2.5) without prevalent fragility fracture. Exclusion criteria were i) BMI > 30; (ii) presence of metal implants or edema at the lower extremity; (iii) no allowance for X-ray exposure; or iv) the inability to understand the nature of the study and follow the instructions. In addition to the measurements on patients, repeated ultrasound readings were performed on three healthy volunteers. The study was registered in the German Clinical Trial Register (DRKS00022217) and was approved by the local ethics committee of the Charité -University Hospital Berlin (reference number: EA4/068/19) and the German Federal Office for Radiation Protection (reference number: Z5-22464/2019-090-G). All participants provided their informed written consent prior to participation.
DXA bone densitometry. DXA (Lunar Prodigy Advance EnCore Software v13.4 or Lunar iDXA EnCore Software v 16.1, GE Medical Systems, Wisconsin, USA) lumbar spine (L1-4) and proximal femur scans were performed as part of the clinical routine examination according to the standard GE LUNAR operator manual. The leg (left or right side) with the lowest aBMD at the proximal femur was defined as index leg for subsequent ultrasound and HR-pQCT readings. aBMD values were assessed at the femoral neck (aBMDFemur(Neck)), the total proximal femur area (aBMDFemur(Total)), and at the spine (aBMDSpine). T-score were calculated, whereas the minimum scores at femur (T-scoreFemur), spine T-scoreSpine, and the lowest value of both anatomical regions (T-scoreTotal) were used for further analyses.
If the T-score from one region could not be assessed, the other one was defined as T-scoreTotal.
Cortical Backscatter (CortBS). The measurement principle has been described in detail previously 24 .
Briefly, a medical ultrasound scanner SonixTOUCH equipped with a SonixDAQ single-channel data acquisition system and a 4DL14-5/38 3-D linear array transducer (Ultrasonix, Richmond, Canada) was used. The system was controlled through a custom-developed user interface. Measurements were performed at the central antero-medial tibia region. The tibia length (LTibia) was assessed as the distance between the medial knee joint cleft and the medial malleolus. Both landmarks were palpated manually. Between these two points, the level of 50% LTibia was marked with a skin marker pencil. The ultrasound transducer was coupled to the skin at this position using an ultrasound coupling pad (aquaflex®, Parker Laboratories, Inc., Fairfield, NJ). Conventional B-mode images were used to position the probe such that a cross-sectional image of the periosteal tibia bone interface appeared in the center of the image. The probe was then manually tilted until the bone surface was approximately normal to the sound beam direction and the focus position Fz was adjusted to be approximately 1 mm below the periosteal bone surface (Fig. 2a). For the measurement, a compound B-mode volume scan sequence was used, which i) scanned a focused beam produced by a 16-element aperture across the array from element position 1 to 128 with an increment of 1, ii) repeated the scan three times with different beam steering angles (-10°, 0°, 10°), and iii) repeated this scan with variable sweep motor positions and respective array tilt angles in the direction perpendicular to the B-mode image plane between  7° with an increment of 1° (Fig. 2a). For each transmit beam, pre-beamformed pulse-echo data were acquired from all 128 elements of the probe. The scan duration was less than 3 seconds. The signal analysis consists of 1) reconstruction of beamformed compound images for all sweep motor positions, i.e., spatial compounding of all three beam steering angles, (Fig. 2b); 2) manual selection of a region of interest covering the bone region to be analyzed (Fig. 2b); 3) automatic detection of the periosteal bone surface within the entire image volume (Fig. 2b); 4) calculation of local beam inclinations, inclination-corrected surface reflection spectrum and an inclination-controlled depthdependent normalized difference spectrum NDS from the pre-beamformed channel data (Fig. 2c); 5), estimation of the frequency-dependent cortical bone attenuation and backscatter coefficients (f) and BSC(f), respectively (Figs. 2d-e); and 6) the estimation of the cortical pore diameter distribution Ct.Po.Dm.D (Fig. 2f). The latter is obtained by minimizing the error between the measured and theoretical BSCs, which are modelled from arbitrary pore size distributions. Except for step 2), all analysis steps were processed fully automatically. A quality parameter, which provides a relative measure of usable data within the selected ROI (i.e., data were not discarded by inclination, signallevel, and other thresholds) was used as an objective criterion to either accept or reject a measurement. Based on repeated measurement with variable probe tilt, a quality score threshold of 77 % was found to produce reproducible results (data not shown). From Ct.(f), slope and intercept values Ct.f and Ct.o were obtained by linear regression (Fig. 2d). Characteristic parameters describing the asymmetric pore diameter distribution (i.e., 10% and 90% quantiles; width, minimum and maximum crossing points of full-width half-maximum (FWHM) values; and peak position were derived. CortBS short-term precision. The short-term precision was evaluated according to 49 by performing 10 repeated measurements with repositioning between each measurement on three healthy volunteers. Absolute and relative precision values were calculated using Eq. (4a) and (5) in 49 , respectively.
High-resolution peripheral computed tomography (HR-pQCT). Immediately following the CortBS measurement, a site matched HR-pQCT scan was performed (XtremeCT II, Scanco Medical AG, Bassersdorf, Switzerland). Subjects were seated in a comfortable, height-adjustable chair. The lower leg of the subjects was positioned carefully in a carbon-fiber cast and fixated in the gantry. Subjects were instructed to sit as still as possible and do not talk or move to avoid motion artefacts. The gantry was moved into the scanner until the skin mark and the laser position indicator were aligned. A total scan length of 10.2 mm in the axial direction divided into 168 cross-sectional images was measured with an isotropic voxel size of 60.7 µm with a scan time of 2 minutes. The total effective dose was less than 5 mSv per scan. A representative reconstructed cross-sectional image is shown in Fig. 3a. Cortical and trabecular properties of the tibia were evaluated using the "3D Density and Structure Analysis" software of the scanner as described elsewhere 30 . Moreover, cortical properties of i) the entire tibia cross-section (Full) and ii) a manually selected antero-medial region of interest (ROI, see Fig. 3) were evaluated using a custom protocol adapted from 18 . This analysis estimates cortical porosity Ct.PoBH using the algorithm proposed by Burghardt 50 , cortical thickness Ct.Th, pore density Ct.Po.Dn, and distributions of porosity Ct.Po.D, pore diameter Ct.Po.Dm, and bone mineral density Ct.BMD.D 18 . From these distributions, characteristic properties, i.e., mean, standard deviation variance, skewness, kurtosis, as well as 10 and 90% quantile values were derived. Statistics. Continuous variables are presented as means and standard deviations (SD). Nonparametric tests were used. Wilcoxon rank sum tests were performed to determine whether parameters were significantly different between the non-fractured and fractured groups. The correlation between ultrasound and x-ray parameters was assessed using Spearman's rank sum correlation coefficient .
To evaluate the correlation of multiple QUS parameters with HR-pQCT and DXA parameters, partial least squares (PLS) regression with three-fold cross validation was used. The coefficient of correlation R² and root mean square error (RMSE) between the predicted QUS based parameter and those measured by HR-pQCT or DXA were computed. The fragility fracture discrimination performance of CortBS, HR-pQCT, and DXA was assessed by means of multivariate PLS discrimination analyses with Leave-One-Out Cross-Validation (PLS-LOOCV) using the libPLS library 51 . For variable selection, a Subwindow Permutation Analysis (SPA) using 10.000 Monte Carlo samplings was repeated until a stable set of significant model variables was found. To avoid overfitting, the final models were created using 3 PLS components. Different discrimination models were developed to predict vertebral, other, and all fragility fractures from DXA based T-scores, HR-pQCT, or CortBS parameters and for combinations with subject's anthropometric data (weight, height, BMI) and age. The area under the curve (AUC) of the receiver operation characteristics (ROC), accuracy, sensitivity, sensitivity, and Odds Ratio (OR) with 95 % Confidence Intervals (CI) were calculated. Except for the PLS-LOOC and SPA analyses, all statistical tests were performed using the Statistics Toolbox of Matlab R2019b (MathWorks, Natick, MA). Statistical results were considered significant for p-values < 0.05.