Generic finite element models of human ribs, developed and validated for stiffness and strain prediction-To be used in rib fracture risk evaluation for the human population in vehicle crashes

To enable analysis of the risk of occupants sustaining rib fractures in a crash, generic finite element models of human ribs, one through twelve, were developed. The generic ribs representing an average sized male, were created based on data from several sources and publications. The generic ribs were validated for stiffness and strain predictions in anterior-posterior bending. Essentially, both predicted rib stiffness and rib strain, measured at six locations, were within one standard deviation of the average result in the physical tests. These generic finite elements ribs are suitable for strain-based rib fracture risk predictions, when loaded in anterior-posterior


Introduction
For the development of occupant restraint systems, anthropometric test devices (ATDs) are generally used. To complement the ATDs, usage of detailed finite element (FE) human body models (HBMs) is increasing rapidly. The HBMs can be used for detailed injury analysis, such as evaluation of rib fracture risk for an occupant in a crash. However, for such analyses detailed and validated models of each rib, a model that combines these rib models into a ribcage, and a model representing the entire thorax, are required. Human ribs are complex. Each rib consists of two types of bone, cortical bone and trabecular bone. The cortical bone is the dense outer surface, while the trabecular bone is the porous internal structure of each rib. In addition, the overall rib dimensions as well as the cortical bone thickness vary along the length and around the circumference of each rib.
Today, there are several HBMs representing average sized male occupants, of which the two major HBMs are referred to as the Total Human Model for Safety (THUMS) AM50 (JSOL, Tokyo, Japan) and the Global Human Body Model Consortium (GHBMC) M50-O (Elemance, Clemmons, NC, USA). These HBMs were developed based on geometrical data from medical imaging and measurements of specific individuals. For THUMS AM50 and GHBMC M50-O, individuals representing the average sized male were chosen based on external measurements, e.g., stature, seated height, shoulder breadth, and body mass [Shigeta et al. (2009) ;Iwamoto et al. (2015); Gayzik et al. (2010); Gayzik et al. (2012); Kindig (2009)]. It was not reported if the overall rib dimensions or rib cross-sectional properties were taken into consideration during the selection process. Therefore, it remains unknown if the dimension of each rib has been based on averages derived from population studies, which presumably would be important for estimating the average rib fracture risk in the population. However, varying cortical bone thickness was implemented in the GHBMC M50 O but not in any of the THUMS AM50 model versions.
In recent years, FE morphing has been used to modify the geometry of THUMS AM50 [Hwang et al. (2016) and Zhang et al. (2017b)] and the GHBMC M50-O [Zhang et al. (2017a) and Hu et al. (2017)] models to represent a larger portion of the population. This enables optimisation of safety systems to humans deviating from the averaged sized male. In these studies, stature, body mass index (BMI), age and sex have been used as independent variables to control the morphed geometry. It was shown in Zhang et al. (2017a) that all these parameters also significantly influences the chest force deflection response in a longitudinal pendulum impact. However, in a study by Wang et al. (2016) it was shown that these four parameters only accounted for some of the ribcage geometry variance (R 2 ¼ 0.51). Thus, to further improve the morphing methods for the ribcage, it is important to investigate which rib parameters that influence the ribs the most.
Two sources of data are available for creating ribs based on averaged properties; in-vivo data from Computer Tomography (CT) and in-vitro rib measurements. To enable creation of a whole ribcage at a level of detail necessary for rib fracture analysis, the data must be defined for all ribs and at several locations along each rib. To develop a statistical ribcage, clinical CT data from 63 individuals was used in [Gayzik et al. (2008)]. In this model, the whole ribcage was described using 106 landmarks. In a similar study using clinical CT data from 89 individuals, Shi et al. (2014) showed that more landmarks were required to capture details such as rib curvature and cross-sectional dimensions of the rib. Shi et al. (2014) suggested defining the rib by means of 28-44 landmarks, grouped in quadruplets, where each quadruple defines the height and width of the rib cross section. In another study, Weaver et al. (2014) collected data from 339 individuals and used between 2700 and 10 400 landmarks to describe each rib. Local rib features were captured with this resolution. Studies based on clinical CT, similar to the three mentioned above, are useful for evaluating the overall shape of the ribs but less useful for the evaluation of the rib cross-sectional dimensions. Perz et al. (2013), as well as Holcombe et al. (2019) showed that the cross section will generally be overestimated when based on clinical CT. Instead, micro-CT or in-vitro studies can be used to evaluate rib cross-sectional dimensions. A detailed data set was developed by Choi and Kwak (2011) based on seven male subjects. Both cross-sectional dimensions and cortical bone thickness, at nine locations along the rib, was measured for all twelve ribs using micro-CT by the authors. This is the only data set in which rib cross-sectional dimensions and rib cortical thickness have been measured in the same individuals. This may be of importance as it has not been shown if rib cross-sectional dimension and cortical bone thickness correlate. This data set is, to the best knowledge of the authors, the most detailed data set for rib cross-sectional dimensions and cortical thicknesses suitable for creation of ribs with average shape and cortical thickness.
In THUMS AM50 and GHBMC M50-O, rib fractures are evaluated based on strain-based element erosion. Forman et al. (2012) proposed an alternative to element erosion to predict rib fracture risk, using a probabilistic framework based on predicted strain. Using strain as a predictor for rib fracture is in line with the results of Trosseille et al. (2008). In their study a "strong correlation between the axial strain of the rib external cortical bone and the fracture location" for longitudinal, oblique and lateral chest loading was observed. However, published HBM validations have mainly focused on stiffness rather than strain responses. For example, the thorax of THUMS AM50 Version 4 was validated with regard to stiffness on whole body level [Shigeta et al. (2009)]. Thorax stiffness validations were published for the THUMS AM50 Version 3 and the GHBMC M50-O at different structural levels, from single ribs up to the whole thorax [e.g., Mendoza-Vazquez et al. (2013); Li et al. (2010c); Poulard et al. (2015); Vavalle et al. (2013); Vavalle et al. (2015)]. Although rib strain validations were not published for any of these models, rib strains have been evaluated indirectly by comparing the number of predicted fractured ribs to corresponding experimental PMHS data for a few load scenarios [Shigeta et al. (2009);Schoell et al. (2015a); Schoell et al. (2015b)]. Shigeta et al. (2009) reported less predicted rib fractures with THUMS AM50 Version 4 compared to the average number fractured PMHS ribs in Kroell type frontal pendulum impacts. Also [Schoell et al. (2015a) and Schoell et al. (2015b)], reported fewer predicted rib fractures with the GHBMC M50-O compared to PMHSs in frontal and lateral impacts. Conversely, reconstructions of real-life crashes using an updated THUMS AM50 Version 3 showed that this model over-predicted the rib fracture risk [Iraeus and Lindquist, 2016;Mendoza-V� azquez et al., 2014]. Hence, to enable biofidelic rib fracture risk prediction with HBMs, validating the model for both stiffness and strain prediction seems necessary.
Validation data, including strain measurements, for single ribs loaded longitudinally have been presented in [Charpail et al. (2005), Kindig (2009) and Agnew et al. (2018)] and laterally loaded in de Dios et al. (2011). Of the three studies comprising longitudinally loaded ribs, the Kindig (2009) study was described in such level of detail that it was the easiest to replicate in a simulation.  validation data was discarded in this study due to uncertainties in the boundary conditions.
The aim of this study is to develop generic models of ribs one to 12, based on averaged geometrical and material data, and to validate the models for both stiffness and strain to enable biofidelic strain-based rib fracture risk predictions. The aim has been extended to include definition of the rib dimensions and material properties that influence the variability in the anterior-posterior bending response. This is an initial step in the development of a complete ribcage model that can be morphed to represent variability in different occupants beyond the average male.

