A new finite element based parameter to predict bone fracture

Dual Energy X-Ray Absorptiometry (DXA) is currently the most widely adopted non-invasive clinical technique to assess bone mineral density and bone mineral content in human research and represents the primary tool for the diagnosis of osteoporosis. DXA measures areal bone mineral density, BMD, which does not account for the three-dimensional structure of the vertebrae and for the distribution of bone mass. The result is that longitudinal DXA can only predict about 70% of vertebral fractures. This study proposes a complementary tool, based on Finite Element (FE) models, to improve the DXA accuracy. Bone is simulated as elastic and inhomogeneous material, with stiffness distribution derived from DXA greyscale images of density. The numerical procedure simulates a compressive load on each vertebra to evaluate the local minimum principal strain values. From these values, both the local average and the maximum strains are computed over the cross sections and along the height of the analysed bone region, to provide a parameter, named Strain Index of Bone (SIB), which could be considered as a bone fragility index. The procedure is initially validated on 33 cylindrical trabecular bone samples obtained from porcine lumbar vertebrae, experimentally tested under static compressive loading. Comparing the experimental mechanical parameters with the SIB, we could find a higher correlation of the ultimate stress, σULT, with the SIB values (R2adj = 0.63) than that observed with the conventional DXA-based clinical parameters, i.e. Bone Mineral Density, BMD (R2adj = 0.34) and Trabecular Bone Score, TBS (R2adj = -0.03). The paper finally presents a few case studies of numerical simulations carried out on human lumbar vertebrae. If our results are confirmed in prospective studies, SIB could be used—together with BMD and TBS—to improve the fracture risk assessment and support the clinical decision to assume specific drugs for metabolic bone diseases.


Introduction
Osteoporosis is a skeletal disorder that causes a reduction in bone strength and an increase in fracture risk [1]. It represents today a major public health issue, which affects more than 200 Although FE models based on bi-planar DXA could be a good alternative to replace qCTbased models in the prediction of vertebral strength, they still require a certain dose of radiation and an additional procedure.
Micro-finite element models based on micro-computed tomography (micro-CT) images play a fundamental role in studying and understanding bone failure processes and bone remodeling [23][24][25], however this technique is not currently applied for clinical measurements because images are obtained by scanning ex-vivo specimens.
Literature reports some DXA-based finite element models to evaluate a fracture risk index, but mainly focused on hip fracture [26][27][28]. In particular, Naylor et al. [26] propose a DXAbased FE model tool for the assessment of bone strength, able to identify patients at high risk of hip fracture who may benefit from treatment to reduce fracture risk. Mancuso et al. [27] show that higher levels of adult physical activity, grip strength, and body mass result in a favorable bone microstructure structure and lower numerical strains. Moreover, they also show that areal BMD does not capture the wide range of strains experienced during typical physiologic loading. Moreover, they also show that areal BMD does not capture the wide range of strains experienced during typical physiologic loading. Yang et al. [28] propose an approach based on the 2D DXA images of femur; these authors adopt a mesh with triangular elements and assign homogeneous properties to the tissue.
In the present study, we propose a novel numerical approach based on the Finite Element Method (FEM), derived from DXA greyscale images of density distribution measured on trabecular bone tissue from porcine lumbar vertebrae. The numerical procedure simulates a compressive load on each vertebra to evaluate the local strain values. The greyscale-based mesh and local properties are fundamental and aim to feed the model with additional information, when compared to [28], regarding the local stiffness distribution of the bone. We believe that providing a local distribution of mechanical properties enriches the numerical model, provides a better representation of the reality and ensures a better estimation of fracture, especially considering that it is a local phenomenon. From these values, a novel parameter, named Strain Index of Bone (SIB), is developed and correlated with bone strength. The procedure is first validated on cylindrical trabecular bone samples obtained from porcine vertebrae, experimentally tested under static compressive loading. Finally, we present a few numerical case studies carried out on human lumbar vertebrae. The proposed procedure could represent a complementary numerical tool, based on a patient-specific FE model, to improve the DXA accuracy. The proposed FEstrain based parameter, SIB, incorporates information on the geometry, density distribution from DXA measurements, and loadings and is defined as a local quantity, differently from BMD and TBS, which are averaged over the scanned region. This choice is based on previous literature findings [29][30][31] that showed that failures occurred in local bands and how local discontinuities, porosities, and microstructure could be considered as a source of damage initiation.

