1. Imaging Techniques for Multiscale Damage Assessment and Prediction

The use of imaging techniques enables researchers to understand bone damage at different hierarchical scales. It is particularly relevant in the comprehension of the implications of fracture processes in the deterioration of bone quality. This section presents an overview of the available imaging techniques (Figure 1) to visualize bone morphologies, to assess bone fractures and to predict the fracture risk from the macro-scale to the nano-scale. The macro-architecture is currently evaluated by means of common clinical images. At lower hierarchical levels, the identification of damage processes is more complex and requires higher-resolution techniques.

Figure 1. An overview of the main imaging techniques to assess bone damage at different scales. Macro- and meso-scale techniques are depicted in a darker color, while micro- and nano-scale imaging techniques are represented in a brighter color. VFA: Vertebral Fracture Assessment; QCT: Quantitative Computed Tomography; MRI: Magnetic Resonance Imaging; pQCT: peripheral Quantitative Computed Tomography; Micro-CT: Micro Computed Tomography; LSCM: Laser Scanning Confocal Microscopy; SEM: Scanning Electron Microscopy; SR imaging: Synchrotron Radiation imaging; AFM: Atomic Force Microscopy.

Different techniques are compared in terms of outcomes, in vitro or in vivo applications, resolution, two- or three-dimensional features and the main advantages and disadvantages.

2. Bone Damage Physical Principle

Imaging techniques are a valuable tool for the visualization of bone features and multiscale damage and could help in the comprehension of bone damage mechanisms. Due to the fact that bone is a hierarchical material with multiple linked scales, bone fracture mechanisms also occur at the multiscale. Those mechanisms are particularly complex and not completely understood.

At the macro-scale, there are essentially two main bone fracture modalities: impact fracture and stress fracture. The former is a fracture caused by an overload, while the latter is associated with cyclic loading and it causes material failure by excessive damage accumulation ^[30][31]. Stress fractures include both fatigue fractures, which result from strong physical activity, and fragility fractures, which are consequences of everyday common activities. The latter are often linked to age-related decreases in bone’s ability to self-repair. This often happens in age-related bone pathologies such as osteoporosis. At the macro-scale, the shape, size and density of bone strongly affect fracture behavior ^[29].

At the meso-scale, trabecular bone microarchitecture influences the damage. Trabecular bone is subjected to a variety of loads during activities of daily life and the orientation of the applied load with respect to the trabecular distribution plays an important role in a possible trabecular tissue degeneration. Indeed, trabeculae generate a three-dimensional, open porous space. Pathologies such as osteoporosis may lead to a conversion from a plate-like to a rod-like trabecular morphology, which contributes to increasing the fracture risk ^[32]. Typically, trabecular failure initiates at the weakest trabecula or at the weakest trabecular region. From there on, the failure will progress and failure bands will develop. During a tensile test, there will be only one failure band; multiple failure bands, however, occur in compressive tests ^[33]. Gibson ^[34] investigated how the basic failure modes (bending, compression, tension, shear) can explain trabecular bone failure (Figure 2). It has been also observed that complete fractures are present in trabeculae that are oriented transversally to the direction of the applied load ^[35].

Figure 2. Different meso-scale failure modes of a trabecular network. (a) An unloaded trabecular cell. (b) Brittle crushing. (c) Bending of horizontal rods, compression of vertical rods. (d) Buckling of vertical rods, bending of horizontal rods. (e) Plastic yielding. Readapted from ^[34].

At the micro-scale, a complete comprehension of the bone damage physical principle is still lacking. Bone is able to sustain many micro-cracks, if they are not critical; this suggests that bone can be classified as a “damage-tolerant material” ^[36]. This concept is particularly interesting when applied to bone microstructure. At the micro-scale, the damage assumes two frequent morphologies: linear and diffuse damage. The linear damage appears as a sharp line of 100–200 µm ^[37] separating the bone matrix, and it is typical of regions subjected to compression ^[30]. When external compression boundary conditions are applied, linear cracks propagate parallel to the compression loading axis ^[38]. On the contrary, diffuse damage is characterized by a cloud of tiny cracks ^[39], whose dimensions are less than 1 µm ^[37], and this is typical of regions subjected to tension. Diffuse micro-cracks are distributed transversally to the applied tensile stress ^[40]. Voide et al ^[41] suggest that diffuse damage could be considered as linear damage when oriented differently with respect to a global reference system. Micro-cracks occur on multiple planes when subjected to torsion or mixed-mode loading (Figure 3).

Figure 3. Orientation of micro-cracks in a bone sample subjected to different loading conditions.

