Biomechanical assessment of vertebrae with lytic metastases with subject-specific finite element models

The assessment of risk of vertebral fracture in patients with lytic metastases is challenging, due to the complexity in modelling the mechanical properties of this heterogeneous material. Currently clinical assessment of patients at high risk of fracture is based on the Spinal Instability Neoplastic Score (SINS), which however in many cases does not provide clear guidelines. The goal of this study was to develop a computational approach to provide a comparative biomechanical assessment of vertebrae with lytic lesions with respect to the adjacent controls and highlight the critical vertebrae. The computed tomography images of the thoracolumbar spine of eight patients suffering of vertebral lytic metastases with SINS between 7 and 12 (indeterminate unstable) were analysed. For each patient one or two vertebrae with lytic lesions were modelled and the closest vertebrae without lesions were considered as control. Metastatic vertebrae (N=12) and controls (N=18) were converted to subject- specific, heterogeneous, isotropic, nonlinear finite element models for simulating uniaxial compression. Densitometric and mechanical properties were computed for each vertebra. In average, similar mechanical properties were found for vertebrae with lytic lesions and controls (e.g. ultimate force equal to 6.2 ± 2.7kN for vertebrae with lytic lesions and to 6.2 ± 3.0kN for control vertebrae). Only in three patients the vertebrae with lytic lesions were found to be mechanically weaker (−19% to −75% difference for ultimate stress) than the controls. In conclusion, in this study we presented an approach to estimate the mechanical competence of vertebrae with lytic metastases. It remains to be investigated in a clinical study if this method, together with the SINS, can better classify patients with vertebrae with lytic lesions at high risk of fracture.