Materials and methods
The proposed numerical procedure is applied to porcine samples for a first validation and then to the human vertebrae. This study reports anecdotal examples for a new numerical procedure analysis of bone quality applied to lumbar DXA scans of three patients (Approval of the Local Ethical Committee: Comitato Etico Milano Area 2. Protocol N 2.0 BQ. 265_2017, 13th June 2017). All the three patients provided their written informed consent.

Porcine bone samples
Porcine samples were experimentally tested by a combination of characterization techniques and then were analysed by the proposed numerical procedure.

Experimental procedure.
The experimental procedure is described in detail in a previous work [32] and includes the following steps (Fig 1a-1c): 1. Sample preparation: 40 cylindrical trabecular bone specimens were cored from six different porcine vertebral columns. The ends of bone samples were glued in aluminium endcaps. The nominal diameter of the cylindrical samples was 13.8 mm and the nominal height 30 mm. However, for the following analyses and as in [32], the duplicates, i.e. specimens extracted from the same lumbar vertebra, were removed. This leads to 33 analysed specimens; 2. Analysis of undamaged samples: DXA images were collected adding 20mm layer of solid water provided by Gammex on undamaged samples, by using a Hologic Discovery A system (Hologic Inc, Marlborough, Massachusetts, USA) installed at the Bone Metabolic Unit of the Nuclear Medicine of the Fondazione IRCCS Ca' Granda-Ospedale Maggiore Policlinico, Milan, Italy. Hologic densitometers are dual-energy pulsed systems with pulsed kV between 100 and 140kV. Image segmentation is performed automatically by SW Apex 3.3 installed on Hologic Discovery™. The collected images are projections of the scans anteroposterior, with a DXA resolution of 0.5mm. They were collected on the cylindrical trabecular samples between the endcaps. The purpose of these pre-damage DXA measurements, and the related numerical models described in the next sub-section, is to estimate the failure cross-sections and to quantify damage propensity of the specimens; 3. Mechanical testing: monotonic compressive tests were performed in displacement control (strain rate of 0.0002 s -1 ; constant stroke rate of 0.05 mm/s). Three preconditioning compression cycles up to 0.1% axial strain were performed, followed by monotonic loading until certain strain levels. After that, the specimens were unloaded and loaded three times to the same strain level, to obtain mechanical damage. Samples were divided into four groups, each set loaded until reaching a specific engineering strain value: • Group G1%, with specimens loaded until 1% strain; • Group G2%, with specimens loaded until 2% strain; • Group G3.5%, with specimens loaded until 3.5% strain; • Group G5%, with specimens loaded until 5% strain. From now on, we will also refer to these groups as 'strain levels' or 'damage levels'.
4. Analysis of damaged samples: DXA scanning was performed in air on damaged samples (i.e. after performing the experimental tests). The purpose of these post-damage DXA measurements and the related numerical models is to verify the estimations from the pre-damage DXA and models and to quantify the damage level induced by given mechanical loading.
More details about the experimental testing protocol are provided as Supporting Information.
2.1.2 Numerical procedure. The Finite Element numerical tool used for the simulations is a commercial software, ABAQUS (v. 2017, SimuliaTM, Dassault Systèmes 1 ). The numerical procedure includes the following steps: 1. Model geometry. The numerical model is a simplification of the real 3D cylindrical geometry of the sample. We created a 2D model with a rectangular geometry based on the greyscale matrix of the density, obtained from the DXA scanning (Fig 1d). To build a 2D plane strain FE model, we accounted for a fictitious constant thickness t � , evaluated from the equivalence between the cross-section of the cylindrical real specimen, A s , and the crosssection of the 2D-FE model, A FEM : where r is the radius of the cylindrical specimen, experimentally measured on each sample. Assuming that the two cross-sections are equivalent (A s = A FEM ), we obtained the equivalent thickness t � of the 2D-FE model as in Eq (2): The geometrical model is then discretized into square elements corresponding to the pixels of the DXA image. The DXA output image is imported in Matlab 1 (v. R2018a, The Math-Works, Inc.) as a matrix, whose entries are numbers that indicate the brightness of the pixel (i.e. 0-255 code number, where 0 stands for black colour and absence of bone, and 255 stands for white colour and maximum bone density). Each element of the matrix, corresponding to a pixel of the DXA image, represents a square element of the mesh. This allows one working with a limited number of elements, with a constant size equal to the DXA pixel resolution (500x500 μm). The used elements are 2D continuum plane strain element with four integration points (type CPE4 in ABAQUS).
2. Material characteristics and boundary conditions. The material behaviour is linear elastic. The implemented numerical procedure is not aimed at simulating the damage. Rather the focus is on the correlation between parameters estimated before damage onset and the effective quantities experimentally determined at the beginning of the damage (yielding stress and strain) and at failure (ultimate stress and strain). The choice of a linear elastic model is also motivated by the possibility of reducing the computational time, enabling easy implementation of the numerical tool into the computer used for the DXA in clinical routine. To assign a local value of elastic modulus to each element, we considered a linear proportion between the values of greyscale, G i , and the values of local stiffness, E i , by associating the value of global stiffness, i.e. the elastic modulus E 0 experimentally measured from the compressive experimental tests [32], with the mean grey value � G: where: • � G is the mean greyscale value of the scanned DXA matrix, defined in Eq (3); • G i is the greyscale level corresponding to i-th pixel; • N is the total number of pixels of the scanned DXA matrix, which corresponds to the total number of elements of the 2D model; • i is the index scanning rows and columns of the greyscale matrix; • E i is the elastic modulus of the i-th element, defined in Eq (4); • E 0 is the initial global elastic modulus measured from the experimental compressive test. With this simple procedure, we can obtain a matrix containing the scaled local values of stiffness, E i , according to E 0 and � G. The Poisson ratio is set constant and equal to 0.3. We assigned the same stiffness, E i , to elements corresponding to pixels with the same greyscale value. Fig 2 shows an example of a DXA scan (Fig 2a) and the corresponding FE-model (Fig 2b), where the colour bar indicates the different local elastic moduli assigned to the bone tissue. This procedure allows all the information, obtained from the DXA image (i.e. 2D specimen dimension, distribution of the bone mass, and material properties of each pixel or element), to be included in the FE-model. Boundary conditions are applied to reproduce the uniaxial compressive tests (Fig 2b), by constraining all the degrees of freedom at the bottom edge and applying a vertical force, F, including the mesh, the loading and boundary conditions. Each colour corresponds to a specific stiffness value, as shown in the colour bar. c) Greyscale image obtained from a DXA of a spine from a patient fractured not in the lumbar area (corresponding to P3 in Table 2). b) Example of minimum principal strain (ε pmin ) field in vertebra L1. uniformly distributed to the upper edge. To apply a uniformly distributed force, we used a kinematic coupling. We applied the same load value, F = 60N, to all the samples, by considering that all the animals have similar age and weight. The value of the applied load is selected to obtain stresses and strains of the same magnitude order of the values determined for humans in standing position.
3. Numerical analyses and data post-processing. Linear elastic analyses, aimed at simulating the uniaxial compressive tests, were performed on each model corresponding to a bone sample. From the displacements of the nodes, a numerical elastic modulus, E FEM is calculated, for each sample, by Eq (5): where: • F [N] is the applied compressive force, equal to 60 N; is the height of the model; • u y [mm] is the vertical displacement of the upper surface of sample obtained by numerical analyses.
From the strain values we defined SIB as the minimum principal strain, ε pmin , (maximum in modulus) calculated at each integration point. For each j-th element, having four integration points, we calculated the mean, SIB mean,j , and maximum (in modulus), SIB max,j . The mean and the maximum value over the sample are calculated using Eqs (6) and (7), respectively: where k is the variable scanning the integration points of the whole model, ranging from 1 to 4N, and N is the total number of elements in each model, which corresponds to the total number of pixels in the DXA image. 10 4 is a multiplicative factor to obtain more readable values.
The SIB values are calculated, for each sample, considering the scanned bone region and excluding 5 rows of pixels, at the upper and lower regions, i.e. 2.5mm per side, where the strain field could be affected by the applied boundary conditions. An example of the SIB trend along the specimen height is given in Fig 1e, where the red line highlights the specimen cross section with the highest SIB.
Besides, we also numerically calculated the values of the Strain Index of Bone in correspondence of the yielding strain, SIB Y . SIB Y values are obtained applying to each model the specific experimental yield force, F Y , evaluated at the linear-elastic limit of the stress-strain curve obtained from the corresponding sample.
We performed a statistical analysis to estimate the correlations between the numerical results (SIB max and SIB mean ), and the experimental clinical (BMD and TBS) and mechanical (initial elastic modulus E 0 , yield strain ε Y and ultimate strain at failure ε UTS , and yield stress σ Y and ultimate strength σ UTS ) parameters measured on the same batch of porcine samples and presented in [32]. The software used for the statistical data analysis is Matlab 1 (v. R2018a, The MathWorks, Inc.), and specifically, the function fitlm, fitting a linear model and calculating the main statistical parameters. Correlations are estimated based on the ordinary (unadjusted) determination coefficient R 2 and the adjusted determination coefficient R 2 adj , defined in [33,34].

