Development of a finite element biomechanical whole spine model for analyzing lumbar spine loads under caudocephalad acceleration

Background: Spine injury risk due to military conflict is an ongoing concern among defense organizations throughout the world. A better understanding of spine biomechanics could assist in developing protection devices to reduce injuries caused by caudocephalad acceleration (+Gz) in under-body blasts (UBB). Although some finite element (FE) human models have demonstrated reasonable lumbar spine biofidelity, they were either partial spine models or not validated for UBB-type loading modes at the lumbar functional spinal unit (FSU) level, thus limiting their ability to analyze UBB-associated occupant kinematics. Methods: An FE functional representation of the human spine with simplified geometry was developed to study the lumbar spine responses under +Gz loading. Fifty-seven load curves obtained from post mortem human subject experiments were used to optimize the model. Results: The model was cumulatively validated for compression, flexion, extension, and anterior-, posterior-, and lateral-shears of the lumbar spine and flexion and extension of the cervical spine. The thoracic spine was optimized for flexion and compression. The cumulative CORrelation and Analysis (CORA) rating for the lumbar spine was 0.766 and the cervical spine was 0.818; both surpassed the 0.7 objective goal. The model’s element size was confirmed as converged. Conclusions: An FE functional representation of the human spine was developed for +Gz lumbar load analysis. The lumbar and cervical spines were demonstrated to be quantitatively biofidelic to the FSU level for multi-directional loading and bending typically experienced in +Gz loading, filling the capability gap in current models.