Micro-cracks usually initiate at regions characterized by high stresses. As soon as the crack initiates, its growth behavior is influenced by micro-scale heterogeneities. In this context, particular interest has recently been devoted to cellular-level porosity represented by the lacunar network, whose role is still under debate. A variation in the lacunar distribution or in the lacunar shape, as happens with aging, significantly affects the bone resistance to fracture ^[42]. Lacunae play a dual mechanical role, having an effect on both strength and toughness. In the first instance, lacunae are stress concentrators ^[43], by representing bone discontinuities able to locally amplify stresses and strains. The average strain around lacunae is 1.5–4.5 times higher than the remote strain applied to the surrounding tissue ^[44][45]. In this sense, lacunar system should contribute to strength decreases.

However, considering bone as a damage-tolerant material, experimental investigations ^[41] show that lacunae are not the starting point for the micro-cracks that ultimately lead to bone failure. Micro-cracks nucleate generally at canals ^[41] or at cement lines and inter-lamellar zones ^[46][47], where the stress amplification is greater than at the lacunae ^[48][49]. According to those investigations, lacunae make a beneficial contribution to toughness ^[38]. Voide et al. suggest that lacunae exert an attraction upon the existing micro-crack: the deviation of the crack path reduces the energy of its progression, slowing down crack propagation ^[50]. The two effects of the lacunar network are schematized in Figure 4. The role of the lacunar network still needs to be clarified by an effective experimental validation of the proposed hypotheses.

Figure 4. The dual effect of the lacunar system on microdamage: lacunae as sites for crack initiation (a) and lacunae as micro-crack deviators (b); (b) is adapted from ^[38].

Finally, at the nano-scale, bone failure and fracture are influenced by bone’s composition in terms of hydroxyapatite and collagen fibers. Hydroxyapatite minerals show a brittle behavior and are more resistant to compression, while collagen fibers are more resistant to tension. When the crack propagates, the hydroxyapatite will fail first and the complete failure only occurs when collagen fibers are fully stretched ^[51].

It is understood how the comprehension of the bone damage physical principle at the multiscale level is a crucial point for the implementation of computational damage models. Those models have recently attracted a high degree of interest ^[52], due to the fact that patient-specific simulations, in particular in this context of age-related bone fractures, allow a more accurate prognosis ^[53][54]. The ability to quantify the effects of aging and pathologies such as osteoporosis is a challenging issue that computational damage models should consider.

3. Multiscale Computational Damage Models

Given that macro- and micro-architecture play an important role in bone damage and crack propagation, computational analyses can be interesting tools to investigate the effect of their variation due to aging on bone fractures. Using finite element (FE) models, generated directly from micro-CT or synchrotron images, it is possible to perform a “virtual experiment”, able to simulate mechanical loading of the sample and to evaluate damage formation. The main advantage of FE models is that, once the geometry of the structure is implemented, the same model can be used for multiple analyses, i.e., applying different loading conditions that may mimic everyday life loadings (e.g., walking, running).

Furthermore, bone FE models provide new insight into the relationship between damage propagation and macro- and micro-scale architecture, by allowing the localization of zones that are more prone to fracture. While in the past FE models were limited by input image resolution and computational power ^[53], nowadays, accurate simulations are feasible, thanks to the recent developments in imaging and in optimized computational software. The increased accuracy of bone FE damage models demonstrates that computational simulations are a reliable way to improve the choice of more effective experimental tests, potentially leading to a reduction in the number of performed tests ^[55][56].

In this section, both macro- and micro-scale models are presented, with an attempt to investigate bone fracture initiation and propagation at different scales ^[57].

Macro-scale organ-level models primarily focus on the effects of mechanical stimuli on bone resistance to fracture. Those models, that start mainly from CT reconstructions, try to implement the observation of J. Wolff ^[58], related to bone’s capability to adapt to the external mechanical loading conditions ^[59]. At this scale, micro-architectural features are neglected. The outer shape object of the problem is simplified, and filled with different meshing strategies, for instance, geometry-based with tetrahedral elements, or voxel-based with hexahedral elements. While the first method, that requires a smaller number of degrees of freedom, is not fully automated, the voxel-based hexahedron meshing has no geometrical limitations and it is a fast and completely automated mesh generation technique ^[60]. However, in the second case, outer surfaces or sharp geometrical discontinuities may generate a lack of convergence. Generally, the voxel-based meshing strategy allows a simpler and more effective interface with CT scans. The bone material is typically considered as a continuum solid ^[61]. The model is loaded with defined boundary conditions (compression, torsion or mixed-mode loading). The external load produces local stresses and strains in each considered element of the mesh (Figure 5), which could be correlated by means of the constitutive equation ^[62]. The constitutive equation, in the condition of linear elasticity, provides a linear relation between stresses and strains by means of the mechanical properties of the bone tissue, assuming the strains to be small or infinitesimal ^[63].