Human vertebrae: Case studies
We propose the application of a numerical procedure, based on experimental DXA scanning, to human vertebrae, similarly to the porcine specimens described in the previous section. The validation of this procedure can be performed based on a wider experimental campaign with statistical evidence, which is beyond the aim of this work. The case studies we present will allow discussing the methodology and the potential applicability in the clinical routine.
2.2.1 Experimental procedure. Three spines, belonging to i) a non-fractured patient, ii) a lumbar-fractured patient, and iii) a non-lumbar fractured patient, i.e. a patient with osteopenia who presented with an osteoporotic fracture at another site than the lumbar spine, were scanned by the same DXA machine, Hologic Discovery A system, adopted for the porcine bone specimens. 1. Model geometry. Human lumbar vertebrae have different complex morphologies and load capacities as a function of their anatomical position. We assumed that the vertebrae are similar to cylinders, thus excluding the vertebral arc and focusing the attention only on the vertebral body, mainly made of trabecular tissue. Even if this hypothesis can limit the precision of the suggested model, we adapted the procedure used for porcine samples to human vertebrae, to have a rapid and easy tool to be implemented in the post-processing of the DXA measurements. The thickness of the vertebral body, t � , is calculated from the width of each cross section of the vertebra. Assuming a cylindrical geometry of the vertebra, as in [35] the average width represents the mean vertebral diameter. Thus, we can calculate t � or the vertebra as from Eq 2: where: • j is the variable scanning the rows of the greyscale matrix of the processed region of the vertebra; • w(j) is the width in correspondence of the j-th cross section; • P n j¼1 wðjÞ n is the average width, corresponding to the average diameter. Then, in Eq (9) we converted the areal bone mineral density (BMD) into a volumetric bone mineral density (vBMD), based on t � , as suggested by Yang et al. [36]: From the values of vBMD, we estimated by Eq (10) the apparent density, ρ app , as in [37]: 2. Material characteristics and boundary conditions. We derived the stiffness E 0v of each vertebra from Eq (11), which is specific for human bones with a standard error of prediction of 23%, according to previous studies [38][39][40]: To assign a specific elastic modulus to each pixel, according to the greyscale level of the DXA image (Fig 2c), we adopted a similar procedure to the one used for the porcine specimens. We implemented a linear proportion between the elastic modulus, calculated in Eq (11), the mean values of the greyscale � G, calculated as in Eq (3), and the considered pixel, similarly as in Eq (4). The Matlab code, implemented for porcine specimens, was modified and adapted to the case of human vertebrae. Particular attention was paid to the automatic selection of the vertebral area from the DXA output images. The DXA scanning machine evidences automatically the vertebra flanks with white borders, one pixel thick. This is done on the image of the whole spine, which is then cropped. Then, borderlines are separately saved in a binary matrix. In particular, the irregularity of the sides is automatically selected masking the background with specific Matlab commands (i.e. bwselect to create the mask and imshowpair to compare the differences between the original DXA image and the mask). The boundary conditions are the same adopted for the cylindrical samples: a compressive force is applied at a point vertically aligned with the centroid of each vertebral area and located in correspondence of the highest ordinate of the vertebra, and all the degrees of freedom of the bottom edge are constrained. For the loading and boundary conditions, we considered only the elements in common between the studied vertebra and the upper or lower ones, i.e. the intervertebral disk is not considered. This loading case wants to simulate the standing position of the patient. Of course, this is only one of the possible loading modes of the spine, since bending and torsion can occur as well. The proposed numerical model represents a first attempt to simulate the mechanical behaviour of the human and, owing to its two-dimensional nature, we limited the simulations to the simple uniaxial compressive loading. The applied load, F, is specific for each patient and each vertebra of the spine. The calculation of F is based on the work by Han et al. [41], which analysed women with weight 50-120 kg and height from 150-200 cm, covering 99% of the world population. According to [41,42], the load acting on the lumbar vertebrae in the standing posture is approximately linearly dependant on the variation of body weight, rather than on the body height itself.
Thus, the value of F is estimated by interpolating the values in [41] as in Eq (12): where: • W [kg] is the patient's total body weight; • W 1 and W 2 [kg] are the reference extreme weights, respectively of 50 and 120; • F 1 [N] is the resultant force at the considered lumbar vertebra for a person of 50 kg and 150 cm, according to [41]; is the resultant force at the considered lumbar vertebra for a person of 120 kg and 150 cm according to [41].
3. Numerical analyses and data post-processing. Linear elastic analyses, aimed at simulating the uniaxial compression were performed on each model. Similarly to the case of the porcine specimens, we excluded from the post-processing of numerical outputs 10 rows of elements (corresponding to 5mm) from both the upper and from the lower regions (Fig 2d), owing to the strain concentration caused by the boundary conditions. From the FE-simulations, we obtained the strain field and the SIB. In particular, SIB mean and SIB max , using Eqs (6) and (7), respectively. The results were processed considering each vertebra, row by row.  the yield strain, ε Y , versus the numerical Strain Index of Bone calculated in correspondence of the yielding strain, (SIB Y ). The yield strain was experimentally measured during the tests as the linear-elastic limit of the stress-strain curve, whereas the SIB Y values are obtained numerically. The significant linear correlations (R 2 = 0.980) between the numerical and the experimental values of elastic modulus and for both the maximum SIB Y value, SIB Ymax , and the mean SIB Y value, SIB Ymean , with R 2 = 0.834 and R 2 = 0.820, respectively, validate the adopted numerical framework, also confirming the validity of the weighting procedure used to assign the local elastic moduli to each element of the model. This means that the distribution of the local stiffness according to the linear proportion between the elastic modulus and the local bone mineral density given in Eq (4) is representative of the real one.

