1 Introduction

Bone metastasis generally results from primary tumor invasion to the bone, and femur is one of the most commonly affected sites by metastatic diseases [1,2,3]. A metastatic femur is usually related to a clinical condition characterized by a compromised femur strength, resulting in a reduction of its load-bearing capacity and thus in an increased risk of fracture. Nowadays, the evaluation of fracture risk of a metastatic femur is mainly assessed by the orthopaedic oncologists, on the basis of their experience and sometimes by referring to the Mirels Scoring System (MSS) [4, 5]. This latter furnishes a quantification of the fracture risk by introducing a score based on size, site and nature of the metastatic lesion, as well as on the pain level declared by the patient. However, it has been demonstrated that the MSS is characterized by a low specificity (about 35%), often resulting in an overestimate of the effective fracture risk, and thereby leading to a surgical over-treatment [6,7,8]. This aspect may be mainly explained by the fact that the Mirels scoring entirely neglects the mechanical determinants governing the bone fracture [9]. In this framework, the development of an accurate tool for failure risk assessment in metastatic femurs has to be still considered as an open and challenging task.

In the last two decades, computational approaches based on finite element (FE) method, and accounting for subject-specific mechanical and geometrical features obtained via Computed Tomography (CT) diagnostic techniques, have been successfully employed to predict femur fracture load [10,11,12,13,14,15]. Moreover, in classifying patients at risk of osteoporotic fractures [16,17,18,19], FE modelling revealed a significant improvement in estimating the fracture risk assessment over the available clinical standard (e.g., over the bone mineral density in the case of osteoporosis). As such, recent studies [20,21,22,23] confirm the growing interest in CT-based FE modelling strategies, as potentially effective for overcoming the MSS limitations by allowing to obtain useful quantitative indications of fracture risk in metastatic femurs. Nevertheless, though the healthy bone modeling from CT images has been well established, the proper modelling of the mechanical behavior of a metastatic lesion and of its interaction with the surrounding tissue, is still challenging. In this context, two different strategies have been mainly employed. In particular, the metastasis has been generally modelled i) as a void region [24,25,26] or ii) as a bone tissue, whose constitutive features are roughly estimated by referring to standard correlations for healthy tissues by considering to CT gray-scale values [13, 22, 23, 27]. Accordingly, in the first case, the presence of metastasis mechanics is fully neglected [28, 29], as well as the bone-metastasis interaction or bone degradation at the bone-metastasis interface are not considered. Despite, extrapolating density information for the metastasis by using CT-based grayscale standard correlations, in the second case, leads to treating the lesion as a pseudo-healthy region, disregarding, in this case too, specific constitutive features of the metastatic lesion. In both cases, the complex mechanical microenvironment of metastasis is not properly accounted for, preventing also the possibility to effectively describe mechano-biological stimuli affecting remodelling processes in healthy and diseased tissues [30, 31].

Contributions to the mechanical modeling of tumor mechanics generally rely on poroelastic models, based on the Biot theory [32, 33], in which the tumoral tissue can be mainly described by two phases: the solid matrix and the interstitial fluids [28]. In this context, the characterization and the measurement of mechanical properties of human lytic tumors are crucial aspects due to also the strict connection between the mechanics of the diseased tissue and the malignancy level of the tumor [29]. As a result, the suitable integration in FE modeling strategies of a proper mechanical description of tumoral lesions can significantly affect the corresponding effectiveness in supporting appropriate therapeutic treatments [29]. Some attempts in measuring the mechanical properties of tumor domains have been performed as reported, for instance, in [28, 34,35,36].

A first attempt on developing a FE model of a metastatic vertebra [37] demonstrated that the use of a poroelastic approach, in combination with specific material properties for the tumor, had greater ability in predicting the vertebral burst fracture mechanisms compared to a linear elastic model, and thus it has been proven to be more effective in predicting the corresponding fracture risk. Starting from this result and considering the extensive effort in developing specific models for tumors, we aim at investigating innovative FE modelling strategies to better describe the femoral regions affected by the metastasis compared to current approaches, and thus to furnish a contribute towards a more accurate prediction of fracture risk in metastatic femurs. In such a scenario, present paper aims to introduce a strain-based fracture risk index, evaluated through patient-specific FE analyses and useful to quantify fracture risk in femurs possibly affected by multiple metastases. The driven hypothesis assumes that the metastasis locally induces changes in bone material properties such to detect local phenomena that predispose the femur to rupture. The constitutive modeling is built on the basis of the poroelastic material theory, both for bone and metastasis, and it considers bone-metastasis interaction leading to a more realistic estimate of fracture occurrence. An extensive parametric analysis of poroelastic material properties is performed, highlighting their effect on the fracture risk and providing indications in comparison with the MSS.