Figure 5. Starting from a human bone (e.g., femoral head), by means of CT acquisition, it is possible to obtain a macro-scale finite element (FE) model of the sample. (a) Geometry-based FE models with tetrahedral elements. (b) Voxel-based FE model with hexahedral elements.

The initiation and propagation of damage in the macro-scale bone model is defined according to a proper failure criterion ^[64]. Based on the theory of elasticity, those models generally assume that bone resorption occurs when local mechanical stress overcomes a homeostatic stress state ^[65]. In the definition of failure at the macro-scale, it should be considered that local tissue failure could provide a loss in the structural integrity. The main results of macro-scale models, as reported in the scheme of Figure 6, show that subject-specific FE models predict values of strain with an accuracy of 90%, obtained by comparison with experimental measurements on cadaver bones (the average errors on surface strains are lower than 10%) ^[66]. The encouraging results are due to the coexistence of a precise geometric reconstruction, an appropriate choice of the density–elasticity relationship, an accurate application of boundary conditions and a suitable algorithm for material property evaluation.

Figure 6. Scheme for the implementation of a whole-bone model. Results, potentialities and limitations are highlighted for each step. EMG: Electromyography.

A relevant aspect of macro-scale models is that they can consider an accurate map of the tissue elastic modulus distribution, based on the density distribution over the continuum bone model, as performed by Viceconti ^[64] and Taddei ^[67]. These models show low average error in predicting failure load in vitro and have a high precision in identifying the location of failure ^[64].

Further attempts in macro-scale FE modes have been performed in order to capture the non-linear behavior of bone that depends both on the anatomic site and the loading mode. In particular, Imai et al. ^[68] observe that the prediction of vertebral fracture is intricate due to the complex geometry of the vertebra and its elastoplasticity. They implement non-linear FE models of the whole vertebra; in order to consider bone heterogeneity, the mechanical properties of each element are computed from the Hounsfield unit value. The correlation between the measured value of fracture strength and the predicted values with the non-linear FE models shows significant improvements ^[68] with respect to previous linear and simplified FE models. Other attempts in modeling bone failure have been performed by Harrison et al. ^[69], considering that the tissue failure consists of two phases: damage and fracture. This study develops a computational model consisting of an explicit representation of complete failure, incorporating non-linear damage criteria, fracture criteria, cohesive forces, asymmetry and large deformation capabilities.

However, those models are not able to capture patient-specific cortical nor trabecular micro-architecture, which play a crucial role in bone strength ^[69][70][71].

Micro-scale models provide further insights into bone fracture prediction at smaller scales. Traditionally, micro-architectural features were assessed by means of histological sections ^[72] that lack in three-dimensionality. Nowadays, the input of micro-scale models often comes from micro-CT or SR images, that have a resolution able to non-destructively capture not only trabecular architecture, but even ultrastructural porosities, such as lacunae. However, this corresponds to higher computational costs, calculated by van Rietbergen et al. ^[73] through an element-by-element method that implies the use of uniformly shaped elements to reduce memory allocation and optimize computational times. Some issues arise when dealing with micro-CT or SR imaging techniques: filtration and segmentation. Reconstructed image data include noise that should be removed or at least reduced by filtering. The choice of an adequate filter (the Gaussian filter is the most suitable one, which is easily implementable and fast in computation ^[74]) is essential in order to obtain an input image for the model, as close as possible to the original sample. Another relevant aspect is the correct segmentation process that selects those voxels that are below or above a defined threshold, so as to separate those elements that are bone and those that can be considered as voids. This distinction is particularly relevant when dealing with the identification of micro-cracks and lacunae, whose typical size is significantly smaller (about 10 µm ^[38]) with respect to the trabecular architecture. Micro-FE models are often used for the prediction of bone fracture. For this purpose, specific damage criteria are implemented. In the work of Pistoia et al. ^[75], a micro-FE linear analysis is performed on the human distal radius. In order to predict bone failure, the chosen criterion is that the bone failure initiation occurs when bone is strained beyond a critical value, defined as the yield strain. This model is able to give a better prediction of bone failure loads (R² = 0.75 correlation with experimental testing) with respect to macro-scale investigation techniques such as DXA measurements. Linear micro-FE models are also proposed for the study of the interaction between micro-cracks and micro-porosities, such as lacunae or surface discontinuities in bones. The role of the lacunar network in damage initiation and propagation is still unclear, due to the difficulty in obtaining dynamic visualization of the advancing crack. Micro-FE models shed some light on micro damage initiation and propagation. In this context, it is necessary to define a proper descriptor of the stress distribution, such as the strain energy density (SED) or the maximum principal stress (σ₁) criteria. SED is a good predictor for the mechanical environment sensed by the osteocytes that reside in lacunae ^[76] and σ₁ is a mechanical quantity often used to assess the failure in brittle materials ^[77]. Those two criteria show different sites of crack initiation and different directions for crack propagation (Figure 7).