Introduction
Vertebral metastatic lesions are very common in cancer patients and over 70% of metastases are located in the spine (Sutcliffe et al., 2013). Breast, lung, prostate, and renal cancers are the most common malignant conditions that lead to the development of lytic lesions on bone (Vialle et al., 2015), which appear in radiological images as focal regions with very low volumetric bone mineral density (BMD) (Sánchez and Sistal, 2014). Other cancers as multiple myeloma and spinal haemangiomas cause a widespread of smaller lytic lesions on bone (Sánchez and Sistal, 2014). The lytic lesions were found to decrease bone strength and increase the risk of fracture (Hardisty et al., 2012;Ebihara et al., 2004). Scoring systems as the Spinal Instability Neoplastic Score (SINS) have been developed to identify patients who need surgical intervention due to the high fracture risk of vertebrae with metastasis (Fourney et al., 2011). However, the SINS fail to identify true negative cases (i.e. specificity equal to 79.5%) (Fisher et al., 2014), leading to the overtreatment of patients already weakened due to the radio-and/or chemo-therapies they need to face against the primary cancer. Moreover, for SINS between 7 and 12 no clear guidelines are reported, making the decision of the clinicians more difficult and based on their experience.
To the authors' knowledge the effect of lytic lesions on the mechanical competence of human vertebrae has been tested experimentally only for mechanically induced defects (i.e. drilled holes in the healthy tissue) Palanca et al., 2018;Alkalay and Harrigan, 2016;Alkalay, 2015;Windhagen et al., 1997;McGowan et al., 1993;Silva et al., 1993). Alkalay (2015)  correlation between stiffness and strength of vertebrae with drilled mechanical defects. Other studies show weak to fair correlations between the size of the induced lytic lesion and the vertebral strength (0.26 ≤ R 2 ≤ 0.52) Windhagen et al., 1997;McGowan et al., 1993;Silva et al., 1993). Palanca et al. (2018) showed how artificially induced defects affect the strain distribution in the external anterior surface of the bone. Alkalay (2015) demonstrated how the mechanisms of vertebral deformation and fracture are affected by mechanically induced lytic lesions. Nonetheless, more studies are needed to better understand the effect of real lytic lesions on the whole vertebral strength. In fact, the tissue around the lesion may remodel over time due to the change in the local tissue environment, changing dramatically the effect of the lesion on the vertebral mechanics.
Subject-specific finite element (FE) models applied to clinical Quantitative Computed Tomography (QCT) scans have been found to be accurate in estimating the structural properties of human vertebrae measured ex vivo (0.28 ≤ R 2 ≤ 0.82 for stiffness and 0.78 ≤ R 2 ≤ 0.86 for ultimate load) (Wang et al., 2012;Dall'Ara et al., 2012;Buckley et al., 2007;Crawford et al., 2003). This approach has been applied to study the structural response of functional spinal units, composed of three vertebrae and two intervertebral discs Alkalay and Harrigan, 2016;Alkalay, 2015;Groenen et al., 2018). QCTbased FE models of two vertebrae one of which included artificially induced lesions have been developed (Alkalay and Harrigan, 2016) and were found to well predict global deformation of the vertebrae with the lesion (R (Vialle et al., 2015) = 0.91) tested under axial compression. Moreover, they showed that the lesion affected the loading transfer between vertebrae, with significant changes in the strain distribution among vertebrae with and without lesions (Alkalay and Harrigan, 2016). Recently, subject-specific QCT based FE models of human vertebral segments with mechanically induced lesions in the middle vertebra have been shown to be able to predict vertebral stiffness (0.64 ≤ R 2 ≤ 0.69) but unable to accurately predict failure loads (0.22 ≤ R 2 ≤ 0.25) estimated from compression tests (Groenen et al., 2018). On the other hand, subject-specific FE models of individual cadaveric vertebrae with lytic lesions based on high-resolution scans resampled to clinical CT resolution have shown to be accurate in predicting vertebral strength (R 2 =0.73 (Vialle et al., 2015)). Despite the clear local microstructural changes induced by the metastases on the bone tissues, it was observed by Nazarian et al. (2008) that the relationship between structural properties of trabecular bone with lytic lesions, measured with experimental compressive tests, and their BMD, computed with a combination of micro-CT scanning and gravimetrical measurements, is similar to that observed for healthy tissue. Therefore, that study suggests that tissue with lytic lesions could be modelled as low-BMD bone tissues, with similar constitutive behaviour. Recently, Lenherr et al. (2018) have confirmed this assumption by showing that there is no significant differences between the material properties of normal and lytic trabecular tissue extracted from human vertebrae and subjected to micro-indentation experiments (N = 14).
FE models of generalised and subject-specific geometries of human vertebrae have been used to study the effect of lesion size, position, bone quality, and other parameters on vertebral strength (Whyne et al., 2001(Whyne et al., , 2003Tschirhart et al., 2004Tschirhart et al., , 2007Galbusera et al., 2018). However, considering the complexity of the vertebral structure, individualised models are required to classify patients at high risk of fracture. Subject-specific heterogeneous QCT-based FE models have been shown to improve the assessment of vertebral stability of patients with multiple myeloma (Campbell et al., 2017). In that study vertebral compressive strength predicted with FE models better classified multiple myeloma patients who suffered of a fracture compared to densitometric or microstructural parameters (Campbell et al., 2017). Nevertheless, the potential of these models in predicting the mechanical stability of vertebrae with lytic lesions remains to be investigated. Moreover, a comprehensive analysis of the mechanical competence of vertebrae with or without lytic lesions in the same patient has not been reported in the literature.
The aim of this study was to use subject-specific QCT-based FE models to evaluate the stability of vertebrae with lytic lesions in patients with vertebral metastases. While a direct validation of the outputs of the models for this specific application was not the goal of this study, the proposed comparative computational approach enables to highlight for each patient if the vertebrae with lytic lesions should be considered critical with respect to the adjacent vertebrae without lesions. After proper validation of the outcomes in clinical studies, this biomechanical approach can be used to support the decision of clinicians when the SINS does not provide clear guidelines.

Patients and scanning
Eight QCT scans were collected from patients (following referred to as P1-P8) with vertebrae with lytic lesions in one or more vertebral bodies classified as indeterminate unstable by the SINS (i.e. scores between 7 and 12) (Fourney et al., 2011). All research approaches were performed in compliance with the local ethical committee ("A multicenter prospective registry for the management and outcome of metastatic spine tumors", 10848-5/2018/EKU). None of the patients was subjected to any radiotherapy session in the six months previous to the QCT scanning. From the eight QCT datasets of the considered patients (three males and five females, 60 ± 12 years old, 70 ± 16 kg weight, and 168 ± 12 cm height), vertebrae with lytic lesions were identified with the help of an experienced orthopaedic surgeon, who assessed the vertebral stability based on the SINS (SINS criteria in Appendix A and results in Table 1).
The QCT scans were acquired with a Hitachi Presto CT machine using an in-line calibration phantom, and a protocol previously defined in the MySpine project (ICT-2009.5.3 VPH) with voltage of 120 kV and intensity of 225 mA (Rijsbergen et al., 2018). Images were reconstructed with a voxel size of 0.6 × 0.6 × 0.6 mm 3 . For each patient, at least one vertebra with a lytic lesion and the two most adjacent control vertebrae (i.e. without lesions) were reconstructed. One (P2, P3, P5 and P8) or two (P1, P4, P6 and P7) vertebrae with lytic lesions were identified and modelled (Table 1).

Image processing
The Hounsfield Units values of the QCT images were converted into BMD equivalent values by using a densitometric calibration obtained with an in-line phantom (Hitachi Presto, Hitachi Medical Corporation, Tokyo, Japan) with five cylindrical insertion with known mean equivalent BMD values (0, 0.5, 0.1, 0.15, and 0.2 g/cm 3 ) (Fig. 1a). Each vertebra was segmented (Amira v6.0.1, Thermo Fisher Scientific, Oregon, USA, Fig. 1b) and aligned with a global reference system, based on an in silico approach described previously in detail (Danesi et al., 2014). Briefly, the vertebrae were aligned with respect to pre-defined anatomical planes, by the selection of 10 virtual landmarks, placed over the most anterior-posterior and medio-lateral regions of the top and bottom surfaces of the endplates, and on the right (LR) and left (LL) corners of the posterior wall of the top endplate (Fig. 1). As the top and bottom endplates of the vertebrae were not parallel to each other the alignment over the sagittal and frontal planes were performed based on the orientation of the respective bisector planes. The local reference coordinate system used for the alignment was set with origin in the LR landmark, x-direction pointing towards the LL landmark, and the xyplane with the same orientation as the transverse bisector plane found between endplates (Fig. 1d).
The cross section area (CSA) of each vertebra was calculated from the binary QCT images as the mean values of the cross section areas calculated in each slice of the vertebral body, excluding the endplates, after manual segmentation (Amira v6.0.1, Thermo Fisher Scientific,  Oregon, USA). The minimum height (Hm) of each vertebra was computed as the axial distance between the most concave points miraof the vertebral endplates (Fig. 1c). The equivalent BMD was estimated for a sub-region of interest of each QCT image, which included the cortical and trabecular bone tissues, in the middle portion of each vertebral body (50% of Hm), excluding the posterior elements. Bone mineral content (BMC) was then estimated as the equivalent BMD times the volume (V) of each considered sub-region of interest.

Finite element models
Each vertebral geometry was meshed using quadratic tetrahedral elements with a maximum edge size of 1.0 mm (more details about the choice of the element size are reported in Appendix B) (Fig. 1e).
Under the assumption that lytic lesions affect only local BMD (Nazarian et al., 2008), bone and lytic tissues were modelled similarly as heterogeneous, isotropic, and elastic-plastic materials, by using the same constitutive law. Elastic properties of the tissue were estimated using a set of density to elasticity relationships from the literature to convert the BMD equivalent value at each element to apparent density [Eq (1)] (Les et al., 1994;Schileo et al., 2008) and then to the elastic modulus [Eq (2)] (Morgan et al., 2003) (Bonemat software, Bologna, Italy). Bone plasticity was modelled using an isotropic Von Mises yield criterion, based on a density-strength relationship [Eq (3)] (Morgan and Keaveny, 2001), and including an isotropic hardening rule, with a 95% reduction in the post-yield modulus [Eq (4)] (Morgan et al., 2003;Bayraktar et al., 2004;Morgan and Keaveny, 2001).   Models were loaded by applying a displacement along the axial direction to the surface nodes of the most cranial endplate until a compression equal to the 1.9% of Hm, assumed as failure of the vertebra (Wang et al., 2012;Keaveny et al., 2014). Nodes from the caudal endplates were fixed in all directions.
To reduce the influence of boundary effects, local properties were analysed for a sub-region of interest of the vertebral models, which included the middle 50% of Hm and excluded all the posterior elements that were farer than 15% of the vertebral anterior-posterior width from the most posterior point of the caudal endplate. Vertebral models had on average 3 million of degrees of freedom and took approximately 2 h to solve in the finite element software Mechanical APDL (ANSYS ® Academic Research, Release 15.0), using parallel distributed memory over a maximum of 32 cores on Iceberg, the High-Performance Computing cluster of the University of Sheffield (3440 cores, 31.8 TB of RAM).
For each vertebral model, the resultant load was computed as the sum of the axial reaction forces calculated at the nodes of the bottom endplate. The applied displacement along the axial direction was computed for the node of the top endplate, which fell closest to its centroid. Such measurements were taken for each iteration of the nonlinear models. Spring stiffness (K) was estimated as the slope of the linear range of the force-displacement curves and ultimate force (F U ) was estimated as the resultant axial reaction force at 1.9% apparent strain (Wang et al., 2012;Keaveny et al., 2014). Work-to-failure (W) was calculated as the area under the load-displacement curve until 1.9% global deformation. The normalized stiffness (E) and strength (σ U ) were computed by normalizing K and Fu by CSA and Hm (Dall'Ara et al., 2012) as following: The yield stress (σ Y ) was estimated from the stress-strain curves with the 0.2% offset method (Morgan and Keaveny, 2001). The energyto-failure (U) was calculated as the area under the stress-strain curve.
Distribution of principal compressive strain (ε p3 ) and stress (σ p3 ) were analysed for each FE model to evaluate how the lytic lesions affect the loading distribution and deformation within each vertebra.
Distributions of equivalent elastic (ε eq el ) and plastic (ε eq pl ) strains were analysed for each FE model to evaluate where the plastic deformation localised within the vertebra.

Estimation of the effect of lesions
For each patient, the effects of lytic lesions were evaluated as percentage differences between the densitometric or mechanical properties of the vertebrae with the lesions, with respect to those calculated for the control vertebrae (i.e. the two most adjacent vertebrae without lesions). Each vertebra with lytic lesion with strength lower than the mean strength value calculated for the controls identified for that vertebra was considered at risk and a biomechanical warning ("FE warning" in the reports) was associated to them. Differences in densitometric and mechanical properties (structural or normalised) predicted between the groups of vertebrae with and without lytic lesions were tested with unpaired two-tails t-test with significance level of 0.05. Linear regressions were used to analyse the relationships between mechanical and densitometric properties or between stiffness and failure load. Moreover, linear regressions between the mechanical properties and the SINS score were also investigated for the vertebrae with lytic lesions. Slope, intercept and coefficient of determination (R 2 ) were reported. Finally, in order to evaluate if a critical threshold in densitometric parameters could be found, for each vertebra with lytic lesion the BMC or BMD were plotted versus the difference in its mechanical properties and those of the respective control.

Results
The values of the geometrical and densitometric properties for each modelled vertebra are reported in Table 2. In particular, the mean equivalent BMD was for all vertebrae (pooled data) ranged from 0.10 g/ cm 3 to 0.33 g/cm 3 (0.20 ± 0.05 g/cm 3 ) ( Table 2) and the mean BMC was 3.28 ± 1.57 g for vertebrae with lytic lesions and 3.20 ± 1.61 g for controls.
The load-displacement and stress-strain curves predicted for the vertebrae with or without lesions show a wide range of mechanical properties (Fig. 2). Predicted ultimate force ranged between 1.7 kN and 12.3 kN (6.2 ± 2.7 kN for vertebrae with lytic lesions, 6.2 ± 3.0 kN for control vertebrae) and vertebral strength ranged between 1.4 MPa and 7.2 MPa (4.5 ± 1.6 MPa for vertebrae with lytic lesions, 4.6 ± 1.3 MPa for control vertebrae). Similar average values of densitometric (volumetric BMD and BMC) or mechanical (K, F Y , F U , W, E, σ Y , σ U , and U) properties were found for vertebrae with or without lesions (p > 0.56) ( Table 3).
Predicted structural stiffness, yield force and ultimate force  correlated well with the mean BMC measured within the vertebral body for pooled data (R 2 =0.75 for structural stiffness, R 2 = 0.84 for yield force and R 2 = 0.82 for ultimate force) (Fig. 3). Slightly better correlations were found between normalized structural properties and BMD for pooled data (R 2 = 0.82 for normalized stiffness, R 2 = 0.85 for yield stress and R 2 = 0.84 for ultimate stress) (Fig. 3). In most cases similar or better correlations between mechanical and densitometric properties were found for vertebrae with metastases (Fig. 3). No significant correlations were found between the SINS values and any of the structural or normalised mechanical properties (p-value between 0.057 and 0.121) (Fig. 3). No clear threshold in densitometric parameters could be identified to classify between vertebrae considered at high risk of fracture (Fig. 4). As expected, predicted normalized structural properties as apparent stiffness and strength were highly correlated for pooled or split data (R 2 = 0.98 for all three cases, Fig. 5). This result is also supported by the localization of the plastic strains in the internal low BMD region of the vertebral bodies for vertebrae with lytic lesions and controls (example in Fig. 7). Considering that homogenous compression has been applied to the superior endplate of each vertebra, the regions (elements) with lower BMD will be the first subjected to higher plastic strains.
The effects of bone lesions on the mechanical properties of the vertebrae were different for different patients (Table 4). Only in three patients (P1, P3, P6) the vertebrae with lytic lesions were found to be mechanically weaker (−19% to −75% difference for ultimate stress) than the controls. In four cases (P2, P4, P5, P8) the vertebrae with the lesions were found to be stronger (8%-88% difference for ultimate stress) than the adjacent control vertebrae without lesions. In one case (P7) one vertebra with lesions was found weaker than controls (−27% difference for ultimate stress), and the other vertebra with lesions was found to be stronger (28% difference for ultimate stress) than controls.
Therefore, for every patient a report was prepared for providing a biomechanical assessment for each vertebra with lesions. The report includes: a sagittal mid-section view of the QCT images, a cross section highlighting the vertebra(e) with lytic lesions, the distribution of BMD within the sub-region of interest of vertebrae with and without lesions, the stress-strain curves predicted with subject-specific FE models of vertebrae with and without lesions and a table with the percentage differences of the normalized structural properties between the vertebrae with lytic lesions and the controls. An example of a report for a patient with vertebrae with critical lesions (P6) is reported in Fig. 8. All the other reports are included in Appendix C.
The distributions of external and internal principal compressive strain and stress were much more homogeneous in the control vertebrae compared to those with lytic lesions, in both cases in which the vertebrae with lesion was found to be stronger or weaker than the control (Fig. 6). However, the denser bone found around the lesion in some cases (e.g. P8) had a shielding effect on the localization of strains which probably lead to increase stiffness and strength.

Discussion
The aim of this study was to develop a computational approach to Fig. 4. Relationship between densitometric paramters and mechanical properties for vertebrae with lytic lesions that showed a lower or higher ultimate force with respect to the controls. Top: linear regressions between BMC or BMD and ultimate force (Fu) for the vertebrae with lower (crosses) or higher (full circles) predicted ultimate force with respect to the adjecent controls. Bottom: scatter plots between BMC or BMD and the difference in ultimate force (Fu) calculated for each vertebra with lytic lesion with respect to adjacent controls.

Fig. 5. Relationship between predicted stiffness and strength.
Linear regression analysis between predicted normalised stiffness and strength for vertebrae with lytic lesions (red circles), and vertebrae without lytic lesions (black circles). The equation for pooled data is reported in gray. evaluate the biomechanical properties of vertebrae with lytic lesions in cancer patients.
None of the vertebrae evaluated in this study had BMD values below the osteoporotic threshold (BMD lower than 0.08 g/cm 3 ), with only two of the vertebrae with lesions (T5 of P3 and L3 of P6) classified as osteopenic (BMD between 0.08 g/cm 3 and 0.12 g/cm 3 ) (Zysset et al., 2015) (Fig. 9). Values of ultimate force obtained from the pooled data were slightly higher than those reported in the literature for human vertebrae without metastatic lesions (Wang et al., 2012;Dall'Ara et al., 2012;Melton et al., 2010) (Fig. 9). Nevertheless, higher variability among the predicted ultimate forces was found compared to previous studies (i.e. coefficient of variation, CV, in this study equal to 46% for pooled data, 49% for healthy and 43% for metastatic vertebrae, versus a CV up to 38% in other studies (Wang et al., 2012)). The higher variability is probably due to the fact that in this study different spine levels (from T4 to L5) were considered. In fact, in the literature only healthy, osteopenic or osteoporotic vertebrae (Wang et al., 2012;Dall'Ara et al., 2012;Melton et al., 2010) from mainly lumbar levels were analysed (T12-L5; L1-L2; or L1-L3).
In line with the literature on vertebral mechanics, structural properties were well correlated with the mean BMC of the vertebral bodies (R 2 = 0.85 for ultimate force versus R 2 = 0.70 in (Dall'Ara et al., 2012); R 2 = 0.80 for stiffness versus R 2 = 0.62 in (Dall'Ara et al., 2012)). As expected, the equivalent BMD correlated better with normalized structural properties (R 2 = 0.86 for ultimate stress versus R 2 = 0.74 in (Dall'Ara et al., 2012); R 2 = 0.82 for normalised stiffness versus R 2 = 0.71 in (Dall'Ara et al., 2012)). These results were extended also for vertebrae with lytic lesions (R 2 = 0.78 for BMC vs ultimate force, R 2 = 0.69 for BMC vs stiffness, R 2 = 0.83 for BMD vs ultimate stress, R 2 = 0.84 for BMD vs normalised stiffness). The similar correlation between BMD and vertebral strength estimated from the FE models for vertebrae with lytic lesions or controls, suggests that the vertebral mechanical properties are driven by the geometrical and densitometric properties of the bone. Therefore, is seems reasonable to use approaches developed to estimate the vertebral strength in osteoporotic subjects for estimating the mechanical properties of vertebrae with lesions. Nevertheless, this is not necessarily true for modelling bones with blastic lesions or primary tumours (e.g. osteosarcoma), for which a similar assessment is required and further experimental analyses to characterize the properties of the tumoral tissue are needed. It should be noted that no correlations were found between the SINS values and the predicted mechanical properties for the vertebrae with lytic lesions. This result suggests that this scoring system is not related to the compressive strength of the vertebra with lesions. The optimal correlation between predicted vertebral strength and normalized stiffness (R 2 = 0.98) is probably due to the simple material model used to describe the post-yield behaviour of bone. In fact, the distribution of the plastic strains showed a localization in the central, low BMD region of the vertebral body, systematically for the different specimens. Similar strong correlations between stiffness and strength have been reported in the literature for human vertebrae without lesions by analysing experimental measurements (R 2 = 0.90) (Pahr et al., 2011) or predictions from FE models (R 2 = 0.92) (Pahr et al., 2011).
The results of this study suggest that a detailed subject-specific biomechanical assessment with computational models can improve our understanding of the criticality of a lytic lesion on the vertebral strength. Moreover, the different results obtained for different patients suggest that patient-specific analyses are recommended for a better estimation of the effect of the lesion on the mechanical stability of the vertebrae. For example, patient P1 had two vertebrae with lytic lesions with SINS equal to 9 for L1 and equal to 7 for L4. Both vertebrae showed slightly lower BMD compared to the controls. The lytic lesion in L1 was larger than that in L4 and it was located in the anterior left region of the vertebral body causing a disruption of the cortical shell ( Figure C1 (Appendix C)). Both vertebrae have shown a lower normalized stiffness and strength compared to the adjacent controls (i.e. reduction of approximately 20%) and, therefore, both were classified as being at risk of fracture. The vertebrae of patient P4 were associated with high values of SINS (10 for L2 and 8 for L3), suggesting that the vertebrae were close to an unstable condition. Nevertheless, from the computational analyses the vertebrae with lytic lesions were considered to be 8%-25% stronger than the adjacent controls, and therefore considered at similar risk of fracture. Another patient (P7), who suffers from an aggressive spinal hemangioma, showed even higher SINS (8 for Table 4 Predicted normalized structural properties (E, σ U , and U) from the subject-specific FE models of a patient set of vertebrae with and without lytic lesions. Percentage differences for each parameter between the estimated values computed for the vertebrae with the lytic lesions and those computed for the control vertebrae were reported for each patient. In the last two columns SINS and indication of critical conditions from the FE models (FE warning) are reported.  Fig. 6. Principal strains. Distribution of minimum principal strain and stress computed for patient P1 for the vertebra with a lytic lesion (L1) and the two adjaent control (T12 and L2). The properties are reported for the 3D distribution (bottom) and for a sagittal section (top). On the left hand side a sagittal section from the QCT image is reported.
T12 and 11 for L1). The lytic lesions, in particular in L1, were widely spread within the vertebral body, which was composed of trabeculae thicker than those observed in other vertebrae, whereas for T12 the lesion was less developed. Only one of the vertebrae was considered at high risk from the computational assessment, the one with the lower SINS (T12 was 27% weaker than the controls, L1 was 24% stronger than the controls). These results highlight that with computational approaches further information about the biomechanical status of the vertebrae can be estimated, something that is not possible with the current clinical tools. There are some limitations in this study. First, a direct validation of the outputs of the FE models for predicting the mechanical properties of the vertebrae with and without lytic lesions was not performed. In fact, the comparative approach used in this study is based on the assumption that the lytic lesions can be modelled as low density bone tissues instead of more complicated constitutive models (Tschirhart et al., 2004;Whyne et al., 2000). Currently, little is known about the properties of the material within the lesions (Whyne et al., 2000) and the composition of the lytic tumour tissues can vary greatly according to the type of primary cancer (Vialle et al., 2015). Nevertheless, the assumption that the lytic lesions can be approximated to low density bone tissues under axial compression seems to be supported by a number of studies Nazarian et al., 2008;Lenherr et al., 2018). In particular, Nazarian et al. (2008), showed that the phenomenological relationships between BMD and elastic properties for trabecular bone tissue with and without lesions are similar, by combining results obtained from micro-CT scans, mechanical tests under compression and gravimetric analyses. More recently Lenherr et al. (2018) have shown that the trabeculae extracted from the human vertebrae in regions affected by lytic lesions or in regions without lesions have material properties, estimated with micro-indentation experiments (N = 14), not significant different (elastic modulus equal to 14.8 GPa for normal tissues and equal to 15.2 GPa for tissue with lytic lesions). Furthermore, Stadelmann et al. (2018) have shown that FE models based on micro-CT images of the human vertebrae with lytic lesions or with small blastic lesions predicted the experimentally measured compressive stiffness and strength similarly to vertebrae without lesions. The three abovementioned studies support the assumption that the material properties of bone tissue with or without lytic lesions are similar and that the same constitutive law can be used to model healthy tissue or tissue with lytic lesions in finite element models of the human vertebrae. Therefore, in our study the subject-specific heterogeneous FE models could be used to evaluate the relative differences in mechanical properties between the vertebrae with or without lesions from the same patients. Nevertheless, before this computational tool can be directly applied in clinics, a direct validation of the approach against ex vivo experiments should be performed. Another limitation of the modelling approach is that the bone tissue was assumed to be isotropic. While vertebral bone has been shown to be anisotropic, this assumption was considered acceptable considering the comparative nature of this study. Furthermore, in this study loading was applied to single vertebral bodies, not considering structures as intervertebral discs and articular contacts between facetjoints, which contribute for the physiologic loading transfer and distribution to the vertebrae (Groenen et al., 2018;Hussein et al., 2012). Finally, only vertebral compressive strength was estimated. In order to account for the possible physiological loading scenarios torsion, bending and multi-axial loading should be simulated. In particular, testing different loading conditions becomes critical when localised lesions may increase the local stresses and strains in critical regions of the bone. Nonetheless, compression is one of the most important loading conditions of the spine which often relates to vertebral fractures (Wang et al., 2012;Buckley et al., 2007;Crawford et al., 2003;Jackman et al., 2015).
In conclusion, in this study we presented an approach to estimate the mechanical competence of vertebrae with lytic metastases in cancer patients. This method has a great potential in improving current approaches to diagnose the risk of fracture of patients with vertebrae with lytic lesions. Nevertheless, it remains to be investigated in a clinical study if this method, together with the SINS, can better classify patients at high risk of fracture compared to the SINS score alone.   Fig. 8. Typical report for biomechanical assessment of vertebral stability. Report for patient P6. The report includes: a sagittal mid-section view of the QCT images of the patient, higlhigting both vertebrae with lytic lesions (a); Distribution of BMD in the chosen region of interest within the vertebrae with lesions (shown by solid and dahsed red lines) and without lesions (shown by solid black and blue lines) (b); Stress-strain curves computed from the subject-specific FE models of vertbrae with lesions (solid and dashed red lines) and without lesions (solid black and blue lines) (c). Percentage differences observed in predicted normalised structural properties (E, σ U , and U) (d).

Conflict of interest
None declared. References. Fisher CG, Charles G. et al., 2010. A novel classification system for spinal instability in neoplastic disease: an evidence-based approach and expert consensus from the Spine Oncology Study Group. Spine 35 (22): E1221-E1229.

Appendix B. Choice of the mesh in the FE models
In this sub-study a mesh size analysis was performed over the heterogeneous QCT-based FE models. The heterogeneous nature of these models implies that there is a dependency between the mesh size and the assignment of the material properties. Therefore, as it is impossible to perform a standard mesh-refinement test, the goal of this sub-study was to evaluate the effect of different mesh refinements on the predictions of both local and structural properties of the vertebrae estimated with the FE models.

Materials and methods
The dataset from one patient (P3) was used. Three vertebrae were modelled: one with a lytic lesion (T5) and two adjacent controls without lesions (T4 and T6). Each vertebra was segmented and aligned as previously reported.
Quadratic (10 nodes) tetrahedral elements were used to discretize each vertebral volume (ICEM-CFD v15.0, Ansys, Pennsylvania, USA). Through a meshing sub-step, the surface mesh of the endplates was isolated from the discretized vertebral volume. In this step, the surface of the endplates was split from the overall vertebral surface object through the manual definition of a cloud of points contouring the endplates. The 3D mesh attached to the surfaces of the endplates was then extracted. The maximum edge size (esize) of the elements used in the mesh refinement study was chosen equal to 0.607 mm, in plane image resolution of the dataset, and by multiplying it to 1.65, to obtain coarser meshes (i.e. esize = 1.00 mm, 1.65 mm and 2.73 mm).
The material properties and the boundary conditions were assigned as described above. The same post-processing approach was used to estimate local and structural properties for each model with different mesh size. Distributions of third principal strains were analysed for each vertebral model to check uniformity in the strain gradients among the different refined models of each vertebra. From each of the most refined vertebral models, the location of the node with the peak value of EPEL3 was considered (i.e excluding surface nodes) and the values in the coarser mesh in that coordinate (displacements and strains) were interpolated using the element shape. The percentage differences (%diff) between the properties estimated from the models at different refinements with respect to those estimated from the most refined model (i.e. reference model) were computed. Convergence was assumed for percentage differences of predicted local and structural properties inferior than 10% (Chen et al., 2014). Stress-strain curves were also computed for each model.
The change in the distribution of material properties within the different refined models of each vertebra were also analysed within the subregion of interest of each model. The tissue elastic modulus and the BMD were calculated as mean among the elements connected to the reference node found in the most refined model, or the values found for the element which contained the coordinates of the reference node in the courser models.

Results
Changes in material properties of the vertebra with a lytic lesion (T5) and the controls (T4 and T6) caused by changes in the size of the maximum edge length were on average smaller than 2.7% for equivalent density and elasticity. As expected, the increase in the size of the elements resulted in a decreased variability of the material properties within the models. Fig. B1. Distribution of the elastic tissue modulus within the region of analysis (left) and stress-strain curves (right) computed for each mesh refinement for the vertebra with the lytic lesion (T5) and for controls (T4 and T6).
Percentage differences in local displacements and compressive strains found between second most refined models (esize = 1.00 mm) and the other models were smaller than 10% for both control vertebrae T4 and T6 (Table B1). However, for the lytic vertebra a percentage difference of 17% was observed for interpolated displacements whereas a percentage difference of 5% was found for interpolated third principal strains (Table B1). Larger differences were found for the courser (Table B1). Table B1 Specifications of the models in terms of size (element size and NDOF), computational expenses (elapsed time and memory usage), direct and indirect convergence of results (EPEL3 at peak nodal location and Fu and K at the structural level). 20 ± 7e-15 −0.099 −0.015 3.8 648 * Elapsed time computed for models running in Iceberg taking 32processors shared among 4 nodes. * * Material properties reported as avg ± std for the most refined model as all the elements containing the peak node of EPEL3. For the coarser models it is only reported the stiffness for the element where the location of the peak node found in the most refined models falls.
Small differences were observed among the stress-strain curves of the different models, with slightly larger difference for the vertebrae with a lytic lesion, T5 ( Figure B1). Predictions of vertebral strength, vertebral stiffness and nodal displacements (%diff lower than 17%) showed in general to be less sensitive than predictions of strains (%diff ranging from 1% to 73% for EPEL3) to the different mesh refinements (Table B1, Figure B1). The percentage differences between the most refined and the second most refined models for stiffness and strength were less or equal to 1% for all vertebrae (Table B1).

Discussion
The aim of this sub-study was to choose the mesh size for the homogenised subject-specific QCT-based heterogeneous models. A small variation was observed for predictions of structural properties, as apparent stiffness and strength, obtained between the refined and the reference models of each vertebra (i.e. %diff < 10%). On the other hand, the differences found for predictions of local properties between different refined modes and the reference models were higher (maximum %diff of 73%).
From the analysis of local stiffness, found within the elements containing the interpolation point in the coarser models, it was observed that in the vertebrae with the lytic lesion the peak location of third principal strains fall in a softer and low density region compared to T4 and T6. This explains the higher local deformation observed for the peak node in T5, of approximately 4%, compared to T4 and T6, peak EPEL3 of approximately 2%. Moreover, it was observed a higher variability within the local elasticity found for the elements containing the peak node in the most refined model of T5 compared to the controls (CV = 1.3 against ∼0-0.6 for T6 and T4 respectively). This finding highlights the high gradient in the material properties within T5 close to the lesion, which explains the higher variability observed in the local properties of this model compared to the controls.
Considering that the modelling from clinical data will be based on the analysis of structural properties, we consider that the models using 1.0 mm element size provide reliable predictions of structural properties, being 73% faster than the most refined models.