Prediction of the Inelastic Behaviour of Radius Segments: Damage-based Nonlinear Micro Finite Element Simulation vs Pistoia Criterion Biomechanics

The Pistoia criterion (PC) is widely used to estimate the failure load of distal radius segments based on linear micro Finite Element ( l FE) analyses. The advantage of the PC is that a simple strain-threshold and a tissue volume fraction can be used to predict failure properties. In this study, the PC is compared to materially nonlinear l FE analyses, where the bone tissue is modelled as an elastic, damageable material. The goal was to investigate for which outcomes the PC is sufﬁcient and when a nonlinear (NL) sim- ulation is required. Three types of simulation results were compared: (1) prediction of the failure load, (2) load sharing of cortical and trabecular regions, and (3) distribution of local damaged/overstrained tissue at the maximum sustainable load. The failure load obtained experimentally could be predicted well with both the PC and the NL simulations using linear regression. Although the PC strongly overestimated the failure load, it was sufﬁcient to predict adequately normalized apparent results. An optimised PC (oPC) was proposed which uses experimental data to calibrate the individual volume of overstrained tissue. The main areas of local over-straining predicted by the oPC were the same as estimated by the NL simulation, although the oPC predicted more diffuse regions. However, the oPC relied on an individual cal- ibration requiring the experimental failure load while the NL simulation required no a priori knowledge of the experimental failure load. (cid:1) 2021 The Authors. Published by Elsevier


Introduction
Nonlinear micro Finite Element (lFE) simulations of bone can lead to unique insights into the failure behaviour of bone. Based on high-resolution micro computed tomography (lCT) scans the over-straining of bone tissue can be estimated through-out the structure on the level of individual trabeculae. However, local insights at comparable high resolutions are not yet available experimentally.
Nonlinear modelling of bone tissue is essential to directly obtain the failure behaviour from simulations. The increased availability of high performance computing (HPC) resources and parallel solvers enabled a broader community to include material nonlinearity when modelling bone tissue on the microscopic scale (Zhou et al., 2016;Hosseini et al., 2017;Arias-Moreno et al., 2019). Nevertheless, the computational requirements for a nonlinear simulation are hundred-fold larger than for a linear-elastic simulation and necessitate special-purpose solvers (Fields et al., 2012;Stipsitz et al., 2019). However, so far only small improvements in the prediction of the maximum sustainable loads were reported compared to a linear-elastic material behaviour (MacNeil and Boyd, 2008;Christen et al., 2013). Thus, it is important to know for which objectives nonlinear material behaviour is required and when a linear-elastic constitutive model is sufficient.
Besides studying the predictability of the maximum sustainable load, lFE is mainly applied for basic research questions where local estimates are required. For instance, the internal stress distribution obtained from lFE simulations has been used to study the role of the cortical-trabecular load sharing (Ulrich et al., 1999;Pistoia et al., 2003;Eswaran et al., 2006;Johnson and Troy, 2018). In most of these lFE studies a linear-elastic material model was applied without investigating if a nonlinear tissue behaviour would alter the results. Radius segments are a good choice to study the benefit and limitations of a nonlinear material model: As a comparison to nonlinear simulation results the failure load can be estimated using the Pistoia criterion (PC) based on linear-elastic lFE simulations (Pistoia et al., 2002;Varga et al., 2010;Mueller et al., 2011). Radius segments allow direct comparisons of experimental and simulated https://doi.org/10.1016/j.jbiomech.2020.110205 0021-9290/Ó 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). data since the experimental set-up and the boundary conditions for lFE models leading to Colles' fractures are largely standardized (van Rietbergen and Ito, 2015).
In this study, a damage-based material model was compared to a linear-elastic model for compression tests of radius segments. A damage-based material model was chosen because of the availability of an efficient implementation which permits the large-scale simulations required for this study (Stipsitz et al., 2019). To summarize, the purpose of this study is to determine when it is appropriate to utilize a nonlinear FE approach versus the linear Pistoia criterion. We hope that this study will help other biomechanical researchers to decide if a nonlinear analysis is required, especially considering the hugely increased computational costs.
Three questions were chosen to qualitatively and quantitatively compare linear and nonlinear predictions: (Q1) Are the failure loads predicted either with the Pistoia criterion (PC) or the nonlinear simulations (NL) different from the experimental measurements? (Q2) Is the load sharing between the cortical and trabecular regions evolving differently over the length of the radius segments with the two methods? (Q3) To what extend do the local damage patterns from NL simulations differ from the overstrained volume predicted by the PC?
For points (2) and (3) an optimised Pistoia criterion (oPC) was used to enable a comparison of the potentials of linear and nonlinear lFE simulations.