Figure 7. Schematic output of models of a bone sample with a lacuna subjected to compression. Crack initiation (dots) and crack propagation (lines) are highlighted, according to strain energy density (SED) damage criterion and to σ₁ criterion.

The σ₁ criterion locates peak stresses and strains at regions matching with the experimental initiation stresses and strains found by Voide et al. ^[41]. Donaldson et al. ^[77] compare different linear micro-FE models that aim at the identification of crack initiation sites and of the role of the lacunar network. They present a relevant novelty by implementing a stress gradient model, schematized in Figure 8, that shows a clear directionality in the advancing crack (instead of the less realistic stress limit algorithm, that is able to simulate only accumulated damage around voids or surface discontinuities).

Figure 8. Stress gradient algorithm. “Critical voxels” are the voxels that exceed a defined stress threshold. “Tip voxels” are the center of the spherical sensation region of the osteocytes and are eligible sites for damage propagation. In the stress gradient algorithm, a key role is played by the deletion direction vector (DDV), that indicates a direction for crack propagation. The iteration stops when the total failure of the sample is reached, according to the Pistoia criterion ^[75].

The stress gradient algorithm clarifies that damage starts from surface discontinuities or blood vessel (Figure 9) and not from lacunae, confirming from a computational point of view the hypothesis.

Figure 9. Stress gradient algorithm (on SR CT data from a cortical bone sample) results: damage starts from big porosities and its propagation follows the stress gradient. Micro-damage is deviated by the presence of lacunae. Readapted from ^[77].

Despite the explained potentialities of linear micro-FE damage models, micro-damage predictions could be improved by implementing the bone non-linear behavior. The experimental tests focused on non-linear bone properties indicate ductile failure modes ^[78], in which damage is combined with the plasticity component of collagen fibers ^[79]. Hammond et al. ^[80] perform non-linear FE models (in this case, material non-linearity is considered) on trabecular samples from the distal femur of a human cadaver. They consider two micro-crack formation criteria, one for isotropic tissue models and one for anisotropic tissue models. For the isotropic tissue models, micro-cracks initiate if the ratio between maximum tensile principal stress and micro-crack initiation strength is equal to one. For anisotropic tissue models, micro-damage initiates if the maximum ratio between the normal of shear tractions with respect to the micro-crack initiation strength is equal to one. This non-linear micro-FE model is able to demonstrate that the anisotropy of bone tissue significantly contributes to bone fracture resistance ^[81]. Stipsitz et al. ^[82] try to implement an efficient micro-FE solver for large-scale non-linear simulations. However, the non-linearity of the system makes the results highly dependent on small structural deviations introduced by coarsening the structure. Additionally, one of the main drawbacks of the use of non-linear models is the increased computational cost, almost 10 times higher than linear simulations ^[53].

An increased degree of complexity can be found when analyzing nano-scale damage models that require further understanding of the variations in chemical pathways due to aging or disease. The first attempts in nano-scale bone fracture modeling were performed by Dubey et al. ^[83]: in their study, bone tissue fracture properties are based on the atomistic strength analyses of type I collagen in combination with hydroxyapatite interfacial arrangement, using molecular dynamics. Additional nano-scale studies ^[84] reveal that a decrease in the hydroxyapatite crystal size may change the mechanical behavior of the whole bone and the failure mode: from brittle-like crack-driven failure in larger crystals to a more widespread failure mode in the smaller ones. However, heterogeneities at the nano-scale remain difficult to model in a constitutive law and their role and effects are still unknown ^[85]. The impossibility to perform direct experimental tests on bone samples at the nano-scale is also a reason why those models are less considered for the purpose of this review.

To summarize, the increased image resolution and the optimization of parallel computational architecture for the solving of FE damage models seem promising tools for the improvement of the clinical understanding of fracture risk prediction. It is necessary to point out that some attempts in performing multiscale analyses have been performed by implementing bottom-up multiscale approaches, starting from the nano-scale up to the meso-scale ^[86] and analogously by Kwon et al. ^[87], who apply the multiscale analysis with variable geometrical parameters in order to determine its effect on the bone properties.

Before translating in silico damage models to clinical practice, the essential step is the experimental validation of the mentioned models that are still in their infancy, especially at smaller scales.