Porcine bone samples
The strain field is post-processed by computing the average stiffness for each row of elements, to identify the weakest section of each specimen. Fig 4a shows an example of the distribution of the minimum principal strain in the same sample model of Fig 2. Fig 4b and 4c show the values of the estimated elastic modulus, E, and the trends of SIB max and SIB mean per row. The plots in Fig 4b and 4c identify a section where the average stiffness per row is minimum and the corresponding strains are maximum. Since the model is linear elastic and the applied load is the same for all the specimens, there is a net correlation between the minimum value of E, corresponding to the maximum compliance, and the maximum values of SIB. Table 1 summarizes all the adjusted linear determination coefficients among the main clinical parameters (i.e. BMD and TBS), the mechanical properties (i.e. E 0 , ε Y , σ Y , ε ULT , σ ULT ), and the numerical quantities (i.e. SIB). We can consider these determination coefficients as meaningful when the p-value is smaller than 0.05, as given in Table 1. To compute these correlations, we considered all the available specimens, by pooling the data from all the samples belonging to different damage groups. Given the results of Table 1, we can evidence that the two SIB parameters have similar determination coefficients. We will focus on SIB max in the following discussion, because the authors suppose that the damage is a local phenomenon, according to [30,31].  By comparing the analyses performed on pre and post-damaged samples, we noticed that the sections characterized by the highest values of SIB max are the same in pre-and post-damaged samples. More specifically, we found that the row where SIB max is maximum (i.e. the cross-section distance from the bottom grip) was equal or very close (± 1mm) between preand post-damaged specimens.
Fig 5c provides a direct comparison between TBS and SIB max . There is a net separation between pre-and post-damaged samples in terms of TBS, which could be represented by a vertical line at TBS = 1. On the other hand, it is possible to sketch a horizontal line, even if less  Strain index of bone as a parameter for predicting bone fracture evident due to the higher scatter of some post-damaged groups, that could state the presence/ absence of damage at about SIB max = 4. Table 2 shows a comparison between the conventional clinical parameters (BMD, T-score, and TBS) and the proposed one, SIB max , for three patients, including fractured and nonfractured.