The paper is organized as follows. Section 2 describes the CT-based finite element strategy and the constitutive description adopted for bone and metastasis. The parametric analysis based on model parameters is also presented, as well as the FE-based fracture risk indices are defined. Section 3 furnishes the main results in terms of fracture risk levels, strain energy density and Von Mises stress patterns. The results are discussed in Sect. 4. Finally, Sect. 5 presents some concluding remarks.

2 Materials and methods

2.1 Patient imaging data

A bilateral Computed Tomography (CT) scan of a female patient (54 years old, type of primary tumor breast cancer) was retrieved from the Departement of Orthopaedics and Trauma Surgery, at the Campus Bio-Medico University of Rome, Italy. The image analysis revealed the presence of multiple osteolytic metastatic tumors distributed on both (left and right) femurs. As shown in Fig. 1a, three metastases have been identified on both femurs. Specifically, three osteolytic lesions localize in the proximal third of the right femur, whereas two metastases appear in the proximal third of the left femur and one in the middle third. The CT scan was performed with the equipment Brilliance 64 (Philips Medical System, The Netherlands) and involved the entire femurs except for distal condyles (120 kVp, 437 mA, pixel size 0.883mm \(\times\) 0.883 mm).

With the help of an expert physician, the Mirels score has been computed for both femurs of the patient, obtaining a score of 11 and 8 for left and right femurs, respectively. Based on this score and mainly on the judgement of the surgeon, the sole left femur was surgically treated through a procedure of internal nail fixation.

2.2 Image processing

The CT images were segmented to reconstruct both femur geometries (ITK-Snap 3.2.0, University of Pennsylvania), using thresholding followed by manual refinement. The image-derived geometries of each femur are shown in Fig. 1b (light grey overimposed to the DICOM image). The CT scans were further manipulated to extract the surface of each metastasis, by using a manual segmentation process based on the judgement of an expert clinician. In particular, the metastases revealed, with good approximation, a nearly ellipsoidal shape (see Fig. 1c blue regions). Accordingly, each metastasis was described by fitting the segmentation-based surface of the m-th metastasis to an elliptic domain, described by the following inequality:

$$\begin{aligned} \Lambda _{m}(x,y,z)&= A_{m}x^{2}+B_{m}y^{2}+C_{m}z^{2}+D_{m}xy+E_{m}xz \,+F_{m}yz+G_{m}x+H_{m}y+I_{m}z \le 0\,. \end{aligned}$$
(1)

Coefficients \(A_{m},\ldots ,I_{m}\) were numerically determined by a least-square fitting procedure (see e.g. [38, 39]), and thereby resulting in the coordinates of the ellipsoid centroid \(x_{0m},y_{0m},z_{0m}\), and semiaxes \(a_{m}\), \(b_{m}\), \(c_{m}\), for each metastatic region. Finally, the CT images were employed to extract the maximum Hounsfield Unit value (HU), by considering the HU histogram in CT voxels corresponding to the segmented femur and taking the value above which there are only 0.1% of the voxels in the entire femur, similar to [23]. This procedure led to identify a maximum value of HU (denoted as \(\hbox {HU}_{{max}}\)) equal to 1642.

Fig. 1
figure 1

a CT images of metastatic femurs analyzed in this study. Red arrows identify the detected metastatic regions. b Solid models of left and right femurs derived by segmenting the CT datasets. c In blue, representative examples of ellipsoidal domains fitting the CT-based metastatic regions

Full size image

2.3 Finite element computational model

Segmented bone geometries were discretized via unstructured 10-node displacement-based tetrahedral finite elements generated within the Comsol environment (Comsol with Matlab, v.5.2 COMSOL, Stockholm, Sweden) by using an automatic algorithm based on the Delaunay method (Fig. 2a). The average element size d was set equal to 3 mm, as the results of a preliminary convergence analysis (Sect. 3.1). The resulting meshes were characterized by approximatively \(2.5\cdot 10^5\) elements and \(1\cdot 10^6\) degrees of freedom for each femur model. Once the meshes have been generated, their quality has been verified within Comsol. In particular, the algorithm used to create the computational discretization assigns a quality value to each mesh element in the range of 0 to 1. A regular well-shaped tetrahedron has a quality index equal to 1, whereas a totally flat/distorted tetrahedron has a quality index near to 0. In the present work, adopted meshes were characterized by an average element quality of 0.8 with minimum values of the quality index always greater than 0.7, and thereby ensuring accurate and stable numerical solutions.

Inhomogeneous isotropic linearly poroelastic constitutive descriptions were assumed for modeling bone material properties. In detail, according to the Biot theory [40], the second-order Cauchy stress tensor is expressed by:

$$\varvec{\sigma }= {\mathbb {C}}:\varvec{\epsilon }-\alpha p \varvec{I} ,$$
(2)