Materials & Methods
CT scans and experimental load-displacement curves until failure were taken from a previous study (Hosseini et al., 2017): 21 radius segments of 20 mm height were obtained from 12 human donors. They were taken from the most distal parts of the radii (starting 5 mm below the distal subchondral plate, including the clinical 9 mm section used for Colles' fracture investigations). The segments were scanned at 16:4 lm resolution.
To reduce the computational costs, the scans were coarsened to 32:8 lm and segmented into bone and background voxels. All disconnected parts were removed. The two most distal and most proximal layers were removed to eliminate artificial damage localization caused by rough cut surfaces.
To obtain a lFE mesh, each voxel representing bone tissue was converted to a linear hexahedral element. The axial compression experiments were mimicked by platen-like boundary conditions: The most distal nodes were fixed. The most proximal nodes were axially displaced, and fixed in off-axis directions. Due to the large model sizes, all simulations were conducted with the highlyparallel lFE framework ParOSolNL (Stipsitz et al., 2019).

Nonlinear lFE simulations
A damage-based material model with previously identified material parameters was used (Stipsitz et al., 2019). It consists of (1) an isotropic, linear-elastic region (initial Young's modulus E 0 ¼ 10 GPa, Poisson's ratio m ¼ 0:3), (2) a nonlinear region where the material degrades based on a scalar damage quantity D, and (3) a failure region where the modulus is set to a small fraction of E 0 . The transition from the linear to the nonlinear regime is determined by an isotropic, quadric damage onset surface (adapted from (Schwiedrzik et al., 2013), shape parameter f 0 ¼ 0:3, critical damage D c ¼ 0:9). In region (2), isotropic hardening (hardening modulus E h ¼ 0:05 E 0 ) and a tension-compression asymmetry (via asymmetric tensile and compressive damage-onset strains e þ 0 ¼ 0:68%; e À 0 ¼ 0:89%) is included. The compressive displacements were applied in 0:05% strain increments. The simulations were stopped at the first drop in the apparent force-strain curve. The apparent strain was computed as the applied displacement in z-direction on the top layer divided by the initial height of the sample. The stiffness was computed from the first (linear) increment. The ultimate load was identified by the maximum in the apparent force-strain curve. The strains are displayed instead of the displacements to normalize to the actual length of the various wrists and better represent the effect of a strain-based NL simulation. For the estimation of the error in the ultimate load F nl max compared to the experimental ultimate load F exp max , the NL simulation results were stiffness adjusted using an individual initial modulus E adj 0 to match the experimental stiffness. This led to a scaled apparent force, while other quantities (strain, number of damaged elements, local damage distributions) were not affected. An element was classified as damaged if the damage D > 0.
The modelling strategy applied for the NL simulations has been previously validated for trabecular biopsies on the macroscopic level (Stipsitz et al., 2019). Good agreement was obtained for the radius segments ( Fig. 1), for further details see Appendix A and (Stipsitz, 2020).

Linear lFE with the Pistoia criterion
The bone tissue was modelled as a linear-elastic material with homogeneous isotropic material properties (E 0 ¼ 10 GPa; m ¼ 0:3).
The failure load was estimated using the PC (Pistoia et al., 2002): Based on the strain energy density U, a local effective strain is computed as e eff ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2U=E 0 p . The ultimate load is predicted by scaling the apparent load such that a certain percentage V c of the elements is strained beyond a critical strain e c .
In general, the appropriate values for V c and e c depend on the scanning technique and resolution of the lFE model, the segment size, and the initial modulus E 0 (Mueller et al., 2011;Arias-Moreno et al., 2019). Two different parameter sets were used in this study (Fig. 2): 1. Standard Pistoia: For question (1), the aim was the investigation of the predictive power of lFE simulations for the apparent ultimate load. Thus, standard parameters (e c ¼ 0:7%; V c ¼ 2%) were used for the evaluation of the PC (Pistoia et al., 2002;Varga et al., 2010). E 0 ¼ 10 GPa was used for all samples. Based on a linear simulation, an element was defined as ''damaged" (i.e. overstrained) if the effective strain e eff exceeded e c ¼ 0:7%. The evolution of damaged elements with the applied strain was computed in virtual increments (Fig. 2). The apparent strain obtained with the optimal Pistoia parameters was taken as maximum strain. Increments were obtained by dividing the maximum strain into 20 equally sized increments, and scaling the local effective strains accordingly.

Q1 -Statistics of the failure load
For statistical tests, all samples were treated as independent although most samples were taken from the left and right hand of the same donors. Previous studies using the same samples have treated them in a mixed statistics model (Hosseini et al., 2017) or have treated them independently (Arias-Moreno et al., 2019).
The statistical significance of the difference between values obtained with the PC and NL simulations were evaluated with the paired t-test with p < 0:05 (two-sided): To evaluate if the PC or the NL simulations can predict the failure load without calibrations, we tested if the failure loads from either the PC or NL simulations were significantly different (s.d.) from the measured failure loads. Additionally, linear regression analyses were performed between the failure load predicted by either the PC or NL simulation and the experimental values using python (package sklearn). T-tests were performed to check if the slope and intercept of the linear regressions differed from one and zero, respectively.

Q2 -Load sharing
To answer question (2), the progression of the load sharing between the cortical and the trabecular region over the length of the radii was evaluated as proposed by Johnson et al. (Johnson and Troy, 2018): The load was computed for each layer in z-direction of the radius segments as where hei r zz is the zz-component of the stress tensor of element hei, and hei A ¼ A 0 is the initial area of an element. Following (Johnson and Troy, 2018), the two most distal and most proximal layers were not included. The force contributions by the trabecular and cortical part were determined and normalized by the total force. The load sharing from the oPC and NL simulations were compared qualitatively at the ultimate load point over the length of the radius segments. Additionally, the average relative deviation between the load-sharing curves from oPC and NL simulations was evaluated as meanðjðF oPC z À F nl z Þ=F nl z jÞ, where the mean was evaluated over all samples, i.e. over all the data points of all the load contributionlength curves (see Fig. 4). A paired t-test was conducted to test whether this relative deviation is s.d. from zero, which would mean that there is on average a s.d. between the curves predicted by the oPC and NL simulations.

Q3 -Local damage patterns
For the qualitative local comparison of damage patterns between the oPC and the NL simulations (question (3)), only a binary value (damaged/not damaged) was compared element-wise. In the NL simulations, the amount of damage (number of damaged elements and damage values) strongly increased if the ultimate load was overestimated, and damage decreased if the load was underestimated. Thus, only samples for which the NL simulation resulted in an error (jðF nl max À F exp max Þ=F exp max j) in the ultimate load F max of less than 5% compared to the experimental value F exp max were included in the comparison. This was done to ensure a fair comparison between the oPC, which exactly reproduced the individual failure load by definition, and the NL simulations.

Q3 -Damage in cortical and trabecular region
The separation of images into cortical and trabecular regions was performed on segmented images with a coarsened resolution of 65:6 lm using the script manager medtool (v4.3, Dr. Pahr Ingenieurs e.U., Pfaffstätten, Austria) following the procedure in (Pahr and Zysset, 2009) (see Appendix B). For the linear simulations, the sum of the percentage of damaged elements in the cortical N cort D and trabecular region N trab D is equivalent to the critical volume V c .

Evolution of damaged and overstrained elements
Since the PC is defined in the number of overstrained elements N D , an investigation of the progression of N D with the applied strain was performed. Specifically, the evolution of N D in the trabecular region N trab D and in the cortical region N cort D determined via the oPC and from a NL simulation were compared. In both simulations, the number of damaged or overstrained elements did not represent a physical quantity but was a result of the numerical modelling.
The absolute values are not meaningful, thus, they were only compared qualitatively: ''Qualitative similarity" in the evolution of damaged elements from the oPC and the NL simulation was identified primarily via the transition point (e t ; N D;t ) of the cortical and trabecular damage evolution (marked by a cross in , simulation results with a uniform initial modulus (E0 ¼ 10 GPa) (blue, filled points), and individually scaled simulation results with E adj 0 (grey, open points) for a selected radius segment (for further samples see Fig. A2). The damage distribution at the last increment is shown, from white (no damage) to yellow (damaged) to dark red (fractured). The red square marks Fmax, the green diamond the 0.2%-offset yield strain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 7). At this point, the main damage location shifted from the trabecular to the cortical region. The evolution curves of the oPC and the NL simulation were categorized as similar if both curves showed a transition point (or both showed no transition point), and the location of the transition point e t was approximately at the same strain ratio (e t =e max ). The strain ratio was used instead of the absolute strain since the absolute strain in the inelastic region estimated from a linear simulation is not reliable.
To quantify the different local distribution of damage predicted either by the oPC or NL simulation for question (3), the proportion of damaged elements in the trabecular region (V trab ¼ N trab D =N tot D ) was evaluated at the failure point. We tested if the predicted V trab of the oPC and NL simulation was s.d. (paired t-test) and investigated systematic differences using a Bland-Altman plot.

Elasticity limit
To compare the role of the different inelastic onset criteria (PC: effective strain threshold vs NL: quadric damage onset surface), the inelastic limit e ie was computed analogously to (Stipsitz et al., 2019). In short, a piece-wise quadratic function is fitted to the individual N D -e curves using e ie as a fitting parameter. Thus, e ie is the strain at which a simulation predicts a decisive amount of damage, i.e. after which non-linearity plays a role locally. Here, e ie is a structural property since the radius segments cannot be considered representative volume elements.
More commonly the 0:2% strain-offset criterion is applied to estimate the yield load. The difference between the 0:2% yield load and the elasticity limit applied here is that the 0:2% yield load is determined purely on the apparent level while the elasticity limit is based on the amount of local damage. For the radii segments, the 0:2% yield load is close to the maximum load where already a decisive amount of local damage is present ( Fig. 1 and Fig. A2).

Q1 -Failure load
The linear regression of the ultimate load between simulations and experiments was comparably good using the standard PC and NL simulations with an unadjusted E 0 ¼ 10 GPa (Fig. 3). However, the standard PC (F stdPC max ¼ À1:9 AE 0:88 kN) largely underestimated the experimental ultimate load (F exp max ¼ À6:22 AE 3:49 kN). The difference was confirmed to be statistically significant by a paired ttest with p < 0:05.
The regression function of the standard PC was far from a oneto-one correspondence (F exp max ¼ À1:48 þ 2:2F sim max , t-tests confirmed that the slope was significantly different from one and the intercept significantly different from zero). To obtain a good agreement between the experimental results and the PC much higher critical volumes V c were required. This becomes obvious when comparing V c obtained with the oPC with the standard value from the literature, V c ¼ 2% (Table 1).
The prediction obtained in the NL simulations with standard parameters was closer to the values from experiments (slightly lower intercept, slope much closer to one: F exp max ¼ À1:29þ 1:32F sim max ). Nevertheless, the slope was significantly different from one and the intercept s.d. from zero. When NL results were stiffness adjusted, the ultimate loads were not significantly different from experimentally obtained values (Appendix A). In this case, slope and intercept were not s.d. from a one-to-one relation.

Q2 -Load sharing
The load sharing of the trabecular and cortical region obtained from a linear simulation (with the oPC) and a NL simulation with individually adjusted E adj 0 at the ultimate stress point was comparable for all samples (Fig. 4). On average the oPC led to a ð0:3 AE 0:2Þ % relative deviation compared to the NL simulations in the load sharing over length curves (This small difference was reported s.d. from zero by a paired t-test.) The largest differences were observed in sample R195. Generally, the trabecular region carried more load on the distal end than on the proximal end, and the cortical region vice versa. In most samples (17 out of 21), an intersection point existed. This indicates that on the distal end the majority of the load was sustained by the trabecular region, and on the proximal end by the cortical region. Missing intersections were predominantly found in weaker radii (Appendix C).

Q3 -Local damage patterns
The local damage patterns obtained with the oPC and the NL simulations showed small differences (Fig. 5). On average, ð21 AE 8Þ % of the elements where marked by both methods as damaged, while ð2 AE 1Þ % (only NL) and ð15 AE 4Þ % (only oPC) were identified by only one method. This indicates that the oPC   Fig. 4. Load sharing between the cortical and trabecular region over the length of the radius segments for six selected samples (0% relative length is most distal, 100% most proximal). Curves obtained from linear simulations with the optimal Pistoia criterion (LIN) and materially nonlinear simulations with individual E adj 0 (NL) were similar.
successfully identified damaged regions but while NL simulations allowed the yielding of elements in these regions, linear simulations led to more diffuse damage since no stress redistribution was possible. This is especially pronounced at the proximal end. In some samples, NL simulations led to a connection between two damage regions which was not present in the linear simulations (e.g. R191, to a limited extent also in R187). The bottom views indicate that in the NL simulations, the trabecular regions contained a larger number of damaged elements than in the linear simulations (e.g. R181, R191). However, this observation did not apply to all samples (for instance not for R188).

Q3 -Damage in cortical and trabecular region
The proportion of damage in the trabecular region at the ultimate load (V trab ¼ N trab D =N tot D ) predicted by the oPC and the NL simulations was significantly different. A Bland-Altman plot revealed that in most samples more damage was located in the cortical than in the trabecular region, i.e. V trab < 0:5 (Fig. 6). The average of the difference between the oPC and NL simulations V lin trab À V nl trab ¼ À0:07 was negative. This implies that on average the oPC led to slightly less damage in the trabecular region than the NL simulation.
N trab D and N cort D obtained with the oPC and the NL simulations showed qualitatively high agreements at the ultimate load point as well as in the progression curves over the applied strain (Fig. 7). Most evolutions (17 out of 21) were qualitatively similar. In three cases (R184, R185, R195) only either the NL simulations or the oPC displayed a transition point. In one additional sample (R194) the transition point obtained in the NL simulation occurred at a much smaller strain ratio than estimated with the oPC. The absolute values (strain range, number of damaged/overstrained elements, etc.) were not comparable.

Discussion
This study investigated whether apparent and local results obtained from linear simulations with two types of the PC, and NL lFE simulations are comparable for radius segments. The predictions of the maximum force using linear regressions were comparable irrespective of the parameters used. Locally, the PC led only to comparable results when an optimal individual calibration, which is practically not possible, was used.

Q1 -Failure load
The PC and NL simulations showed comparable linear regressions for the ultimate load with experimental data. However, computational demands were quite different: 424-3958 core hours for a NL simulation compared to 2-12 core hours for a linear simulation (85-380 mio DOF).
NL simulations without stiffness calibration achieved better absolute predictions of the ultimate loads than linear simulations. The PC with standard parameters largely underestimated the failure loads, as also reported in the literature (Zhou et al., 2016).  For practical applications, where one is interested in assessing the fracture risk, the PC is sufficient since the systematic error can be easily determined via a linear regression. Thus, Pistoia parameters are typically calibrated with the aim of achieving the best R 2 in a linear regression (Mueller et al., 2011;Arias-Moreno et al., 2019). Using this approach, one study also reported a regression function that was not significantly different from a one-to-one correspondence (Arias-Moreno et al., 2019). However, for the present data set, this calibration approach did not improve the prediction of absolute ultimate loads compared to the standard parameters. For this study the absolute failure load was important because it was the only means of estimating the quality of the local damage/ overstraining regions which were studied in question (Q3). Thus, an individual oPC was proposed which exactly reproduced the failure load. This approach led to substantially higher overstrained volumes than the standard V c ¼ 2%. However, the mesh resolution was higher than in the original study (Pistoia et al., 2002) and the literature suggests that, generally, a finer image resolution requires a larger V c (Mueller et al., 2011). In a previous study with comparable resolutions, a V c ¼ 19% was used (Zhou et al., 2016).
In both simulation types, the accuracy of the predicted ultimate loads depended on the material parameters. The experimental ultimate force was needed to evaluate the oPC. In contrast, the NL simulations did not necessarily require experimental data but it was possible to improve the failure load predicted by the NL simulations by using the experimental stiffness, which can be obtained in a non-destructive test. Nevertheless, some samples showed high deviations in the predicted ultimate load, but a calibration of the NL simulation beyond the initial modulus was not easily possible.

Q2 -Load sharing
The load sharing between the cortical and the trabecular region depended much more on the individual microstructure than on the tissue model (Fig. 4). Thus, for simple structures as the radius a materially linear simulation seems to be adequate to estimate the normalized load sharing or evolution of the apparent forces.
The general trend to more load in the trabecular region at the distal end, and more load in the cortical region at the proximal end was confirmed (Johnson and Troy, 2018). However, some samples exhibited no transition point. It seemed that weaker samples (lower load at failure) exhibit earlier or no transition compared to stronger samples (see Appendix C). This could be evidence of advancing osteoporosis where it is known that trabecular struts resorb and the load transfer is taken up by the cortical compartment (Rüegsegger et al., 1991).

Q3 -Local damage patterns
Damage patterns obtained from the NL simulations and the oPC were observed at similar locations. This can be explained by the definition of a damaged element as D > 0 or e eff > e c and quite similar values for the damage-onset strains e AE 0 and e c . However, higher damaged areas could not be investigated via the PC. Furthermore, local results with the oPC cannot be obtained independent of experiments since they rely on an accurate individual calibration. In contrast, NL simulations did not require any material parameter calibration to obtain the damage patterns. The local damage distributions are modelling results which were not validated experimentally. In experiments, damage can be visualized as whitening of bone tissue (Thurner et al., 2007) or by sequential staining (Lee et al., 2000). To obtain threedimensional images lCT scans can be performed repeatedly during loading (Müller and van Lenthe, 2006). However, this is very time consuming and the microdamage is not accessible at a subtrabecular resolution (Mueller et al., 2011;Costa et al., 2017;Martelli and Perilli, 2018).

Q3 -Damage in cortical and trabecular region
The extent to which damage was shifted from the trabecular to the cortical region between a materially nonlinear and linear simulation was small (Fig. 6). This was surprising since in a NL simulation, after damage onset the failure of individual trabeculae and the subsequent stress-relocations could have easily led to larger damaged regions in the trabecular region compared to a linear simulation.
The distribution of damaged elements over the trabecular and cortical region (Fig. 7) showed larger differences than the load sharing (Fig. 4). This could be expected since the load distribution is a physical quantity while the percentage of damaged elements is a modelling parameter. Results indicate that the percentage of overstrained elements may not be unique. However, the PC is based on this assumption since it was formulated in V c (Pistoia et al., 2002). This finding could explain why it is necessary to calibrate the PC for most new data sets.
The evolution of the number of damaged elements N D was qualitatively similar between the oPC and NL simulations. The absolute strains from oPC were mostly smaller than from the NL simulations since they were based on a single linear simulation. The absolute values of N D were not comparable since different criteria were used to decide if an element was damaged: a strain-threshold (PC) or a quadric damage-onset surface (NL). The difference between the two criteria is obvious when comparing the inelastic limit e ie , which marks the point after which the simulations predicted that nonlinearity plays a role in the structure (Fig. 7). Most curves suggest that in healthy radii, trabecular bone absorbs the damage first and protects the cortex until the intersection where the cortex receives more damage and causes failure of the structure.

Limitations
As far as possible similar modelling strategies were chosen in the linear and NL simulations so that a comparison of results is legitimate. Most importantly, simplified boundary conditions were used. Spurious damage was found on the top and bottom layer and may lead to differences in the load distributions (Johnson and Troy, 2018). However, these boundary conditions are widely used in the literature for the distal radius (Vilayphiou et al., 2011;Hosseini et al., 2017) and seem to be meaningful (Varga et al., 2010).
Additionally, only the impact of materially nonlinear modelling using a damage formulation was investigated. Geometric nonlinearity and other aspects of the constitutive behavior of bone (viscosity, poroelasticity, plasticity) were neglected. Thus, the reported results have to be seen as a first step. Due to high computational demands and availability of suitable specialized solvers it was not possible to investigate the impact of a wider range of material models.
It is not clear to what extent the results can be transferred to other anatomic sites. The PC was originally used for Colles'-type fractures of the distal radius (Pistoia et al., 2002) but was also successfully applied to the tibia (Vilayphiou et al., 2010). Generally, stress redistributions could play a much larger role in whole bones or under non-trivial loading.

Conclusion
Good agreement between linear (oPC) and NL (damage-based) simulations was found on all levels for distal radius segments. For practical applications, where one is interested in the failure load (as in Q1) or apparent properties like the load distribution (as for Q2) the PC with adequate normalization and calibration of the results was sufficient. If the failure load is known from experiments, the oPC can even estimate local properties (as for Q3) although the predicted overstrained regions were more diffuse compared to a NL simulation. These local damage distributions could be interesting for basic research, for instance, to investigate the impact of structural changes on bone strength. The benefit of NL simulations is that the damage localization is represented and obtainable even if no experimental measurements until failure are available, although, at much higher computational costs.

Declaration of Competing Interest
DP is CEO of Dr. Pahr Ingenieurs e.U. which develops and distributes the software medtool. MS and PZ have no conflicts of interest to declare.

Acknowledgments
The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC). PZ would like to gratefully acknowledge SNF Grant No. 165510.