Discussion
This study proposes a new parameter, the Strain Index of Bone, the aim being to improve the ability of DXA-based clinical analyses, in predicting the strength and fragility of trabecular bones. The adopted numerical procedure is first validated on the experimental testing carried out on ex-vivo samples taken from porcine vertebrae. Then, it is used to compare the SIB max with the clinical and mechanical parameters and to investigate how the SIB max is influenced by mechanically-induced damage. The determined correlations among clinical, experimental and numerical mechanical parameters (Table 1) evidence that SIB mean and SIB max are the best predictors of the compressive strength of trabecular bones, σ ULT , with determination coefficients, R 2 adj = 0.65 for SIB mean and R 2 adj = 0.63 for SIB max , much higher than the correlation of BMD (R 2 adj = 0.34) and that of TBS, which has a non-meaningful correlation (R 2 adj = -0.03). SIB values are correlated with BMD (R 2 adj = 0.26 for SIB mean and 0.24 for SIB max ), but are not correlated with TBS. For the sake of precision, we should mention that ultimate stress, σ ULT , and ultimate strain, ε ULT , were not available for all the considered specimens, because only specimens belonging to the groups G3.5% and G5% (17 specimens in total) reached the ultimate strength and strain, while the elastic modulus and the yielding strain and stress are calculated also from specimens of groups G1% and G2%.
Moreover, the outcome of this study suggests that SIB is a better predictor of the experimental vertebral strength than BMD and TBS, showing a two-fold increase in the correlation coefficient with respect to BMD. This outcome is in agreement with the recent results of [43]; these authors found that 2D DXA-based numerical models of human spines, compressed with a constant displacement, are able to estimate a load well-correlated with the experimental bone strength (R 2 = 0.66), higher than the single BMD estimation (R 2 = 0.56).
This improvement in the estimation of bone strength could be clinically relevant, with potential improvements in the early diagnosis of osteopenia, as the work [44] did by an experimental campaign, even if with a less detailed 2D model. The yielding properties (ε Y and σ Y ) showed lower correlation coefficients with SIB max than the corresponding ultimate properties (ε UTS and σ UTS ). Bone is a rather complex material, and failure of the trabecular structure normally occurs when the ultimate tensile strength is reached. In all the samples analysed, the highest values of SIB max were found in the regions characterized by the lowest stiffness, thus validating the ability of the proposed parameter in predicting, not only the strength but also the most critical zones of trabecular tissue. SIB max has been shown to be largely sensitive to the damage and that it increases with the damage level. TBS, instead, is affected by the damage and always shows a post-damage reduction, but it is not sensitive to the damage level. Indeed, a large difference in applied strain levels does not correspond to a proportional variation of TBS value. TBS has shown to be a good indicator of the bone quality and of the presence of damage, without being able to quantify it though [45][46][47][48][49][50]. Also, it cannot be used to accurately explain the risk of fracture in the bone tissue or to define how healthy bone is, as it does not show a linear relationship with the mechanical strength [32]. The comparison between SIB max and TBS shows how the latter is able to clearly distinguish between damaged and undamaged samples. The ability of SIB max , in assessing the presence of damage is a little less evident. However, SIB max has shown to be more sensitive to the amount of damage (i.e. strain level), also if compared with other parameters such as BV/TV, BS/TV, Tb.Sp., Tb.Th. and DA, analysed in our previous study [32]. This means SIB max has the potential to be a more quantitative parameter. TBS is directly related to the number of trabeculae and their connectivity and inversely related to the trabecular spacing. However, the number of trabeculae and their connectivity is not always an indicator of bone strength. The ability to withstand the load also depends on the local orientation of the trabeculae, which directly affects the local stiffness and strength. The skeletal architecture is paramount, as local stress concentration may arise, being caused by the trabecular thickness and orientation, and leading to premature failure, independently on bone mass. This study introduces the SIB and supports its validity in damage identification and quantification when applied to the case of ex-vivo trabecular bone samples, from similar porcine spines and having low scattering in the mechanical properties, because pigs were of close age, weight and were all healthy. SIB from undamaged samples offers information about the weakest section, where damage is more prone to occur, i.e. failure propensity of bone; this identification is confirmed by the simulations on post-damaged specimens.
However, one should also point out some limitations of the procedure. The main one consists in reducing a 3D problem to a 2D model. This is due, basically, to the nature of 2D DXA images. As a consequence, the intrinsic thickness variability results in a pattern with a denser region at the core of the sample and less dense regions at the lateral edges (Fig 2b). Nevertheless, this simplified model seems able to catch the strain intensification in the weakest sections of samples and vertebrae.
Moreover, no damage law is included in the linear elastic numerical model. Indeed, the analyses are not aimed at the failure simulation itself, rather the focus is on the correlation between parameters estimated before damage manifestation, i.e. from linear elastic loading, and the beginning of the damage, which can be identified by the yielding. This choice is in accordance with [51]; also, it reduces the computational time and allows a quick run on computers installed on DXA systems.
Despite these drawbacks, the proposed methodology evidences potential clinical applications of the Strain Index of Bone, in particular on human lumbar vertebrae, that are the most stressed along the spine. The application of the procedure to humans is more challenging due to the complex shape of vertebrae, the intrinsic variability of the density, also due to the external cortical shells, and the presence of soft tissues, resulting in a wider greyscale for DXA images. Table 2 compares SIB max values with BMD (bone quantitative parameter accounting for density) and its T-score, and with TBS (bone qualitative parameter accounting for bone texture), commonly used to identify bone metabolic diseases such as osteoporosis. There is a consensual trend between T-score, TBS and SIB max . Indeed, as for BMD and T-score, we notice a remarkable difference between SIB max values of the non-fractured patient (P 1 ) and those of the fractured patients (the osteoporotic P 2 and the osteopenic P 3 ). Similarly, but less evident, TBS is the highest for P 1 , decreases for P 2 , and it is the lowest for P 3 .
The lumbar-fractured patient (P 2 ), experiences the highest value of SIB max at the lumbar vertebra L1, which is located in correspondence of the region where the spine changes its curvature and is statistically more prone to failure. The high SIB max value could indicate that refracture is likely to occur, for this patient, in correspondence of vertebra L1. This consideration could not be drawn based only on BMD and T-score neither on TBS.
We can underline that BMD and T-score agree in the classification of patient P 1 as a healthy condition of the spine, and of the patient P 2 as an osteoporotic condition. On the other hand, P 3 patient is an osteopenic patient according to T-score classification, i.e. between -1 and -2.5 [12]. For this patient, also SIB max evidences average and peak values in between P 1 and P 2 , but it seems interesting to note that the vertebra L1 of P 3 experiences a quite high value of SIB max , in accordance with the lowest TBS. This observation could be an indicator of the local strain state in this particular vertebra. More generally, for patient P 3 , SIB max values seem more in line with the values of P 2 than of P 1 . All these comments could be helpful in the future to assume clinical decisions about specific therapy for metabolic bone diseases.
As a further comment, it seems that SIB max accounts for the density information of BMD and T-score, as well as for the texture information of TBS, adding some details, as for L1 of P 2 .
Clearly, we only provided a few case studies to describe the possibilities of fracture prediction through the proposed numerical procedure. Further data need to be collected before fully assessing, with meaningful statistics, the ability of SIB to predict the event and the location of a failure. At present, DXA data are being processed for healthy patients, as a control group, and for osteoporotic patients, in order to support SIB predictions with statistical relevance. Data from periodic check-ups, combined with the temporal clinical histories of the patients, are also being cross-checked to support the capability of SIB to predict bone failures. In the future, this procedure, based on Strain Index of Bone measurements, could also support, together with BMD and TBS, the clinical prediction of fragility fractures and monitor osteoporosis treatments.