Method
The development of generic FE meshes of ribs one through 12 is described in Section 2.1 below. Rib cross-sectional dimensions and cortical thicknesses are based on Choi and Kwak (2011) while the overall rib dimensions and shapes are based on the statistical model developed by Shi et al. (2014). Section 2.2 describes the material properties assigned to the rib models based on experimental data from Kemper et al. (2005) and Kemper et al. (2007). Section 2.3 describes simulations performed with the models to validate model stiffness and predicted strain to experiments performed by Kindig (2009). Finally, Section 2.4 presents the simulations performed to assess influence of geometrical and material parameters on the strain response in rib end-to-end bending tests.

Rib geometry
The twelve generic ribs were created in four steps, as is schematically shown in Fig. 1. Rib dimensions from the data set by Choi and Kwak (2011) and the rib curvature and twist from the statistical model by Shi et al. (2014) were used.
For each rib the following steps were performed using the preprocessor ANSA 16.1.1, Beta CAE Systems (2016): Step 1 A line, the length corresponding to the length of the rib currently being modelled, was created based on Choi and Kwak (2011). Along this line, nine ellipses representing the rib outer dimensions, were created with half axes and positions according to Choi and Kwak (2011). The ellipse dimensions for the two end positions were created by linear extrapolation of the cross-sectional dimensions using data from the two neighbouring cross sections. Thus, in total 11 elliptical cross sections were defined for each rib.
Step 2 The cortical thickness was defined at 16 locations along the perimeter for each of the nine cross sections according to Choi and Kwak (2011). The two cross sections at each end were given the same thicknesses as the adjacent cross sections. The 16 points described by the inner cortical surface of each cross section were connected using a spline.
Step 3 The 11 cross sections were used to create an inner and an outer surface for the cortical bone. To enable modelling the cortical bone as shell elements in the mid layer, the cortical bone was initially meshed using a single layer of solid elements. Nodes were placed at the exact locations of all the points defined in Step 2. The average element side length of the resulting mesh was 2.5 mm. The trabecular bone was meshed using solid elements. Finally, the solid elements representing the cortical bone were converted to thin shell elements, positioned halfway between the top and bottom surfaces of the cortical bone. In this process the thickness for each node (LS-DYNA shell element formulation 16) was defined based on the solid element thickness.
Step 4 The straight rib from Step 3 was finally morphed to the actual shape of the rib. The statistical ribcage model defined in Shi et al. (2014)  Thus, these landmarks were used to define both the rib curvature and twist. The detailed morph procedure used for each rib is illustrated in Fig. 2. Using this morph procedure the ribs were stretched uniformly along the length axis to fit the rib length indirectly given by Shi et al. (2014), while keeping the cross sections defined in Step 1 unaffected.