Appendix A. Apparent level agreement of nonlinear and experimental data
While the PC is widely standardized for radii segments, nonlinear lFE simulations can be performed with different settings depending on how the nonlinear material behaviour is modelled.
To validate the nonlinear lFE approach applied in this study, additional paired t-tests were conducted to test if the stiffness and stiffness adjusted failure load differ from the experimental measurements. This information should facilitate the comparison of different nonlinear lFE approaches.
The apparent stiffness and adjusted ultimate load obtained from the NL simulations showed good correlation with experiments (Fig. A1). The concordance correlation coefficient was 0:87 for the stiffness and 0:92 for the ultimate load. Both stiffness and ultimate load were generally underestimated using E 0 ¼ 10 GPa.
Using E adj 0 , the ultimate load was slightly overestimated. All apparent measures (stiffness, yield load and ultimate load with uniform or individual E adj 0 ) were not significantly different from the experimental results (Table A1). Absolute errors were not notably higher in weaker or stronger radii.
Qualitatively, most force-strain curves matched quite well to the experimentally obtained curves (Fig. A2). However, some samples showed high deviations between experiments and simulations.
The initial tissue modulus had a strong impact on the one-toone predictive power for the ultimate load while not affecting the coefficient of determination. To illustrate this effect, apparent forces were scaled to match the apparent stiffness from experiments. These are also shown in Fig. A2. In the individual simulations, E 0 of 10:56 AE 2:51 GPa (confidence interval: ð9:38; 11:73Þ) were required to reproduce the experimental stiffness.
In some samples an exceptionally good agreement between simulated and experimentally obtained force-strain curves was achieved (e.g. R191 in Fig. A2) while in other samples, no resemblance was found (e.g. R180). The shape of the curve suggests a poor contact between the sample and the loading platten. Experimental artefacts may also be at the origin of the mismatch between the curves. Other possible reasons could be natural variation in the apparent modulus, variations in the sample preparation and testing, or structural details that were not resolved at  the scan resolution which led to damage localization. Also, it was assumed that there was no initial damage in the structures, for instance introduced by the preparation process.
Generally, simulations overestimated the stiffness of low stiffness samples, while under-estimating the stiffness of high stiffness samples. Thus, lFE was not able to reproduce the full variability  (Table A1). The proposed material model did not include all aspects relevant for reproducing the experimental results: First, variations in the stiffness and ultimate load obtained in the simulations were generally lower than variations in the experiments (Table A1). Second, the simulated stress-strain curves showed an abrupt drop after the ultimate load rather than a slow decrease which was found in the experiments (Fig. A2). This abrupt drop of the simulated apparent stress is caused by the discontinuous failure behaviour in the material model. However, in the literature a similar simple material model (an asymmetric bilinear model in principal strain) was able to reproduce the post-ultimate behaviour (MacNeil and Boyd, 2008). Third, the failure strains obtained from simulations was much more constant than from experiments (À1:07 AE 0:09% vs À1:25 AE 0:3%).
Interestingly, the 0:2% yield point was close to the ultimate load point in all samples. Thus, it is not surprising that the yield point has been found to be a good predictor of strength for radius segments (MacNeil and Boyd, 2008). Generally, the yield point is used in combination with an elasto-plastic material behaviour. Since the location of the 0:2% yield point depends on the applied material model, it cannot be concluded that its location is similar using a damage-based or an elasto-plastic material model.

Appendix B. Separation of cortical and trabecular volume
The radius segments were extended on the top and bottom, and large pores were filled using a morphological closing operation. Outer and inner masks were obtained using a ray based filling operation. The cortical mask at the coarse resolution was obtained by subtracting the inner mask from the outer mask. Finally, the fine mask at 32:8 lm resolution was obtained by refining the coarse mask and mapping it onto the original segmented image.
The number of damaged elements in the cortical and trabecular region was obtained by applying the fine masks to the resultant damage distributions. The percentage of damaged elements N D is normalized to the total number of tissue elements in the image. . Part 1: Evolution of the load sharing between the cortical and trabecular region over the length of the radius segments (0% of length is most distal, 100% most proximal). Evolutions obtained from linear simulations with the oPC and materially NL simulations were similar.