Introduction
The North Atlantic Treaty Organization (NATO) identified the cervical and thoraco-lumbar spines, the lower leg and foot/ankle, and non-auditory organs as the primary areas of interest for protection of occupants subjected to under-body blast (UBB) attacks (NATO 2007). This effort is hampered by the lack of an available spine injury research surrogate. The human spine has parallel load paths through the anterior spinal column, consisting of vertebral bodies and intervertebral discs (IVD), and the posterior vertebral structures (figure 1). The extent of the contribution by the posterior structures was quantified by Prasad et al (1974) and estimated at between 28 and 40% of overall lumbar spine load bearing capability (Bogduk 2005). Current physical test surrogates, such the Hybrid III and Warrior Injury Assessment Manikin (WIAMan) anthropomorphic test devices (ATDs), are not capable of representing the two-column response since they are equipped with a single column lumbar spine (Pietsch et al 2016).
FE analysis is a suitable approach for filling the need of a biofidelic spine surrogate. FE models of the lumbar spine have been developed to study various aspects of spine biomechanics such as automotive safety, range-ofmotion (ROM), orthopedic treatment, and disc degeneration (Shirazi-Adl 1994a, 1994b, Chen et al 2001, Ayturk and Puttlitz 2011, Park et al 2013, Xu et al 2016, Xu et al 2017, Amiri et al 2018, Finley et al 2018. These models are, in general, designed for applications related to vibrational, static, or horizontal plane loading, are incomplete for whole body loading analysis and insufficiently validated for caudocephalad acceleration (+Gz) loading. There were some models reported in literatures have limited features in them that are relevant to the current study. An L1-L5 spine was modeled with a discrete mass attached to L1 to represent upper body mass (Zhang et al 2011). The model was exercised using an idealized 20 ms, 102 g peak +Gz to represent UBB loading. Subsequently, a sacrum, pelvis and two femurs were added to the model and it was quasi-statically validated for vertebral body displacements under a 700 N compressive load . This partial spine model was not validated for bending or shear. The Takata human body model (TKHM) was updated with a T12-L5 lumbar spine based on the geometry of the Global Human Body Model Consortium CAD data (Zhang and Zhao 2013). The cervical spine of the TKHM was previously validated (Zhao and Narwani 2007). The lumbar spine model was validated for three separate loading conditions: 3.5 mm of compression, 18 mm of anterior shear, and 9 degrees of flexion. The validation criterion was based on the model's L1-L5 average stiffness responses matching at least one of the two post mortem human subject (PMHS) objective studies (Belwadi andYang 2008, Demetropoulos et al 1998), with matching defined as the simulation average response falling within the PMHS response average plus or minus one standard deviation (SD). The model did not include validations for extension and the validations were whole-lumbar spine so the biofidelity of the functional spinal units (FSU) were unknown. Lei et al (2018) modeled an L1-L5 spine and pelvis to investigate UBBtype +Gz loading. The lumbar spine was validated for force and bending responses under one loading condition. A 30.2 kg mass was attached to L1, both L1 and L5 were constrained to vertical motion and L5 was loaded with an approximate 120 ms, 14.8 g (peak) +Gz pulse. In the PMHS experiments that provided the validation objective responses, T12-L5 spines were preloaded with 5 N m of flexion and the upper mass was asymmetrically applied through a loading cylinder. From the simulation protocol description, it is unclear if or how these boundary conditions were included in the simulation. The peak simulation compressive force (Fz) and peak flexion bending moment at L5 were reported to be within one SD of PMHS experimental responses. The effect on the force and bending responses of the spine with the inclusion of T12, an extra spine level in the experiment versus the simulation, is unknown. The model was not validated for bending or shear loading and since it was a whole-lumbar spine validation, the biofidelity of the FSU, the thoracic and cervical spines were unknown.
The aim of this study was to design and validate an FE whole-spine model (or spine subsystem) that quantitatively matches in vitro experimental PMHS spine force and moment responses, has a posterior loadpath, and is suitable for integration with an FE wholebody model. In order to accomplish this, the designed spine subsystem had to include posterior components representing the apophyseal (facet) joint load-bearing capability. Also, the spine surrogate had to include a biofidelic cervical spine, because data obtained from ejection seat +Gz experiments proved that the head and neck play a critical role in lumbar loads (King and Vulcan 1971). This model is unique because it was validated at both whole-lumbar spine and FSU levels against experimentally obtained stiffness response curves. Validations that match single response characteristics, such as peak load, only capture a single point on the spine's load curve response. With exception of the TKHM, the prior models lacked a validated cervical spine. The significance of this model is that it provides an analytic tool to quantitatively evaluate human body lumbar spine forces and moments, which can be directly compared to experimentally determined spine damage thresholds to predict spine injury.

Methods
The development of the spine subsystem consisted of four primary steps: (1) finite element design and modeling, (2) optimization of material/geometric characteristics, (3) validation against PMHS experimental data and (4) a mesh density sensitivity study to ensure model convergence (figure 2). The convergence study was performed as the last step to ensure convergence was evaluated on the final version of the spine model.   Prasad et al (1974) highlighted the importance of the facet joint when analyzing +Gz lumbar spine loading. Developing representative facet joint models is challenging due to large interpersonal variability. Bogduk (2005) examined incidence of flat versus curved facet joint surface planes from a Horwitz and Smith (1940) study. Of the 284 evaluated facet joints, 48% were flat and 52% curved. Bogduk also analyzed the facet joint plane angles (with respect to sagittal plane) for the lumbar spine from the same study and found high variability. For the L3-L4 joint 12% were 0 degrees, 45% were evenly split at 15 and 30 degrees, 37% were at 40 degrees and 5% at 60 degrees. Because variations were too large, Panjabi et al (1991b) omitted articular facet ligament attachment coordinates from a 3D morphology study.
A study on the biomechanical effects of curved versus flat facet planes determined the joint response was similar for flexion and extension, but deviated by 36% in shear loading (Holzapfel and Stadler 2006). Since off-axis compressive loading will generate a shear component, shear loading can play role in +Gz biomechanics.
Since the purpose of the spine subsystem model was biomechanical biofidelity, the design was not necessarily restricted to an anatomical model. Therefore, both functional and anatomical designs were considered. An upfront concern with choosing an anatomical model was the large interpersonal and intrapersonal variabilities present in the human population. Functional and anatomic modeling methods have their advantages and disadvantages for FE analysis of spine loading. The primary benefits of functional representation, without constraints of representing the detailed geometry, are computationally efficient, simple mesh development and parametric flexibility. The benefits of anatomical representation are increased biofidelity (but only for specimens with similar geometry) and prediction of potential risk of injury based on local stresses and known biological material failure thresholds.
The functional modeling method was chosen for this implementation. Since the model validation objectives involve experiments utilizing specimens from numerous PMHSs, a single anatomical spine model could produce outlier characteristics that prohibit validation under certain loading conditions. In addition, a functional spine with parametrically defined characteristics is well suited for stochastic geometrical spine studies. Even with automated procedures, modifications to anatomical spine geometry typically requires complex morphing as part of simulation preprocessing.
The posture of the new spine subsystem was based on data generated for the WIAMan Soldier Seated Study (SSS) (Reed and Ebert 2013). The spine subsystem was modeled for eventual merging with the Wayne State University Human Body Model (WSUHBM) (Belwadi et al 2012). The WSUHBM thoracic spine was 11 mm longer than SSS thoracic spine. The discrepancy was compensated for by adding 1 mm to each thoracic joint space from the SSS. Simulations were run in explicit mode using LS-DYNA (Livermore Software Technology Corporation, Livermore, CA) FE software (Hallquist 2006).

Geometry
The anterior portion of the functional spine was modeled as a column of elliptical shaped cylinders representing the vertebral bodies and IVDs (figure 3). The vertebral bodies were modeled as elliptical cylinders. The vertebral body cross-sectional profiles were uniform through the height of each vertebra. The depth, width and height of each body was based on measurements from anatomical studies (Panjabi et al 1991a, Panjabi et al 1992, Panjabi et al 1991c. The IVDs were meshed dimensionally continuous with the inferior surface of the superior body and with the superior surface of inferior body with linear tapering in-between the two margins. The height of each disc was dictated by spaces between the vertebral bodies when aligned to the scaled SSS posture. The inferior aspect of the L5-S1 IVD was shaped to mate with the WSUHBM sacrum. The posterior elements were simplified and represented by wedges projecting posteriorly with a slot in the superior surface and a vertical post projecting inferiorly (figure 3). The length of the posterior elements (wedges) was calculated from vertebral measurements by Gilad and Nissan (1985) for the cervical and lumbar spine, while thoracic dimensions were determined from Busscher et al (2010). The length of the posterior elements is meaningful as it defines the supraspinous ligament (SSL) moment arm, as well as providing a representative contact offset for the posterior elements of the spine. The superior surface of the posterior wedge hosts a slot and a post projects in the caudal direction from the inferior wedge surface (C2-T1 and T10-L5 only). The post itself has inferior and anterior protrusions that were faced with 1 mm thick shell elements representing cartilage (figure 3). Both anterior and inferior protrusion lengths can be programmatically adjusted during model preprocessing with LS-PREPOST (Livermore Software Technology Corporation, Livermore, CA) or adjusted parametrically at run-time within LS-DYNA.
The adjustable protrusions provide accommodations for interpersonal variability and interact with the inferior vertebra, providing the contact function of the facet joint. The slot of each vertebra interacts with the post from the superior adjacent vertebra to provide facet joint functions of compressive load bearing and limiting extension, axial rotation, anterior and lateral shear. The anterior protrusion contacts the posterior surface of the inferior vertebral body under sufficient spine extension and/or anterior shear. The inferior protrusion contacts the base of the inferior vertebral slot during compression and extension. The posterior element contacts were defined using CONTACT_ AUTOMATIC_SURFACE_TO_SURFACE methods. The lateral surfaces of the slot (lumbar only) provide a boundary to limit lateral shear and torsion.
The posterior offsets for ligament attachments were derived from a dimensional study of the origin and insertion of lumbar spine ligaments (T12-L5) (Panjabi et al 1991b), with exception of the SSLs, which were attached at the posterior boundary of the posterior wedges. The SSLs and interspinous ligaments (ISL) were positioned in the sagittal plane and the ligamenta flava (LF) and facet joint capsules (JC) were attached at the lateral margins of the posterior wedges. The tensile function of the facet joint capsules was represented by the parallel JC ligaments.

Mesh
The anterior spine, vertebral bodies and IVDs, were modeled with 34 300 nodes and 24 900 solid elements, hexa-or penta-hedra, typically 2-3 mm in size. The superior margin of the IVD attaches to the superior body via a CONSTRAINED_NODAL_RIGID_BODY (CNRB) and inferior margin attaches to inferior body via a TIED_CONTACT. The posterior spine was modeled with 4 000 nodes, 1 650 solid elements (also hexa-or penta-hedra) and 46 shell elements, with typical element sizes from 1-11 mm. The posterior elements were attached to the vertebral body via the same CNRB that attached the body to the inferior IVD. The atlanto-occipital joint (or C0-C1) was defined as a discrete joint. The spinal ligaments were modeled with tension-only seatbelt elements.

Material Properties
The vertebral bodies and posterior elements were initially defined with the material properties of aluminum with a 69 GPa elastic modulus to facilitate potential implementation in an ATD spine. This was deemed acceptable since from biomechanical standpoint, the spine kinematics are dictated by the physical properties of the soft tissues, IVD (Robin et al 1994) and ligaments (Gardner-Morse and Stokes 2004). Considering 'lumped' material properties, that is cortical/cancellous and anulus/nucleus, the elastic modulus of the osseous material is on the order of 50 to 500 times stiffer than the IVD using modulus approximations of 1 GPa and 2-20 MPa, respectively (Lin et al 1978, Stemper et al 2010. In the simulations, the modulus was reduced from 69 GPa to 10 GPa to increase the simulation timestep and throughput by a factor of 2.5 times. The IVD-to-body stiffness ratio was 1:500 (assuming 20 MPa and 10 GPa moduli, respectively), which suggests a 0.2% vertebral body contribution to combined IVD and body compression. However, the body contribution doubles based on a body-to-IVD height ratio of approximately 2:1, resulting in a net vertebral contribution of 0.4% to overall IVD and body compression. The difference with a 69 GPa versus 10 GPa body modulus is the body compression contribution is reduced by a factor of close to seven or down to a 0.06% contribution. Since lowering the vertebral body modulus was expected to change the response variables by less than 0.5%, it made the 250% performance improvement a worthwhile trade-off.
The LS-DYNA material * MAT_PLASTICITY_ COMPRESSION_TENSION (MAT_124) was selected for the IVD material. MAT_124 provides separate stiffness characteristics for compression and tension which was necessary for IVD modeling since IVD have distinctly different compression/tension stiffness characteristics. Stemper et al (2010) determined elastic compression and tension moduli of 17.7 MPa and 2.5 MPa, respectively, in the young male thoracic spine. The key MAT_124 parameters for material response were load curves for compression and tension (LCIDC and LCIDT). The IVD annulus and nucleus were modeled together as homogeneous lumped entities, with the response contributions of the annulus fibrosus and nucleus pulposus approximated by the separate load curves for compression and tension. A reverse engineering method was used to determine material load curve parameters which produced FE model responses that matched the applicable flexion, extension, compression, bending and/or shear PMHS response curves. Nominal elastic moduli, E, of 1 GPa for cervical and lumbar spines and 2 GPa for thoracic were coupled with very low cutoff stresses (8.78 kPa) for parameters PC (pressure in compression) and PT (pressure in tension) to direct the material behavior to the LCIDC and LCIDT load curves at very low strains. Cowper-Symonds strain rate scaling parameters, C and P, were implemented to accommodate varying strain rates experienced in wholebody dynamic simulations (thoracic and lumbar spines only) by scaling the yield stress according to the strain rate, ε′ (equation (1)). The atlanto-occipital joint was defined as a discrete joint with stiffness characteristics to match the PMHS flexion/extension load curves.
Spinal ligaments were modeled as linear elastic seatbelt elements with linear force-deflection properties in tension with no force input under compression. The anterior longitudinal ligament (ALL), posterior longitudinal ligament (PLL), JC, LF, ISL, and SSLs were included in the model while the intertransverse ligaments (ITL) were not represented in the model due to their minimal contribution to posterior stiffness in flexion. The ITL stiffness is on the order of 0.3-1.8 N mm −1 for strains ranging from 0%-23.3% (Rohlmann et al 2006) putting its contribution to combined posterior ligamentous stiffness at a less than 1% of the overall structure. The entire spine utilized the lumbar ligament stiffness coefficients of 33.1, 20.0, 26.2, 30.5, 1.3 and 24.8 N mm −1 for the ALL, PLL, LF, JC, ISL, and SSL, respectively (Yoganandan et al 1988). Those coefficients were derived from quasi-static experiments using a 2.5 mm s −1 displacement rate. In 10 g +Gz whole-body simulations with this spine model merged into the WSUHBM, the average peak ligament strain rates were 2.4±1.9 s −1 for lumbar, 6.8±5.4 s −1 for thoracic and 12.9±14.1 s −1 for cervical spines (evaluated over a 0.8 ms interval from data filtered at channel frequency class (CFC) 60 Hz with an SAE J211 filter). Classifying strain rates of 0.5 s −1 as 'quasi', 20 s −1 as 'medium' and 150 s −1 as 'high' (Mattucci and Cronin 2015) biases these lumbar and thoracic peak strain rates towards a quasi-static rate, mitigating ligament rate insensitivity effects. The ligament's mass per unit length was based on a density (Ashby 2017, 660 pp) of 1 300 kg m −3 and ligament cross-sectional areas described by Eberlein et al (2004). A materials overview of the model components are summarized in table 1.

Objective PMHS studies
The functional spine subsystem characteristic parameters, consisting of material properties (including stiffness curves) and posterior geometric dimensions, were optimized based on data reported in a mix of quasistatic and dynamic PMHS experiments. Loading modes were compression, flexion/extension, and anterior/ posterior/lateral shears for the lumbar spine and flexion/extension for the cervical spine. A requirement for the PMHS experimental studies to be included for model optimization was either the availability of SD or multiple response curves that could be used to calculate the SD. The SD provided the corridors that were necessary for the model validation. The following studies met this inclusion criterion and provided the objective responses for optimizing and validating the spine subsystem.

Lumbar compression
A study comparing the mechanical properties of PMHS lumbar spines versus Hybrid III ATD lumbar spines included force-deflection histories for ten PMHSs (Demetropoulos et al 1998). The aim of the study was to characterize the dynamic properties of PMHS and Hybrid III lumbar spines subjected to tension, compression, anterior shear, posterior shear, left lateral shear, flexion, extension and left lateral bending. The non-destructive tests using T12-L5 lumbar segments were conducted at a rate of 100 mm s −1 . The specimens were potted at T12 and L5. The PMHS compression response curves were averaged to determine the optimization objective response and SD was calculated to provide the validation corridors. Belwadi (2007) tested seven complete PMHS lumbar spines under multiple loading modes. Specimens included T12-S1 whole lumbar spines (LS) and twodisc motion segments: L4-S1 (M1) and T12-L2 (M2). LS specimens were compressed to 3.5 mm of deflection, the two-disc motion segments, M1 and M2, were compressed to 1.75 mm of deflection at a rate of 8 mm s −1 . The compressive force-deflections responses were averaged to determine the objective responses and SDs for optimization and validation. The evaluation of the M2 response was truncated at 0.9 mm due to early subject response drop-out.
In a study to determine stability effects of burst fractures, Panjabi et al (1994a) subjected 13 PMHS thoracolumbar two-disc motion segments (11 T11-L1 and two T12-L2) to 200 N of quasi-static compressive loading. The 'INTACT' response curve provided an optimization objective response along with the defined SDs for validation corridors.

Lumbar flexion/extension
Nine PMHS lumbar spines (five L1-S1 and four L2-S1) were subjected to quasi-static pure flexion/extension bending moments from 2 through 10 N m in 2 N m increments in a study by Panjabi et al (1994b). A 100 N compressive preload was maintained during testing. The average flexion/extension load curves were utilized for model optimization and the SDs along the average curves for model validation. This Panjabi study was chosen over a similar more recent study by Guan et al (2007), which also include statistical corridors, because their maximum loading was limited to 4 N m.

Lumbar anterior/posterior shear
The same Belwadi (2007) study that generated objective lumbar compression load curves was also utilized for anterior/poster shear. The LS specimens were subjected to 18 mm of anterior and posterior shear, while the M1 and M2 specimens were limited to 8 mm. The shear was applied at a rate of 4 mm s −1 . The response curves were averaged and SDs calculated for model optimization and validation. Sundararajan et al (2005) subjected 12 lumbar spine motion segments to quasi-static (0.5 mm s −1 ) lateral shear loading (with 1 000 N compressive preload). Levels tested were eight T12-L1 (reduced to seven as a result of a specimen mounting failure), two L2-L3 and two L4-L5 specimens. The samples were loaded to failure with disc-endplate separation being the most common failure mode. The average displacement at the first sign of injury was 4.7±1.6 mm and the lower bound of the 95% confidence interval for start of injury was 3.5 mm. The researchers determined there were no major differences between the motion segments at different spinal levels, so the study's response data was used as an objective for optimizing and validating T12-S1 under lateral shear loading.

Cervical flexion/extension
The cervical spine has an important role on torso kinematics under +Gz loading. Vulcan et al (1970) determined that head flexion was instrumental in generating torso flexion. Optimization of the cervical spine was limited to flexion/extension since predicting cervical spine loads was not a study objective. Objective responses for cervical spine material optimization were based on flexion/extension load curves to ±1 N m from quasi-static experiments for C0, also known as the occiput (O), through C7 (Panjabi et al 2001) and to ±8 to 10 N m from 500 deg s −1 dynamic experiments for C2 through T1 (Barker et al 2014). Flexion/ extension response curves with ±3.5 N m of loading for O-C2 and C3-T1 were included in the cervical spine validation (Nightingale et al 2002, Nightingale et al 2007. The O-C2 responses were derived by combining the C0-C1 and C1-C2 responses. Although C1-C2 primary ROM is in axial rotation ('no' gesture), that motion was beyond the scope of this study and the joint was optimized for flexion/extension only.

Thoracic flexion/compression
The thoracic spine was not wholly optimized or validated with exception of T11 (compression) and T12, which were included with the lumbar spine model development. The thoracic IVD material characteristics were determined by optimizing compression and flexion responses of the mid-level, T5-T6, IVD. The majority of the thoracic spine has a limited ROM with respect to the cervical spine, largely due to interactions with the rib cage. Supplemental seat belt elements were included in parallel with the SSL to limit the hyper-flexion promoted by +Gz loading to biofidelic response ranges.
Thoracic spine compliance coefficients were determined by Panjabi et al (1976) for T1 through T12, however each coefficient was based on testing of a single specimen. The authors concluded that the stiffness did not vary significantly by spine level. The average stiffnesses reported in the study were 0.76 and 1.25 kN mm −1 for tension and compression, respectively, and 2.22 and 2.82 N m deg −1 for flexion and extension, respectively.
The Panjabi flexion/extension specimens were loaded with a moment of 6 or 9 N m. A 6 to 9 N m moment divided by Panjabi's flexion and extension stiffness coefficients equates to a ROM between 4.8 and 7.2 degrees per spinal level or 48 to 72 degrees for T1 through T10. When compared with a T1-T11 ROM estimate of 52 degrees, (White and Panjabi 1978) suggests Panjabi's flexion/extention experiments covered the bulk of the thoracic spine's fore/aft ROM and Panjabi's flexion coefficient of 2.22 N m deg −1 was a suitable objective for optimizing the FE model spine's flexion response.
A limitation of the Panjabi et al (1976) thoracic study was the peak compressive loads were limited to 160 N, which is a fraction of expected loads under +Gz UBB loading. King and Vulcan (1971) generated force-deflection curves of two-disc functional spinal units (T11-L1 and L1-L3) with dynamic compressive loading of up to 2.2 kN. The force-deflection responses were non-linear. King and Vulcan also reported that a 1.3 kN impact of L1-L3 generated a 0 to 160 N chord stiffness (i.e. force-displacement slope from the origin to the 160 N displacement) of 1.1 kN mm −1 , which was close to Panjabi's 1.25 kN mm −1 . The 2.2 kN impact of L1-L3 generated a higher chord stiffness coefficient of 1.6 kN mm −1 at Panjabi's 160 N peak load point, which was not surprising given the higher loading rate. Since compressive Fz are expected to exceed 2.2 kN, an extrapolated King and Vulcan 2.2 kN response was used to optimize the thoracic T5-T6 IVD material response under compressive loading. A comparison of King and Vulcan (experimental and extrapolated) with Panjabi (extrapolated) is shown in figure 4.
Supplemental LS-DYNA seat belt elements were added to the thoracic posterior elements in parallel with the SSL to function as flexion joint stops for T1-T10 per the combined flexion and extension ROM of 4-6 degrees per motion segment (White and Panjabi 1978). The flexion ROM for each vertebral level was estimated by proportioning the combined ROM by a flexion to extension ratio of 1.4:1 (Panjabi et al 1981). The flexion ROM for T2-T7 IVD joints were limited to approximately 1.5 degrees based on 7.5 N m flexion experiments comparing spines with intact ribcages versus disarticulated ribs (Brasiliense et al 2011).

Optimization
The objective responses (and SDs) were either provided by the PMHS researchers or calculated by the present authors using responses from series of PMHS experiments. All experimental PMHS response curves were digitized using Engauge Digitizer software (Mitchell 2020). If the PMHS study only provided a series of responses, then time-history averages and SDs were calculated from those responses. The SDs were not used in the optimization other than as a visual reference for goodness-of-fit, but were required for the subsequent validation.
The C0-C1 discrete joint characteristics were dictated by joint stiffness constraints at the response curve inflection points described by Panjabi et al (2001). Based on the same PMHS response data, joint stops were included at 10 degrees of flexion and 25 degrees of extension using very high joint stiffness constraints. All other modes of movement and loading were optimized with LS-OPT (LSTC, Livermore, CA). The cervical, thoracic, and lumbar spines were optimized separately.
Separate LS-OPT optimizations were conducted for each spine loading configuration (with exception of flexion and extension conditions which were defined as multi-objective optimizations). Each optimization modeled the same combinations of FSU levels loaded in the same manner as their respective PMHS experiments. Optimization design variables were LCIDC and LCIDT load curve parameters, anterior and inferior protrusion offsets, and the 'C' and 'P' Cowper-Symonds parameters. The lower region of the lumbar compressive load curve was a parameterized second order polynomial of the form y=ax 2 +bx through an abscissa of 0.7, followed by a smooth transition to a very high stiffness to avoid element negative volumes. The cervical compressive load curves were each represented by a separate, linearly scaled variant of the lumbar LCIDC with a shallower stiffness rampup from an abscissa of 0.4. The tensile load curves, LCIDT, were defined as parametric three-point, bi-linear curves. Curves start from the origin (a1= o1=0, where 'a' is abscissa and 'o' is ordinate) to the inflection point at a2, o2, then terminating at a3, o3 (a3 was fixed at 1 or equivalent to 100% strain). The optimized design variables were a2, o2, and o3. The thoracic spine compression and tensile load curves were both parameterized bilinear curves..
The optimization and validation simulations shared the same boundary conditions for each test configuration. The inferior vertebral body was fixed, mid-body, and the superior body was displaced with prescribed motion applied to the nodes on the superior transverse plane of the body. The superior vertebral body displacement was driven by a scaled quarter haversine load curve to avoid instantaneous force spikes at the start of the simulation. Force and moments were measured with a load cell that was coincident with the displacement application. The superior body was detached from cephalic spine at the adjacent IVD junction and superior ligature was removed to avoid artificial loads. The spine levels inferior to the test region were isolated from forces by the fixed inferior body.
To assess displacement rate effects on model response and determine if it was necessary to simulate the anterior and posterior shear experiments in pseudo 'real-time' (i.e. 4.5 s), a loading rate sensitivity study for anterior shear was performed on the lumbar spine (T12-S1) with application of 18 mm of displacement ranging from 100-2 000 ms. The simulations demonstrated a consistent force-displacement response (with the exception of increased amplitude noise at the higher loading rates). The response consistency suggested it was permissible to compress the overall simulation times. The simulation times, displacements and pre-loads, if applicable, are shown in table 2 and 3. The simulations were conducted using motion segments consisting of the vertebrae ranging from the 'Dynamic' level through the 'Fixed' level and any spine segments in between.

Validation
The optimized spine subsystem model was evaluated using CORrelation and Analysis (CORA) 3.61 (pdb, Gaimersheim, Germany). The PMHS experimental responses used to optimize the spine model were also used to validate the model. Simulation forces and moments for all analyses were filtered by LS-OPT at CFC600. CORA recommended default parameters were used to generate the rating with exception of phase, which was weighted to zero since the forcedeflection curves used in this study were monotonic and did not have a distinct phase. The inner and outer corridors were defined as plus or minus one and two SDs, respectively (Barker et al 2014).
There is no formal metric for a CORA rating that describes a 'validated' model. Gehre and Stahlschmidt (2011) have assumed good correlation for a single signal from a simulation if the rating is 'clearly better than 0.8' and 0.7 or higher if a complete test with multiple signals is evaluated. Based on this assessment, a minimum CORA rating of 0.7 was targeted for considering the model as validated.

Mesh density sensitivity study
After the initial model validation was complete, a mesh size study was conducted to determine if the mesh element sizes of the IVD were adequate for model convergence. The study entailed iteratively simulating the lumbar spine subsystem validation series with progressively smaller element sizes and evaluating the responses with CORA after each series. Convergence objective was to have the CORA response ratings within 5% of the previous lower resolution mesh. CORA was chosen as the convergence metric since it utilizes the history of the load curve, rather than a single point on the curve. After the CORA ratings converged, the lower resolution mesh was selected as the final model, to provide improved simulation performance over the higher resolution mesh. A larger mesh size than the initial model was not investigated in order to maintain a minimum of four rows of elements for each lumbar IVD.

Optimization
The final optimization generated LCIDC and LCIDT load curve parameters, anterior and inferior protrusion offsets, and the 'C' and 'P' Cowper-Symonds parameters. The material properties for the IVD are shown in tables 4 and 5. The anterior and inferior protrusions were based on initial post geometry and specific to this model, so their dimensions were omitted; however, the more important optimized gaps were included (table 6). The base FE model units were kg-mm-ms-kN-GPa while other units were derived from this base set.

Validation
The lumbar spine subsystem cumulative CORA rating was 0.766, which was also above the 0.7 objective criteria (table 7), as was the cervical spine overall cumulative CORA rating of 0.818 (table 8). The individual lumbar and cervical spine CORA load case ratings are summarized tables ES-1 and ES-2 and loading response curves in figures ES-1 is available online at stacks.iop.org/BPEX/ 7/015009/mmedia through ES-11 of the electronic supplemental material.

Convergence study
The typical IVD element size of the validated lumbar spine subsystem model was 2 mm. In order to determine if 2 mm element size was too coarse, all the validation simulations were re-run with 1 mm IVD element sizes ( figure 5). Reducing the sizes of the vertebral components was unnecessary due to their high relative stiffness and minimal deformation. For both the 2 mm and 1 mm element size models, the superior margin of the IVDs were attached to the superior bodies via CNRBs. The inferior margins of the 1 mm IVDs also used CNRBs to connect to the inferior bodies since the IVD-to-body node alignment was not conducive to TIED_CONTACTs as was with the 2 mm IVDs. A summary of the CORA response rating differences from the convergence study are shown in table 9. The average CORA rating percentage differences between 1 mm and 2 mm element sizes for all cases combined was 0.77%. All of the load cases with  Panjabi L1-S1-Flexion 0.914 5 Panjabi L1-S1-Extension 0.657 6 Belwadi T12-S1-Anterior Shear 0.802 7 Belwadi T12-S1-Posterior Shear 0.765 8 Sundararajan-Lateral Shear 0.605 Lumbar CORA cumulative rating 0.766

Discussion
The functional representation of the spine described in this study quantitatively demonstrated biofidelic responses over a wide range of loading conditions. The convergences study suggests the lumbar element sizes were near optimal for combined performance and consistency. In general, reducing element size did not substantially effect response. Although the cumulative rating for the model surpassed the overall objective rating, not all of the FSUs reached their individual objectives under all modes of loading.

Extents of model validation
CORA ratings for two of the individual lumbar studies fell below the objective 0.7: the Panjabi extension and the Sundararajan lateral shear. This was not unexpected given there were twenty-seven simulations and six different loading conditions for the lumbar spine. The Panjabi extension series rating was 0.657. Two of the motion segments, L2-L3 and L5-S1, had CORA ratings in the 0.4-0.5 range, which brought the overall series rating below the 0.7 threshold. The FE model relies on the inferior post contact with the inferior wedge to replicate the increased joint stiffness shown in the latter half of the PMHS experimental load  curves. For the L2-L3 and L5-S1 motion segments, the contact occurred several degrees later indicating either a longer post was required or the gap between the post and wedge was too wide at initial configuration. Lengthening the posts on those segments was not fully explored in this study. If gap spacing was the cause, it shows the importance of the initial configuration of the model in motion segment studies. The Sundararajan lateral shear series had a rating of 0.605. The primary reason for the low lateral shear cumulative rating was excessive model stiffness for L1-L2 and L2-L3 joints which rated 0.413 and 0.311, respectively. Both joints departed from the PMHS response corridors at approximately 2.5 mm of lateral displacement which was inside the 0-3.5 mm CORA evaluation interval. The other four vertebral levels, T12-L1, L3-L4, L4-L5, and L5-S1, were above 0.7. To assess the effects the excessive stiffness of the L1-L2 and L2-L3 joints, the peak lateral shear forces of the present study's spine were examined in WSUHBM whole-body simulations at up to 10 g +Gz. The peak observed lateral shear force response was 355 N with less than 1 ms loading duration at that level. In the extreme case of sustained lateral shear, a 355 N load would generate approximately 1 mm of lateral displacement in a spinal motion segment, according to the PMHS response data. Under the same load, the FE spinal motion segments respond with approximately half that displacement. A half-mm displacement differential is a small discrepancy given the overall spine movement, therefore the lower correlation was not expected to adversely affect the model under the whole-body simulated loading conditions. Five of the six cervical studies cumulatively exceeded the 0.7 CORA rating objective, with the sixth having a 0.617 rating. Of the forty cervical loading conditions, no segment fared poorly in more than one rating (table ES-2). The only extension segment to rate poorly was O-C2 (0.377) (table ES-2(f)) compared with the Nightingale response and that was largely due to the wide C0-C1 neutral zone from the Panjabi study. The segments that rated poorly with the Nightingale flexion responses, C3-C4 (0.494) and C5-C6 (0.447) (table ES-2(e)), rated well with the Barker responses, 0.970 and 0.987, respectively (table ES-2 (a)). Likewise, the C4-C5 (0.470), C6-C7 (0.408), and C7-T1 (0.326) segments rated poorly with the Barker responses, but were 0.868, 0.956, and 0.916 with the Nightingale study. The differences in loading rates between Nightingale and Barker likely contributed to these differences.

Mesh convergence study
The mesh convergence objective was a maximum 5% difference between CORA ratings for each load case. There were three load cases that exceeded the 5% threshold. The 1 and 2 mm load response curves for these three 'outliers' are shown in figure 6. The Belwadi T12-L2 compression load curve had a difference of 7.4% (figures 6(a) and (b)). The load responses are very close to each other and closely track the outer (2 SD) corridor. However, the 2 mm element mesh size response falls along and outside the outer corridor sooner than the 1 mm element mesh response, which results in additional zero-weighted scoring. The Panjabi L4 L5 extension load responses were also close (figures 6(c) and (d)), with the 2 mm element mesh size response closer to the outer corridor in the toe region of the load curve which cause its rating to be lower than the 1 mm element mesh size. The Sundararajan L2 L3 lateral shear had the largest differences in load response between the 1 mm and 2 mm element mesh sizes (figures 6(e) and (f)). Since the Sundararajan rating was low, the raw differences in the CORA rating resulted in a larger percentage difference. More importantly, the primary load response difference between the 1 mm and 2 mm element size meshes occurred after 2 mm of lateral shear which is well outside the expected response in the model's current application. As a result of this and minor differences between Belwadi and Panjabi outliers, the mesh was considered converged at a nominal 2 mm element size for the lumbar IVD.

Injury assessment
One of the limitations of functional spine model versus anatomical is the inability to predict injury based on localized strain and known material failure properties. However, a functional model is capable of identifying potential injury by comparing local loads with injury assessment reference values or failure metrics obtained from anatomical injury tolerance studies. Henzel (1967) and Yoganandan et al (2018) describe injury risk curves and fracture thresholds (table 10) for compressive loading. Yoganandan's non-parameterized curves reflect a 5%, 10%, and 50% probability of injury under respective forces of 4 750 N, 5 211 N, and 7 223 N for T12-L1 and 4 710 N, 5 372 N, and 8 545 N for L5-S1. Henzel's data predicts fractures at forces of 7 815 N for T12 and 10 524 N for L5 which are reasonably consistent with Yoganandan's estimates at a higher than 50% confidence level. With regard to axial spine failure, the vertebra is the weak link the load path, typically first failing with the endplate (Adams and Dolan 2011). Therefore, the failure tolerances for IVD discs are unnecessary in prediction of impact-related spine injury. Estimates for flexion spine injury range from 60 N m for the start of ligament damage to 120 N m for gross damage and finally to 140-185 N m for complete failure (Adams and Dolan 1995). Adams and Dolan (1995) estimate the lumbar anterior shear tolerance at approximately 2 000 N due to the resistance of the facet joints. Sundararajan et al (2005) reported an average force of 1 597±654 N at first sign of failure in quasi-static PMHS lumbar lateral shear experiments.

Design assessment
In the course of designing safety systems and comparing designs, differences in model response should provide insight into which design is better. Generally, the design with the lower loads or displacements are better, however, responses in adjacent anatomy should be monitored for unintended consequences. For example, additional restraint to the torso may have desired effect of lower thoracic and lumbar spine loads, but drastically increase cervical spine loads with dire consequences to the occupant.

Conclusions
A functional FE element model of the entire human spine was developed to investigate +Gz spine loading. The spine geometry was compatible with the WSUHBM. Capabilities and limitations include: • Model has parametric geometry adjustment capability to accommodate interpersonal variability.  • The lumbar portion of model was cumulatively validated for compression, flexion, extension, anterior/posterior shear and lateral shears by exceeding the preset 0.7 CORA objective criteria for validation.
• Individually, the lumbar extension and lateral shear series did not meet the objective criteria for validation, but did produce reasonable correlations.
• The cervical spine was validated for flexion and extension only.
• IVD were modeled with a homogenous non-linear material relying on the separate compression and tension load curves to represent the behavior of the combined nucleus and annulus constituents. Separate modeling of the nucleus and annulus may improve the model response.
• The thoracic spine IVD material was derived from optimization of T5-T6 for flexion and compression.
• The model was not validated for lateral bending, torsion, or tension.
• The single post-slot representation of the bony aspects of two facet joints may not respond correctly under asymmetric loading such as torsion or lateral bending, particularly when the spine is in extension.
• Published PMHS injury tolerance studies provide an initial bases for assessing injury potential.

Future work
Future versions of the model could be improved by validating compression and bending of the thoracic spine; compression and shear loading of the cervical spine; lateral bending, torsion, and tension of cervical and lumbar spines. PMHS experimental responses at higher loading rates, particularly compressive, should be added to this optimization matrix to expand the capability envelope to higher g-loading. Develop unified material load curves for the cervical and lumbar IVD materials to improve motion segment consistency across spine level.