Rib material data
Four published studies presenting material properties of human rib tissue were evaluated [Kemper et al. (2005) Albert et al. (2017) presented results from 37 tissue samples from 29 donors. However, as no information on material yielding was presented it was not possible to create a complete stress-strain relationship up to material failure, which is needed for analysing the response of ribs to crash type loadings. Only the two studies published by Kemper et al. [Kemper et al. (2005); Kemper et al. (2007)] included the complete stress-strain histories. Therefore, this data was reanalysed to calculate parameters for a non-linear material model. The material data was parametrised using a minimum number of parameters, i.e., the non-linear stress-strain curves were assumed to be  bi-linear.
All 163 stress-strain curves presented in Kemper et al. (2005) and Kemper et al. (2007) were digitalised, as shown in Fig. 3. For each curve, the initial slope (corresponding to Young's modulus E) was estimated with a linear regression model using all data points below 3000 microstrains. The tangent modulus (TM) was estimated as the slope in a second linear regression model using all data points above 10 000 microstrains. The two curve fitting regions were selected after analysing the shapes of a large number of stress-strain curves. The yield point was estimated as the intersection between the two linear regression models, i.e., where the two curves in Fig. 3 intersect. Then, the average of the material parameters was calculated in two steps. Firstly, the average for all tissue samples from the same subject was calculated. Secondly, the average based on each subject's average was calculated to obtain an average representative of the whole population (12 subjects).
These resulting material parameters were used to develop an isotropic plasticity material definition [LS-DYNA MAT24, Hallquist (2006)]. To complete the material model with all parameters needed for the yield surface, Hansen et al. (2008) data for strain rate scaling of the yield surface was used.
The trabecular bone was also modelled using an isotropic plasticity material definition (MAT24) with Young's modulus, Poison's ratio and yield Stress from Kimpara et al. (2005). Post yield strain hardening parameters were taken from Zhao and Narwani (2005).

Rib model validation
The isolated rib end-to-end bending tests performed by Kindig (2009) (test series RIBBEND 2 and 3) were simulated using LS-DYNA SMP Version R8.0.0, LSTC (2015) to validate the rib stiffness as well as the rib strain distribution for ribs two to 10. The simulation test setup can be seen in Fig. 4. In the simulation, each rib end was fully constrained in a rigid pot of similar size as was used in the physical tests. The end of the rib was positioned at the centre of the potting of the simulation model, as described for the physical tests. The anterior rib end was offset laterally to the posterior end. The magnitude of this offset was taken as the average offset for all physical tests at each rib level in the physical tests. The pots were allowed to rotate freely about the y-axis and were fully constrained in all other degrees of freedom. The anterior potting was moved posteriorly at a velocity of 1 m/s, in the simulations using prescribed motion of the potting. The reaction force was measured at the posterior end, using the boundary forces constraining the posterior potting from movement. The strain was measured in six cortical shell elements, three located at the pleural side and three at the cutaneous side, at equivalent locations to the experimental strain gauges, see Fig. 4. The physical strain gauges were uniaxial, aligned along the long axis of the rib. Correspondingly, in the FE model, the rib strain was measured in the element local coordinate system, in the element direction approximately aligned along the long axis of the rib. Details about exact strain measuring points can be seen in Appendix A.
The resulting model rib stiffness was calculated based on the output reaction force and prescribed displacements. It was compared to the experimental corridors for rib stiffness presented by Poulard et al. (2015), in which the data from Kindig (2009) was processed to create corridors according to the method proposed by Lessley et al. (2004). The model's output strains were compared to the experimental average of the maximum strain measured in all tests for each of the six strain gauges [Kindig (2009)]. However, since the ribs failed at different displacements these strain measurements were extracted from ribs with different deformations. The strain predictions derived from the simulations, were extracted at a rib deformation equal to the average rib deformation at rib failure in the experiments.

Rib parameter study
A parameter study was carried out for the simulation of the rib endto-end bending test of rib seven to explore the sources of variance in rib strain in the experimental data by Kindig (2009). Geometrical parameters, defined in the Kindig (2009) study include: rib arch height and chord length, cross-sectional width and height and the rib's initial rotation around the posterior end, defined according to Fig. 5 with parameter variations according to Table 1. In addition, the peak normalised deformation was varied according to the data for rib seven in Kindig (2009).
Finally, rib material parameters: Young's modulus, tangent modulus and yield stress, were varied according to the results of this study and presented in the Rib material data results in Section 3.2, below. All parameters were assumed to follow a Gaussian distribution.
The parameter study was carried out as a stochastic simulation using Latin hypercube sampling, n ¼ 50, in LS-OPT Version 5.2 Stander et al. (2015), including a pre-processing step, to apply the parameter changes. As many of the parameters require complex geometrical changes, the modification of the rib dimensions was carried out using geometrical morphing in the pre-processor step, in accordance with Fig. 6. Initially, the rib centreline was modified to adjust the rib arch, chord length and rotation, by morphing a curve to the desired dimension and position. Then the centreline of a cylindrical morph box, containing the rib's finite elements, was fitted to the morphed curve. The use of a cylindrical morph box ensures that the rib cross section does not change when the overall dimensions change. Next, the anterior support was translated and rotated to match the updated rib dimensions. Finally, the rib cross-sectional width and height were morphed using rectangular morph boxes contained within the circular morph box and aligned with the principal axes of the rib cross section. When changing the rib cross section, the same cross-sectional dimensions offset was applied along the whole rib length.
For each strain gauge, a multiple linear regression model, excluding interaction terms, was fitted to the results of the 50 stochastic simulations. Independent variables include the four geometrical and two boundary condition parameters in Table 1, and the three material parameters defined in the rib material data section. All variables were kept in the regression models, regardless of being statistically significant or not. However, as the variation of each variable differs, each of the regression model coefficients was multiplied by standard deviation for the corresponding variable, to quantify to which extent each variable contributes to the total variability. A more thorough explanation is provided in Appendix C. To summaries, the result of this analysis revealed the level of variation in strain gauge x, explained by the variation in parameter y.
In addition, the total variances of the predicted rib strains (SG1-SG6) derived from this stochastic analysis were compared to the total variances in the results from the physical experiments. Statistical analysis was conducted using R Version 3.2.3, R Core Team (2017).

Rib geometry
The shape and cortical bone thickness of the generic ribs, one through 12 can be seen in Fig. 7. The thickness of the shells representing cortical bone ranges from 0.2 to 1.65 mm.

Rib material data
The material data re-analysis resulted in an average Young's modulus of 14.7 � 2.0 GPa, an average tangent modulus of 1.94 � 0.5 GPa, and an average yield stress of 100.7 � 12.9 MPa. Additionally, correlation between these parameters were analysed, and the only significant correlation occurred between Young's modulus and the yield stress, with Pearson correlation coefficient r ¼ 0.77 (p ¼ 0.005). Please refer to Appendix B for further details.

Rib model validation
The force-deformation curves for the generic ribs were compared to the corridors created from Kindig (2009), see Fig. 8. The predicted force-deformation curves remain within the corridors, with the exception of ribs six and seven at below 30 mm displacement and rib nine above 40 mm displacement. Fig. 9 compares the simulation and the averaged experimental rib strain measured just prior to fracture. Forty-nine out of 54 (91%) of the simulation strain results were within one S.D. of the average of the physical test results. Of the five strain measurements outside the �1 S.D. range, most were right on the border. The strain prediction deviating the most from the physical test was the predicted strain in gauge 2 on Rib 5. Fig. 10 displays the contribution to the strain gauge variance in SG1 to SG3, i.e., the level of variation in strain gauge x explained by the variation in parameter y. The results have been normalised to the parameter with the greatest contribution, i.e. the rib cross-sectional height. Thus, Fig. 10 shows the parameters relative importance in explaining the total strain variance, namely, the level of importance of the parameter for predicting the correct strain result. After the rib crosssectional height, the deformation at failure was the second most important parameter. Additionally, the rib arch height, rib rotation, rib cross-sectional width and rib chord length parameters were also found to contribute to strain variance, while the contribution of the material parameters was found to be less significant. Please refer to Appendix C for further details.

Rib parameter study
In Fig. 11, the total variance (�1 S.D.) of rib strains in the simulation (red) parameter study was compared to the total variance in the physical test (grey and black) for the six strain gauges of rib seven. As the length of the red line was shorter than the length of the grey and black lines, it can be concluded that the full variance in the results in the physical tests   had not been captured by the parameters in the study. However, for the posterior strain gauges (SG3 and 6), most of the variance was captured.

Discussion
This study indicates that with correct geometrical modelling of the rib cross sections, rib strains can be accurately predicted in anteriorposterior rib bending. The rib parameter study shows that the rib cross-sectional height, defined according to Fig. 5, was the most influential parameter for rib strain. The second most influential parameter was the deformation at failure, the displacement of the rib end. Other geometrical rib properties, i.e. rib arch height, rib chord length and rib rotation (governed by the amount of costal cartilage), also influenced the rib strain, but to a lesser extent. The material parameters did only affect the results marginally. On the one hand, comparing these results to an analytical solution for peak strain in a hinged beam arc with constant cross section and linear elastic material properties, see Equation (1) (variables are explained in Fig. 5), it can be seen that, similar to the parameter study, the strain is only a function of geometrical parameters.

Rib peak strain ¼ Constant⋅
ΔL⋅h H⋅L (Eq. 1) On the other hand, the analytical solution for the stiffness of the same hinged beam arc, subjected to a prescribed displacement, see Equation (2), shows that both geometrical (H, I, L) and material (E) parameters will influence rib stiffness.

Rib stiffness
The test setup used for the validation in the current study, subjected the ribs to a prescribed displacement. However, in a real crash the rib peak deformation will be determined by the actual crash pulse, the interior safety system and the total stiffness of the thorax. According to Kent (2008), in a 48 km/h frontal impact, the ribs together with the intercostal muscles only account for one third to one half of the total stiffness of the thorax. Essentially, to estimate peak rib deformation adequately, established as the second most important parameter for accurate strain prediction, the complete thorax stiffness must be accurately modelled. Additionally, thorax stiffness is mostly influenced by  anatomical structures other than the ribs.
The average material properties derived in the current study, in particular the Young's modulus of 14.7 GPa and the yield stress of 100.7 MPa, were higher than utilised in other current HMBs, e.g., the GHBMC M50-O model with Young's modulus of 11.5 GPa and yield stress of 88 MPa (Li et al. (2010a)), or the THUMS AM50 model with Young's modulus of 9.9 GPa and yield stress of 67 MPa (Shigeta et al. (2009)). As seen in Equations (1) and (2), the material parameters were mainly  Fig. 9. Rib strain at peak force in elements located approximately at SG1 to SG6, compared to experimental results from Kindig (2009). The vertical lines represent the �1 S.D. peak strain measured using strain gauge one through six for ribs two through 12. The red dots represent the strain measured in elements in the simulation model, at similar positions, at the same deformation as when the fracture occurred in the physical tests.
influenced the rib stiffness, and thus the force-deformation response. Most of the generic rib force-deformation curves in Fig. 8 were within the physical test corridors. Two of the curves (ribs 6 and 7) were slightly below the corridor for deformations up to 30 mm while one curve (Rib 9) was above the corridor for significant deformations. Further, comparing the force-deformation curve shape with the shape of the corridors might indicate that the rib models lack certain non-linearity, and thus the yield stress and/or tangent modulus are incorrect. However, revisiting the raw curves used to develop the corridors, which are published in Li et al. (2010c), shows that very few of the ribs displayed the non-linear behaviour suggested by the corridor. In fact, many of the ribs showed a non-linearity similar to the response of the generic ribs in the current study. Altogether, the results of this study indicate that the average material parameters, derived in and used for this study, are reasonable.
In this study, six geometrical and three material parameters were included. Together with the parameter distributions defined in this study, the total variance for the simulated strain was underestimated compared to the physical tests, in particular for the laterally positioned gauges, shown in Fig. 10. This implies that there can be other, excluded, parameters important for accurate strain prediction. An example of this can be seen in the variation of the boundary conditions, such as the distance between the centre of rotation to the rib ends, or the variation in rib potting material stiffness. Variance in the position and orientation of the physical strain gauges may also introduce uncertainty in the strain results. Furthermore, variance in the distribution of trabecular bone will affect the cortical bone strains, (Iraeus et al., 2019). However, one parameter that most likely is very important, also omitted from this study, is the variability of thickness of the cortical bone. In reality, the thickness and distribution of the cortical bone will vary between individuals, as illustrated for one rib in Choi and Lee (2009). Indirect evidence of the importance of cortical thickness distribution can be seen in Fig. 12. This figure shows the thickness distribution around SG2 in Rib 5, the strain gauge with the greatest deviation from the physical measurements (Fig. 9). The thicknesses highlighted are the thicknesses based on the micro-CT measurements from the Choi and Kwak (2011) dataset, from which all thicknesses in the generic ribs have been interpolated. For Rib 5, at this specific location, a local reduction in thickness is present in the Choi and Kwak (2011) dataset. As seen in the right sub-figure, this local reduction in thickness has produced a strain hotspot at this particular location. Future studies should quantify and include this inter-individual variance.
Cortical bone thickness varied between 0.2 and 1.65 mm in this study. The thickness is comparable to the thickness in the Kalra et al. (2015), in which the average thickness for male subjects was 1.09 mm. However, the variability of cortical bone thickness along the rib, had not been evaluated.
As discussed above, accurate geometric modelling of the rib cross section, in particular the rib cross-sectional height, is a prerequisite for strain based rib fracture modelling. The cross-sectional dimensions of the generic ribs created in this study are based on micro-CT and should in theory be more accurate than ribs created using clinical CT. To evaluate this hypothesis, the cross-sectional height of the generic ribs, as well as the ribs from THUMS AM50 and GHBMC M50-O were compared to other published data [weighted average of Mohr et al. (2007) and Kindig (2009)], see Table 2. The cross-sectional height of the generic ribs were within �5% of the weighted average of the compared data. On average, the THUMS AM50 Version 4 model overestimates the rib cross-sectional height by 22%, and the GHBMC M50-O by 18%. These figures agree well with the findings of Perz et al. (2013) and Holcombe et al. (2019), who found that cross-sectional rib dimensions based on simple thresholding of clinical CT will generally be overestimated. Yet, the rib cross-sectional height in the THUMS AM50 Version 3 and 5   11. Range (�1 S.D.) of peak rib strains in Rib 7 (red lines), SG1 to SG6, due to variation in four geometrical and two boundary condition parameters defined in Table 1 and the three material parameters defined in the rib material data section, compared to the variation in physical tests (black and grey lines) of Rib 7 presented in Kindig (2009). models has been found to underpredict by 26% ( Table 2). The modelling of the THUMS AM50 Versions 3 and 5 has predominantly been based on a commercial dataset provided by Viewpoint Datalabs (Iwamoto (2002). However, details explaining the reason for the underprediction has not been provided. Over or underestimation of rib cross-sectional height to such an extent may be one of the reasons the THUMS AM50 and GHBMC M50-O models over or underpredicts the number of fractured ribs.
Using micro-CT or other in-vitro rib cross-sectional data sources is thus beneficial to obtain unbiased cross-sectional dimensions. However, as only gross cross-sectional dimensions (height and width) were presented in the Choi and Kwak (2011) dataset used for the current study, it is impossible to recreate the actual cross-sectional shape. The assumption adopted throughout this study is that the shape of the rib cross section can be approximated using an ellipse, with half axes corresponding to the rib cross-sectional height and width. The elliptical assumption has already been analysed and justified in previous studies, e.g., Roberts and Chen (1970) and Choi and Kwak (2011).
For anterior-posterior bending the most important cross-sectional property is the bending moment of inertia around the weakest axis (2-2 axis in Fig. 5). In Appendix D the true bending moment of inertia in the Choi and Kwak (2011) dataset has been compared to the generic ribs, to other current HBMs, and to other published data. Generally, the bending moment of inertia of the generic ribs closely follows the true values based on the Choi and Kwak dataset for the anterior 70% of the rib length, while for the posterior 30% of the rib length the approximation either underestimates or overestimates the true value. This finding means that rib strain prediction may be less accurate for the 30% most posterior aspect of the ribs. In fact, four out of the five rib strains, not predicted within one S.D., see Fig. 9, was located closest to the posterior end. In a number of recent studies based on real life data, most rib fractures were found in the anterolateral or lateral aspects for frontal impacts, and the anterolateral, lateral or posterolateral aspect for lateral impacts [Lee et al. (2015) and Ejima et al. (2017)]. Hence, using the generic ribs for analysis of frontal impacts seems justifiable, although they should be used with caution when analysing the posterior rib aspect in lateral impacts.
The validation of a simulation model is an essential step in making it biofidelic and useful. Applying the methods of the current study; rib strain in the cortical bone must be validated in order to predict rib fracture risk using this particular metric. As highlighted in the introduction, rib strain has not been validated in the THUMS AM50 and the GHBMC M50-O models. During the development of the GHBMC M50-0 model, different options of modelling human ribs were considered [Li et al. (2010a) and Li et al. (2010b)]. In this process rib strain predictions were validated for certain ribs harvested from a PMHS. However, the rib models in the final GHBMC M50-O model were not equivalent to the ones in the above mentioned studies, instead they were based on a clinical CT data set with a much coarser resolution [Li et al. (2010c); Gayzik et al. (2010)] in contrast to the clinical CT used in the Li et al. (2010a) and Li et al. (2010b) studies.
In recent years, mesh morphing has become a popular technology for modelling populations of HMBs [Hwang et al. (2016), Zhang et al. (2017a), Zhang et al. (2017b) and Hu et al. (2017)]. The parameters driving ribcage morphing are based on regression models developed by Shi et al. (2014) and Wang et al. (2016),. that are based on analysis of clinical CT data. Thus, overprediction of rib cross section dimensions should be expected. According to the parameter study, the most important rib dimensions for accurate strain prediction in anterior-posterior bending includes rib cross-sectional height, rib arch height and rib rotation. Rib arch height and rib rotation prediction using simple thresholding of clinical CT will be accurate as these aspects are involved in the overall rib geometry. However, the rib cross-sectional height, will generally be overpredicted using this technique. As this measure has been established as the most important parameter for accurate strain prediction, future studies should focus on correcting the cross-sectional height estimates in these regression models. One promising technology, which can be used to remove the rib cortical thickness bias in clinical CT, is cortical bone mapping, which has already been applied in ribs by Holcombe et al. (2018).
A recent study by Poulard et al. (2015) suggested a hierarchical approach to assess the biofidelity of a HBM thorax. Initially, the authors began validation on the organ level, and subsequently gradually increased the validation complexity using point loading, pendulum impacts and tabletop tests. The isolated ribs in the current study have been validated for rib stiffness as well as rib cortical strain, in anterior-posterior rib loadings. The next step is to assemble these ribs in a full HBM and validate the rib strain on the full HBM level, using a similar approach as suggested in Poulard et al. (2015). A validated full HBM model can then be used to predict the risk for rib fractures using the probabilistic framework presented in Forman et al. (2012), for instance. This approach would not be required to rely on further Fig. 12. Left: Thickness distribution around SG2 (lateral strain gauge on the cutaneous side) inRib 5. Right: Distribution of effective plastic strain at peak displacement. The black marking shows the element in which the strain was measured corresponding to SG2.

Table 2
Comparison of average rib cross-sectional height. The published data rows represent weighted averages of the Mohr et al. (2007) and Kindig (2009)  validation of the actual progressive fracture, for example by using element erosion. Rather, the model strain would simply be used to estimate risk for fracture initiation, in a post processing step.

Limitations
One limitation of this study is that it is only based on one loading condition, i.e., anterior-posterior bending. It is very likely that the generic ribs developed in this study would not perform as well in other deformation modes, for example in rib shear or torsion. It is also very likely that other rib parameters will become as important, or more important, than the cross-sectional height when considering other deformation modes.
Another limitation is that although the current study is based on averaged values, some properties are based on rather small data sets. The rib cross-sectional data and the cortical thicknesses are based on just seven male subjects. The material data is based on twelve subjects and the validation data are based on another set of seven males. Although some properties are cross-checked with other data sources it remains unknown if these datasets correspond to the true population averages.

Conclusion
Twelve generic ribs, representing the ribs of an average male, were created based on averaged data from several data sources and publications. Rib stiffness and rib strain predictions have been validated for ribs two through 10, loaded in anterior-posterior bending. These generic FE ribs are suitable for strain based rib fracture risk predictions, when loaded in anterior-posterior bending.
A parameter study showed that the rib cross-sectional height, i.e., the smallest of the cross-sectional dimensions, accounted for most of the strain variance during the anterior-posterior loading of the rib. To accurately model this parameter is important in future studies, for example when morphing of the thorax is in focus.

Declaration of competing interest
The authors whose names are listed immediately below certify that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Supplementary data to this article can be found online at https://doi.org/10.1016/j.jmbbm.2020.103742.

APPENDIX A
The normalised position of the strain gauges in the simulation, measured from the anterior end, is presented for each rib in Table A3. This can be compared to the location of strain gauges in the Kindig (2009) study; average of SG 1/4, s ¼ 28.5% � 2.3%, average of SG2/5, s ¼ 56.9% � 4.6%, and average of SG3/6, s ¼ 78.1% � 2.9%.

Table A3
Normalised position of the six strain gauges, measured from the anterior end.

APPENDIX B
The results of the reanalysed rib material data from the experiments presented in Kemper et al. (2005) and Kemper et al. (2007) are presented in Table B1. The values in each row represents the average of all tests for that subject, while the values in the lowest row represents the overall average.

Table B1
Results of re-analysis of rib material data [Kemper et al. (2005)

APPENDIX C
To quantify how much each of the parameters contribute to the variance, linear regression models for the six strain values were created, see Table C1. Each column represents one regression model, i.e., the regression model for the first strain gauge SG1 is shown in Equation (C1).

Table C2
Each parameter's relative contribution to the strain variance. Each column is one regression model, the first for the force and the six others for the six strains. Relative contribution levels below 0.2 are shaded.

APPENDIX D
The bending moment of inertia (around the weakest axis) for ribs four through seven, see Fig. D1, compared for some HBMs and other published data. Fig. D1. Average bending moment of inertia around the rib 2-2 axis (weak axis). Comparison of different HBM models and measurements of real human ribs. The xaxis is defined from posterior (0.0) to anterior (1.0)..