Specimen-specific finite element prediction of surface strain at the distal radius in a fall configuration

Keywords:

• Forward fall
• bone strength
• distal radius
• Non-axial loading
• finite element model

1. Introduction

Distal radius fractures are the most common upper extremity fractures in adults (65 years and older) occurring after a fall from standing height or lower (Vogt et al. 2002). These fractures are associated with age related decline in bone quality and because they occur earlier in life than other osteoporotic fractures, they can be interpreted as a warning signal for later and more deleterious fractures (Melton et al. 2010).

Subject-specific finite elements models have been proposed to evaluate both bone strain and fracture risk. The accuracy of these models is heavily dependent on material properties, boundary conditions and the chosen failure criterion that defines fracture threshold (Edwards and Troy 2012). The main aim of the current study was to assess the surface strains prediction of radius specimen-specific finite element models.

2. Methods

2.1. Specimens and image acquisition

Data were acquired from thirty fresh-frozen cadaveric left radii from elderly donors (50–96 y.o. 15 males, 15 females). The most distal 8 cm of the bones cleaned of soft tissues were imaged with a clinical CT scan (New-Tom 5 G, QR, Verona, Italy, voxel size: 150 µm × 150 µm × 150 µm) along with a forearm calibration phantom (QRM, Moehrendorf, Germany) with calcium hydroxyapatite equivalent concentrations of 0 and 200 mg/cm³. The scans of the thirty specimens were used to establish the following linear relationship between CT Hounsfield units (HU) and calcium hydroxyapatite equivalent density ${\rho }_{HA}$ in g/cm³: $${\rho }_{HA}=0.00038\times HU-0.1262$$

2.2. Ex vivo protocol

The experimental protocol used here has been previously described in details by Zapata et al. (Zapata et al. 2017). Briefly, the radii were potted in a polyurethane resin in a steel cylinder, with an angle of approximately 75° between the anterior face of the radius and the impactor. The specimens were then placed in a horizontal cylinder bar on a rail system which is free to slide along the loading axis. This bar weights 12.5 kg, which is an arbitrary value representing the mass involved in a fall (i.e. a percentage of the body weight). Speckle was painted on the specimens. They were then loaded through a resin mold at 2 m/s using a hydraulic high-speed testing machine (LF technologies, France) recorded by four high-speed cameras (FASTCAM SA3, Photron, Japan) and leading to fractured and non-fractured radii. Video recordings were analyzed using VIC3D stereo-correlation software (Correlated Solution, South Carolina, USA) in order to compute Von Mises surface strain fields on a specified non-fractured area by means of principal strains.

2.3. Finite element models

Models of the radii based on segmented CT data from the open-source platform 3 D Slicer were imported into Ansys Software (v. 2019 R1 Ansys^® Inc., USA) to be meshed. The models consisted of 200 000 10-nodes tetrahedral elements in accordance with a preliminary mesh convergence analysis resulting in an approximate 1 mm edge length.

Because no density–elasticity relationships were established especially for the radius, we investigated three femoral based equations (Eqs. (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) –(3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) ) to assign inhomogeneous linearly isotropic material properties to the finite element models: (1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) (Keller 1994), (1) (2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) (Morgan et al. 2003), (2) (3) $E=10\mathrm{ }095{\rho }_{HA}$(3) (Duchemin et al. 2008), (3) where Young’s modulus E is expressed in MPa, and the densities, ${\rho }_{\mathit{ash}}$ (ash density), ${\rho }_{\mathit{app}}\mathrm{ }$(apparent density) and ${\rho }_{HA}\mathrm{ }$(hydroxyapatite equivalent density) in g/cm³. For Eq (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) , calcium hydroxyapatite equivalent density was converted to ${\rho }_{\mathit{ash}}$ using: ${\rho }_{\mathit{ash}}=0.0698+0.839\mathrm{ }{\rho }_{HA}$ (Kaneko et al. 2004).