where \({\mathbb {C}}\) is the drained fourth-order bone elasticity tensor, characterized by minor and major symmetries and positive definite, \(\varvec{\epsilon }\) is the second-order infinitesimal strain tensor, \(\alpha\) is the Biot coefficient, \(\varvec{I}\) is the second-order identity tensor, and p is the pore pressure assumed to be positive when it is compressive. The symbol ‘ : ’ stands for double tensor contraction.

The tensor \({\mathbb {C}}\) is modelled to account for the mechanical interaction between metastatic lesions and surrounding healthy bone, by assuming that the presence of a metastasis induces a transition of local mechanical properties from healthy to metastatic tissues, in agreement with available histological evidence [41]. As a result, the drained fourth-order elasticity tensor can be expressed as:

$${\mathbb {C}}(x,y,z)={\mathbb {C}}_{CT}(x,y,z)+\sum _{m=1}^N {\mathcal {H}}_{m}(x,y,z) f_{m}(x,y,z)[{\mathbb {C}}_{l}-{\mathbb {C}}_{CT}(x,y,z)] ,$$
(3)

where N identifies the number of metastatic regions clinically detected in each femur, \({\mathbb {C}}_{CT}(x,y,z)\) represents the fourth-order elasticity tensor derived from CT images, completely described in terms of Young’s modulus (\(\hbox {E}_{{CT}}\)) and Poisson’s ratio (\(\nu\)). According to [42], \(\nu =0.3\) was assumed homogeneously distributed through the femur. Conversely, an heterogeneous Young’s modulus distribution was considered by employing the following stiffness-density correlation [43]:

$$\rho _{app}= 1.1\cdot 10^{-3} \mathrm{HU} \qquad [{\mathrm{g/cm}}^{3}]$$
(4a)
$$E_{CT}= 3790\rho _{app}^{3} \qquad \qquad \,\, [\mathrm{MPa}]$$
(4b)

with \(\rho _{app}\) being the apparent density. Equation (4a) denotes the calibration of the CT datasets, that should be properly obtained by using a calibration phantom. In the present work, due to the lack of a specific calibration phantom, the relation between HU and density has been established by considering that:

  • HU=0 corresponds to a null value of \(\rho _{app}\);

  • the maximum value of HU (equal to 1642) corresponds to \(\rho _{app}=1.8 \,{\mathrm{g/cm}}^{3}\), representing the maximum value of the apparent density for cortical bone [44] (see Sect. 2.2).

For the sake of example, Fig. 2b shows the corresponding distribution of \(E_{CT}\) for the left femur model, at the coronal section.

The tensor \({\mathbb {C}}_{l}\) in Eq.( 3) refers to the elastic properties of metastasis, defined through the Young’s modulus \(\hbox {E}_l\) and Poisson’s coefficient \(\nu _l\). These elastic coefficients were set equal to 0.003 MPa and 0.1, respectively [37] .

Moreover, functions \({\mathcal {H}}_{m}(x,y,z)\) and \(f_{m}(x,y,z)\) in Eq. (3) are the Heaviside function and a Gaussian-like shape function, respectively, referred to the m-th metastatic lesion, that is

$${\mathcal {H}}_{m}(x,y,z)=\left\{ \begin{array}{ll} 1 &{} \qquad \text {if } \quad \Lambda _{m}(x,y,z)\le 0\\ 0 &{} \qquad \text {if } \quad \Lambda _{m}(x,y,z) > 0 \end{array}\right.,$$
(5)
$$f_{m}(x,y,z)= \exp \left[ -\left( \frac{(x-x_{0m})^{2}}{a_m^2}+ \frac{(y-y_{0m})^{2}}{b_m^2}+\frac{(z-z_{0m})^{2}}{c_m^2}\right) \right].$$
(6)

It is worth observing that, at the m-th metastatic boundary \(\Lambda _m=0\), the function \(f_m\) assumes the value 0.37. Such a condition corresponds to prescribe a difference in the elastic properties between healthy and metastatic tissues facing on the interface \(\Lambda _{m}=0\) in the order of 37%, in agreement with the experimental evidence reported in [45] in the case of trabecular bone. The prescribed constitutive relations allow to simulate a graded degradation of mechanical properties induced by the metastasis in the bone tissue adjacent to the lesion (Fig. 2c).

It is worth also mentioning that a purely elastic constitutive description not accounting for the bone-metastasis mechanical interaction can be straight retrieved by considering \(\alpha =0\) in Eq. (2) and \(f_m(x,y,z)=0\) in Eq. (3).

Fig. 2
figure 2

Sketch of the image-based FE modeling strategy. a Tetrahedral unstructured mesh derived from the CT-based solid model. b Heterogeneous Young’s modulus distribution (\(\hbox {E}_{{CT}}\)) derived from the CT dataset. c Young’s modulus distribution (labelled as \(\hbox {E}_{{new}}\)) obtained by accounting for the bone-metastasis interaction (see Eq. 3). The red square shows a metastatic region, highlighting the change of the Young’s modulus distribution due to the presence of the metastasis. d Boundary conditions employed for numerical analyses: a compressive force (\(\hbox {F}_{{z}}=1000 \, \hbox {N}\)) distributed on the femoral head, and the distal portion of the femur fully constrained (3.5 cm long in the z-direction)

Full size image

2.4 Boundary conditions

The CT-based femoral computational models were employed to investigate the femoral mechanical response under a physiological quasi-axial load, describing a “stance” condition. In detail, by referring to the Cartesian framework introduced in Fig. 2d, a compressive force (along the z-axis) of 1000 N was distributed over a pseudo-circular surface region (with a radius of about 10 mm) at the top of the femoral head. The value of the applied force was chosen based on the range of admissible values for a stance state, in agreement with [46]. Moreover, displacement boundary conditions were enforced at the distal epiphysis of femoral computational domains. In detail, nodes located at the external surface of the distal epiphysis (about 3.5 cm long) were fully restrained (8000 nodes approximately). Although there is a lack of consensus about the amplitude of the restrained region to mimic an in-vivo loading condition, clamped boundary conditions were widely adopted in well-established numerical studies that simulated the same loading condition [13, 15, 47]. As reported in [48], the amplitude of the constrained region generally affects the femoral bending state. In particular, a reduction of the clamped zone can induce an overestimation of the bending entity at the diaphysis. In the present study, the choice of the amplitude of the restrained surface has been performed in agreement with previous clinical studies available in literature [18, 19, 22, 23, 49] and it can be considered as effective to perform comparative assessments of different constitutive descriptions. The remaining boundaries were assumed stress-free.

The average computational time required for each simulation was about 15 min, by referring to a off-the-shelf desktop PC (Intel® Core™ i7-4790 CPU @3.60 GHz 8.00 GB RAM).

2.5 Parametric analyses

To investigate the effect of model parameters and the influence of constitutive description on the mechanical response of metastatic femur, an extensive parametric analysis was performed. In order to give specific comparative indications, each femur computational domain was divided in three contiguous sub-domains: (C) cortical healthy region, \(\rho _{app} >1 \, {\mathrm{g/cm}}^{3}\), labelled with the subscript c; (T) trabecular healthy region, \(\rho _{app} \le 1 \, {\mathrm{g/cm}}^{3}\) [50], labelled with the subscript t; (L) metastatic region, clinically detected as described in Sect. 2.2, labelled with the subscript l. For each sub-domain, different values of the Biot’s coefficient \(\alpha\) and pore pressure p were considered for the numerical analyses. Accordingly, we compared a purely elastic description (EL) with two poroelastic constitutive models (namely, Homogeneous Biot model, HB; Inhomogeneous Biot model, IB) accounting for bone-metastasis mechanical interaction. We highlight that the EL model, independent on p and \(\alpha\), does not account for the degradation of material properties induced by metastasis, though it represents a modelling strategy widely employed in the specialized literature [22, 23, 27]. Addressing the constitutive description IB, and in agreement with the experimental evidence reported in [51], different levels of pore pressure in the metastatic regions were considered, by comparing the influence of an increase in \(p_{l}\) with respect to \(p_{c}=p_{t}\). In particular, two different scenarios were investigated:

  • \(p_{l}\) 100% greater than \(p_{c}=p_{t}\), as experienced in the initial stage of the tumor invasion;

  • \(p_{l}\) 50% greater than \(p_{c}=p_{t}\), corresponding to a pseudo-steady tumor evolution and resulting from a time dependent mitigation of the tumoral pressure.

It is worth observing that the latter condition fulfils the reduction of pressure induced by the tumor observed after approximatively 21 days from the bone tumoral injection in the study proposed in [51].

In Table 1 the different analysed cases are summarized.

Table 1 Values adopted in the proposed parametric analyses. Biot’s coefficients (\(\alpha _{c}\), \(\alpha _{t}\), \(\alpha _{l}\)), pore pressure (\(p_{c}\), \(p_{t}\), \(p_{l}\), expressed in kPa) for healthy cortical, healthy trabecular, and metastatic regions. EL: linearly elastic model; HB: homogeneous Biot model; IB: inhomogeneous Biot model. (References:\(^\mathrm{a}\) [52, 53]\(^\mathrm{b}\) [54, 55]\(^\mathrm{c}\) [56],\(^\mathrm{d}\) [32, 57])
Full size table

2.6 Post-processing and risk indices

Numerical results obtained via FE simulations were post-processed in terms of stress and strain patterns via a home-made Matlab code (Matlab/Simulink®, v8.1, 2013a, The MathWorks Inc., Natick, USA), by identifying synthetic quantitative indications. In particular, the von Mises stress measure \(\sigma _{VM}\) and the strain energy density \(\Phi\) were locally computed and adopted as synthetic indicators of the material workload on the femoral domain. The von Mises stress measure is widely adopted in many numerical studies addressing bone tissue, for analyzing loading transfer mechanisms and/or for optimizing prosthetic devices and treatments. Nevertheless, failure mechanisms that can be generally detected with a certain accuracy via a von Mises stress assessment properly refer to symmetric failure material behaviour in tension and compression, as well as to failure states strictly dependent only on the second deviatoric stress invariant. As a result, von Mises stress measure does not provide effective indications in terms of fracture risk in materials such bone, characterized by very different limit responses in compression and tension. Besides, bone tissue usually undergoes failure mechanisms as also affected by compressibility effects and loading mode triaxiality, thereby as dependent on all the three isotropic stress invariants [58,59,60,61,62].

Accordingly, aiming at identifying a synthetic and straightforward evaluation strategy of femur fracture risk, more suitable indications can be expected by addressing the strain energy density \(\Phi\), that naturally accounts for compressibility contributions and it can be directly related to the local fracture risk level. As such, in agreement with energy-based fracture criteria inspired to the Griffith theory, crack initiation and propagation in elastic media can be locally associated with the strain energy amount, that acts as the driving force regulating the possible fracture process [63].

A local limit threshold value for \(\Phi\) cannot be straight determined adequately, since the highly patient-dependent variability of tissue strength features (locally depending on bone porosity and microstructure), as well as due to the complex interplay among possible modes that can affect fracture mechanisms. Nevertheless, just as an exemplary and rough indication, by referring to yield-stress correlations proposed in [64] and for physiological ranges of Young modulus in femurs, an indicative range for \(\Phi\) threshold level is in the order of \(0.05 \div 1 \, {\mathrm{MJ/m}}^{3}\). However, the strain energy density (always positive in sign) does not allow for distinguishing between local compression and tension states. Therefore, the local principal strains \(\epsilon _{i}(x,y,z)\) (with \({i}=1,2,3\)) were employed to define local strain measures in tension \(\epsilon ^{+}(x,y,z)\) and compression \(\epsilon ^{-}(x,y,z)\), respectively:

$$\epsilon ^{+}(x,y,z)=\max \{0,\epsilon _{1}, \epsilon _{2}, \epsilon _{3}\}$$
(7a)
$$\epsilon ^{-}(x,y,z)=-\min \{0,\epsilon _{1}, \epsilon _{2}, \epsilon _{3}\}$$
(7b)

As a result, a robust quantitative indication of tissue failure risk is introduced by locally defining the following femur risk factor (RF):

$$\mathrm{RF}(x,y,z) = \max \biggl \{\frac{\epsilon ^{+}}{\epsilon ^{+}_{o}}, \frac{\epsilon ^{-}}{\epsilon ^{-}_{o}}\biggr \} ,$$
(8)

where \(\epsilon ^{+}_{o}\) and \(\epsilon ^{-}_{o}\) denote the bone elastic limit strains in tension and compression, respectively. Accordingly, local tissue failure risk occurs in regions where \(\mathrm{RF}\ge 1\).

It is worth observing that limit strain values \(\epsilon ^{+}_{o}\) and \(\epsilon ^{-}_{o}\) should be considered as dependent on tissue local density [65], as well as on the tissue remodelling and adaptation to biomechanical and biochemical stimuli, and thereby resulting in different and evolving average limit levels in trabecular and cortical healthy/metastatic bone. Nevertheless, due to the lack of refined patient-specific experimental data, and in agreement with evidence provided in [66], a simplified assumption was considered, by enforcing constant elastic limit strain values equal to \(\epsilon ^{+}_{o}=0.73\)% and \(\epsilon ^{-}_{o}=1.04\)%.

As a matter of fact, the choice of a strain-based risk index is in agreement with experimental [67, 68] and numerical [12, 14, 15, 64, 69, 70] findings that suggest that femoral failure mechanisms are generally strain driven, and therefore they can be effectively detected by using a strain-based approach.

3 Results

In the following, in order to establish useful comparisons and indications, numerical results are presented by referring to a sub-partition of each femoral domain \(\Omega\). In particular, concerning the sketch in Fig. 3a, three sub-regions were introduced: proximal region \(\Omega _p\) (femoral segment 10 cm long); diaphyseal-proximal region \(\Omega _d^1\) (segment 12 cm long); diaphyseal-distal region \(\Omega _d^2\) (11 cm long). The bone region (about 3.5 cm long) between \(\Omega _d^2\) and the clamped segment were not accounted for post-processing analysis, to avoid possible artefacts induced by the enforced displacement-based boundary conditions. Previously-introduced three regions are used to perform a comparative analysis involving as synthetic indicators the risk factor (RF, Eq. (8)) and the strain energy density \(\Phi\). In particular, their average and peak values over each sub-domain were computed for different constitutive models and by addressing different choices of model parameters, specifically resulting in a total of 20 different numerical analyses. Moreover, results obtained in terms of the risk index RF were employed to compute the volume fraction (over each subregion) of the bone tissue prone to undergo possible failure mechanisms (i.e., where RF is greater than one).

Fig. 3
figure 3

a Partition of the CT-based computational models adopted in the post-processing phase. \(\Omega _p\): proximal region; \(\Omega _d^1\) diaphyseal-proximal region; \(\Omega _d^2\) diaphyseal-distal region. b Relative energy norm \(E_{\varepsilon }\) versus the average mesh dimension d, computed for the left femur model and for different constitutive formulations (see Table 1). Regarding the IB formulation, convergence results are limited only to the case \(p_c=p_t=50 \, \hbox {kPa}\) and \(p_l=100 \, \hbox {kPa}\), for the sake of compactness. Dashed line represents the theoretical convergence slope

Full size image

3.1 Convergence analysis

A preliminary convergence analysis was performed to set a suitable average mesh-size, to provide accurate numerical results. In detail, computational models have been defined by referring to five different mesh refinement levels, associated to the average mesh-size d varying in the range \(2.5 \div 4.5 \, \hbox {mm}\), and corresponding to a number of tetrahedral elements for each femoral model ranging from about \(7 \cdot 10^4\) (\(d=4.5 \, \hbox {m}\)) to about \(4.1\cdot 10^5\) (\(d=2.5 \, \hbox {mm}\)). To estimate the mesh convergence, the relative energy error norm \(E_{\varepsilon }\) was computed for each mesh refinement in agreement with [71]:

$$E_{\varepsilon }=\left\{ \frac{\int _\Omega \left[ {\mathbb {C}}^{-1}:(\varvec{\sigma }_d- \hat{\varvec{\sigma }} )\right] \cdot (\varvec{\sigma }_d-\hat{\varvec{\sigma }} ) \,d\Omega }{\int _ \Omega \left[ {\mathbb {C}}^{-1}:\hat{\varvec{\sigma }} \right] \cdot \hat{\varvec{\sigma }} \,d\Omega } \right\} ^{1/2}$$
(9)

where \(\varvec{\sigma }_d\) represents the total stress tensor associated to a discretization with average element size d, while \(\hat{\varvec{\sigma }}\) denotes the unknown exact solution, assumed to be effectively approximated via the finest mesh, i.e. for \(d=2.5 \, \hbox {mm}\) (\(\hat{\varvec{\sigma }}\sim {\varvec{\sigma }}_{2.5}\)).

As depicted in Fig. 3b with reference to the left femur model, the relative energy error norm decreased from \(25\%\) to less than \(4\%\) for d passing from 4.5 to 3 mm. Moreover, by analyzing the behaviour of \(E_{\varepsilon }\) versus the average mesh size, it can be noticed that the FE results reached a convergence behaviour that is close to the expected theoretical one for displacement-based quadratic FE formulation in linear elasticity [71]. As a consequence, in the following reference will be made to computational models based on \(d=3 \, \hbox {mm}\), corresponding to about \(2.5\cdot 10^5\) elements, and resulting in a good compromise between computational cost and numerical accuracy.

3.2 Risk indices

Figures 4 and 5 show average and peak values of RF and strain energy density \(\Phi\), respectively, computed for left and right femur computational models over the three regions \(\Omega _{p}\), \(\Omega _{d}^{1}\), \(\Omega _{d}^{2}\), by varying both material model formulation and porosity pressure level (see Table 1). It can be observed that HB and IB constitutive descriptions induce higher values of RF and \(\Phi\) within \(\Omega _{p}\) for all the analyzed pressure levels, compared to the values provided by the EL model. To assess the effect of the pore pressure, the percentage increases in peak values of RF and \(\Phi\) in \(\Omega _{p}\) region have been computed for each constitutive description (see Table 2). In terms of RF, in the HB model an increase of the pore pressure did not produce a further increase of RF (34% vs 35% for the left femur; 5% vs. 6% for the right femur). Conversely, looking at the IB model the combined effect of a pressure increase with the presence of a pressure gradient at the bone-metastasis interface led to a further and significant increase of RF (34% vs. 45% for the left femur; 5% vs. 9% for the right femur). In terms of \(\Phi\), in the HB model a pressure increase induced a slight \(\Phi\) increase (70% vs. 76% for the left femur; 15% vs. 18% for the right femur). Similar to RF results, in the IB model the combined effect of pressure increase with the presence of a pressure gradient at the bone-metastasis interface led to a significant increase of \(\Phi\) (70% vs 105% for the left femur; 15% vs. 30% for the right femur). Conversely, no significant difference of RF and \(\Phi\) is detected in \(\Omega _{d}^{1}\) and \(\Omega _{d}^{2}\) among predictions provided by the four different material descriptions. Moreover, in both femurs and for all the analyzed models, possible failure risk regions (namely, wherein \(\mathrm{RF}\ge 1\)) occur only within \(\Omega _{p}\) with an incidence dependent on the constitutive description and on the metastasis features. In detail, referring to the left femur, HB and IB models predict wider risk region (about 28% in volume of \(\Omega _{p}\)) in comparison to EL model (18% of \(\Omega _{p}\)). On the contrary, in the right femur, where metastases slightly affect the cortical bone, any significant difference in the extension of the risk region was observed by comparing the different models, such a region occupying about 1% of \(\Omega _{p}\).

Fig. 4
figure 4

Average (bars) and peak (lines) values of risk factor (RF) computed for left and right femoral models over three different subregions (see Fig. 3a). EL: linearly elastic model (dark gray); HB: homogeneous Biot model (light gray); IB1: inhomogeneous Biot model with \(p_{l}=1.5p_{c}\) , with \(p_{c}=p_{t}\) (black); IB2: inhomogeneous Biot model with \(p_{l}=2p_{c}\) with \(p_{c}=p_{t}\) (white). Dashed red line indicates the risk threshold \(\mathrm{RF}=1\)

Full size image
Fig. 5
figure 5

Average (bars) and peak (lines) values of strain energy density (\(\Phi\)) computed for left and right femoral models over three different subregions (see Fig. 3a). EL: linearly elastic model (dark gray); HB: homogeneous Biot model (light gray); IB1: inhomogeneous Biot model with \(p_{l}=1.5p_{c}\) , with \(p_{c}=p_{t}\) (black); IB2: inhomogeneous Biot model with \(p_{l}=2p_{c}\) with \(p_{c}=p_{t}\) (white)

Full size image
Table 2 Increases in peak values of RF and \(\Phi\) (with respect to the linear elastic description EL), and extension of failed region (% in volume of \(\Omega _{p}\)), for each constitutive description and computed for left and right femoral models
Full size table

For both femurs, the increase of the pore pressure slightly affect the extension of the failure risk region, as shown in Fig. 6.

Fig. 6
figure 6

Failure risk regions (RF\(\ge\)1) for the left (a) and right (b) femurs. Comparisons between HB and IB models, for different levels of pore pressure (notation is the same as in Table 1)

Full size image

The influence of CT-based geometry and metastasis localization is further highlighted in Figs. 7 and 8, where the spatial distributions of RF and \(\Phi\) within \(\Omega _{p}\) are shown, respectively, for the left and right femur models. The differences that can be straight highlighted are essentially due to the spatial occurrence of metastases and their constitutive description. In particular, the presence of a metastatic region in the greater trochanter of the right femur (see Fig. 1a) tends to induce a strain-energy-density localization detected only by a Biot-based approach, able to take into account the inhomogeneous distribution of pore pressure between metastatic and healthy bone tissues (Fig. 8b).

Fig. 7
figure 7

Spatial distribution of the RF index computed within \(\Omega _{p}\) for a left and b right femoral models comparing four different material descriptions. Notation is the same as in Fig. 4. (Circled) Fracture risk index contours at femur head and (boxed) at the coronal section \(y= 150\,\mathrm{mm}\)

Full size image
Fig. 8
figure 8

Spatial distribution of the strain energy density computed within \(\Omega _{p}\) for b left and b right femoral models comparing four different material descriptions. Notation is the same as in Fig. 4. (Circled) Strain energy contours at femur head and (boxed) at the coronal section \(y= 150\,\mathrm{mm}\)

Full size image

For the sake of completeness, Fig. 9 depicts the spatial distribution of Von Mises stress, clearly highlighting different stress patterns for the left and right femur models.

Fig. 9
figure 9

Spatial distribution of the von Mises stress (\(\sigma _{VM}\)) computed within \(\Omega _{p}\) for a left and b right femoral models comparing four different material descriptions. Notation is the same as in Fig. 4. (Circled) \(\sigma _{VM}\) contours at femur head and (boxed) at the coronal section \(y= 150\,\mathrm{mm}\)

Full size image

4 Discussion

Numerical results clearly showed that a refined geometric-constitutive description of bone-metastasis interaction, based on a heterogeneous biphasic poroelastic formulation, allows identifying strain localization mechanisms, related to possible fracture risk and not ultimately detected by a purely-elastic description, and significantly impacts the RF estimates. Even if it is not possible to conclude which model among HB and IB descriptions is more effective, to the authors opinion this does not affect the main conclusion of the work, i.e. that a refined poroelastic description accounting for an interaction between bone and metastasis is able to localize local mechanical signals that are not completely detected by a linearly elastic model and that are strictly connected to biochemical processes driving the onset and the evolution of metastasis.

The results obtained in terms of RF for the left femur outline its critical mechanical behavior in agreement with the clinical evidence, and thereby confirm the need of the surgical treatment. This occurrence is in accord with the Mirels scoring system (Mirels score clinically computed as equal to 11). On the other hand, the Mirels score clinically attributed to the right femur (equal to 8, and thereby judged as a “dilemma”, but resulting in the specific case in a conservative management) does not reflect the dangerous fracture risk level numerically experienced. Accordingly, proposed result furnishes clear evidence of the potential of CT-based FE modeling in fracture risk assessment of metastatic femurs and of its ability to overcome severe limitations associated with the MSS.

Proposed results also highlight that bone regions mainly prone to fracture occur at the proximal third of the femur, with localization mechanisms depending on the distribution of the metastatic lesions. The presence of more extensive metastatic zones also affecting the cortical bone (as in the case of the left femur for the analyzed case study) seems to induce a transition from a transcervical neck-fracture mechanism (detected for the right femur model) to an intertrochanteric one. It is worth observing that such a piece of evidence, clearly highlighted by analyzing Figs. 7 and 8, cannot be described entirely by referring to a sole stress analysis based on Von Mises equivalent stress measure. Von-Mises-based analyses (widely adopted for investigating bone mechanical performance), in fact, can provide useful insights into the overall load transfer mechanisms. However, they are not able to identify the influence of the constitutive models on the detection of localization effects for the onset and propagation of fracture due to the sharply different tension-compression bone strength limits. On the contrary, the proposed strain-based measures numerically confirmed that accounting for localizations in strain energy density and strain measure, can be more effective for assessing local fracture risks.

5 Conclusions

In the current work, a novel subject-specific CT-based FE modeling strategy was developed to investigate the biomechanical behavior of metastatic femurs. In particular, the study integrates a specific constitutive description of metastasis into a CT-based FE computational model of femurs, by investigating the effects on the macroscopic mechanical response depending on the material parameters. The developed strategy was applied to a clinical case with both left and right femurs affected by multiple metastases. Computational domains were reconstructed from CT images and analyzed under a pseudo-functional compressive load. The proposed approach predicts higher failure risk levels with respect to a formulation that neglects a specific constitutive description of metastasis, with different bone zones prone to failure. Moreover, the developed strategy gives a more precise indication on the fracture risk also in those cases for which the Mirels scoring system (MSS) does not give a support in clinical decision making. The obtained numerical results infer on the importance of a refined description of bone metastasis, enabling to describe the bone-metastasis interaction by accounting for both heterogeneous material properties and the different influence of fluid pore pressure on healthy and metastatic bone matrix. The approach opens towards the possibility of a predictive identification of onset and propagation of fracture mechanisms in femurs possibly affected by multiple metastases, aiming at filling the gap concerning the altered material properties due to the presence of metastatic lesions [8, 22].

We conclude mentioning limitations and perspectives of the present study. One of the significant limitations is related to the clinical study design. Only one patient has been herein investigated. As such, with the aim to further substantiate the reliability of the proposed strategy and perform a more effective comparison among HB and IB models, it will be necessary to apply the proposed approach on a broader retrospective clinical cohort and on cadaveric metastatic femurs for which the experimental measurements and observations are available. Other limitations arise from the CT-based FE modeling procedure adopted, as mainly related to the lack of a phantom for CT calibration and to the consideration of a single loading condition (i.e., a maximum compressive force on the femoral head) to investigate the femoral mechanical response. In particular, regarding the latter aspect, the large variability in the direction of the forces acting on the femur during physiological conditions should be considered to enhance the detection of localised weak features, as demonstrated for instance in [18, 19]. Moreover, in addition to physiological conditions, accidental conditions (e.g., a fall) should be considered in forthcoming analyses, due to the evidence that the majority of fractures occurs during a fall [72, 73]. Finally, we foresee to investigate further: (i) the effect of irregular shape of metastasis on the mechanical response of femur with respect to a regular ellipsoidal shape; (ii) the possibility of employing different imaging modalities, e.g. Magnetic Resonance, to improve the detection of metastatic boundaries; (iii) the accuracy of the adopted constitutive description, by referring to specific measurements on ex-vivo metastatic samples; (iv) the influence of possible remodelling processes on the localization effects of failure risks.