For equation (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) ${\rho }_{\mathit{ash}}$ was divided by 0.6 to obtain ${\rho }_{\mathit{app}}$ (Schileo et al. 2008). This resulted in around 220 materials ranging from 0.05 to 12.1 GPa for Eq. (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) , from 0.03 to 12.0 GPa for Eq. (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) and from 0.01 to 15.7 GPa for Eq. (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) . Each element was assigned a Poisson’s ratio of 0.3.

The models were specifically oriented in the same configuration as the experiments (with an uncertainty of ±2°) in order to apply the maximal peak load observed experimentally along the z axis on the articular surface and to fully constrain the proximal end of the radius at the location of potting. Analyses were performed using Ansys (v. 2019 R1, Ansys^® Inc., USA). Average values of the surface strain field (Von Mises strains) were computed in order to compare the bone numerical response with the experimental data analyzed with stereo-correlation.

2.4. Data analysis

Model predicted and experimentally measured means Von Mises strains were compared using Spearman correlation (ρ), root mean squared error (RMSE), and maximum error (Max err). Scatter was assessed using Bland–Altman plots.

3. Results and discussion

Due to experimental issues, only 27 bones were tested successfully. Sixteen of them were not fractured and were considered for the model evaluation in the current study. Qualitatively, similar results were found in von Mises spatial strain distribution for the three equations. Furthermore, patterns present common characteristics with the experimental fields (Figure 1).

Figure 1. Comparison between experimental (A) and computational (B: Eq. (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) , C; Eq. (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) , D: Eq. (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) ) Strain distribution for a non-fractured specimen.

Nonetheless, the finite element predicted strains vary as a function of Eqs (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) –(3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) (Figure 2). Similar results were found for Eqs (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) and (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) predicting more accurately the experimental values for non-fractured bones whereas Eq. (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) overestimate the measurements. Considering the error due to the ±2° uncertainty, correlation coefficients for experimentally measured strains versus predicted strains ranged from ρ = [0.68-0.71] (p < 0.01) for Eq. (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) to ρ = [0.75-0.82] (p < 0.001) for Eq. (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) and to ρ = [0.78-0.79] for Eq. (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) (Table 1). The RMSE and Max Error were smallest for Eqs (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) and (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) and largest for Eq (1) (Table 1). Moreover, Bland-Altman plots showed a randomly distributed scatter across strain magnitudes for Eqs (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) and (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) with the difference obtained between best predicted and experimental values within the mean of difference ±1.96SD (respectively: mean difference: 0.064%ε and 0.057%ε, ±1.96SD: [-0.12, 0.25]%ε and [-0.14, 0.25]%ε).

Figure 2. Predicted (with an uncertainty of ±2°) and measured mean von Mises strains for each of non-fractured specimen classified from the smallest experimental measured strain to the largest.

Table 1. Statistical parameters.

	Eq (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1)	Eq (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2)	Eq (3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3)
Spearman ρ	[0.68-0.71]*	[0.75-0.82]**	[0.78-0.79]**
RMSE (%ε)	[0.65-0.98]	[0.14-0.19]	[0.15-0.18]
Max err (%ε)	1.63	0.37	0.36

* p-value < 0.01,

** p-value < 0.001.

Overall, predicted strains using the Eqs. (2)(2) $E=6\mathrm{ }850{\rho }_{\mathit{app}}^{1.49}$(2) –(3)(3) $E=10\mathrm{ }095{\rho }_{HA}$(3) were statistically similar and most closely matched to measured strains in terms of correlation coefficients, error, and scatter. The difference with Eq. (1)(1) $E=10\mathrm{ }500{\rho }_{\mathit{ash}}^{2.29}$(1) may be explained by the global lower Young’s modulus found for most materials due to the use of ash density and leading to an overestimation of the von Mises strains.

4. Conclusions

Radius bone strain can be predicted by finite element analysis in forward fall configuration. After the validation in the non-fractured case (elastic), the next step will be the investigation of a fracture strain-based criterion.

Disclosure statement

No potential conflict of interest was reported by the authors.