A Statistical Model of Spine Shape and Material for Population-Oriented Biomechanical Simulation

In population-oriented ergonomics product design and musculoskeletal kinetics analysis, digital spine models of different shape, pose and material property are in great demand. The purpose of this study was to construct a parameterized finite element spine model with adjustable spine shape and material property. We used statistical shape model approach to learn inter-subject shape distribution from 65 CT images of training subjects. Second order polynomial regression was used to model the age-dependent changes in vertebral material property derived from spatially aligned CT images. Finally, a parametric spine generator was developed to create finite element instances of different shapes and material properties. For quantitative analysis, the generalization ability to emulate spine shapes of different people was evaluated by fitting into 17 test CT images. The median fitting accuracy was 0.8 for Dice coefficient and 0.43 mm for average surface distance. The age-dependent bone density regression curve was also proved to well agree with large population statistics data. Finite element simulation was performed to compare how shape parameters influenced the biomechanics distribution of spine. The proposed parametric finite element whole spine model will assist the design process of new devices and biomechanical simulation towards a wide range of population.


I. INTRODUCTION
Spine controls the movement of the trunk, provides mechanical stability for vital organs and hosts major nerve routes through the spinal canal. It is critical to perform a thorough biomechanical assessment before the product design The associate editor coordinating the review of this manuscript and approving it for publication was Orazio Gambino . for spine structure. Nowadays, biomechanical simulation [1]- [3] has become a standard procedure for spine-related ergonomics product design, musculoskeletal kinetics simulation [4], orthopaedic implant development [5], medical image analysis [6], [7].
Considering the anatomical and physiological variations between population of different ages and genders, increasing studies have investigated spine modelling towards different VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ population groups or even individual person. To conduct comprehensive spine biomechanical analysis, computational spine models with accurate geometry and material properties are indispensable. Although patient-specific models can be obtained via segmenting spinal structure from a tomographic image and determining Young's modulus from CT voxel intensity [8], these subject-specific models are not suitable for population-oriented studies. Hence, there is further need to incorporate characteristics of population into generic application. For population-oriented modelling, it has been more than 20 years since the first construction of the generative spine model. Early works mainly focused on local spine modelling [3], [9], while later studies paid attention to global spine modelling [10]. To the authors' knowledge, generative models incorporating adjustable morphology and material property of the whole spine are still lacking.
To conduct accurate biomechanical simulation, precise modelling of spine geometry and material properties is in huge demand. Therefore, a basic requirement of the generative spine model is to represent the distribution of shape and bone density across the population. For shape modelling, the statistical shape model (SSM) method has been proven to effectively describe shape distribution within a training set of subjects. The concept of SSM was first introduced by Cootes and Taylor [11] and has been widely used for medical image analysis [12]- [14]. SSM uses principal component analysis (PCA) to learn a set of primary shape modes from training subjects. Similar to SSM, statistical appearance model (SAM) [15]- [17] characterizes image appearance (i.e., pixel intensity) distribution based on a set of registered training images. Although SSM and SAM methods have been widely used in the medical image analysis field, few studies have applied these techniques for population-oriented spine modelling. Campbell et al. developed a lumbar model generator using the SSM approach [3], but they did not focus on the entire spine or bone material property. It is well-known that bone material property (e.g., Bone Mineral Density (BMD) and Young's modulus) can be derived from CT pixel value using empirical equations [18], [19]. However, this relationship has not yet been used in SAM to model material property distribution for entire spine.
To provide a generative spine model for populationoriented biomechanical simulation, this study developed parametric finite element (FE) modelling of spine shape and bone material property based on a training set of 65 CT images. The SSM approach was used to model spine shape variations across the population, and nonlinear CT intensity regression was performed to correlate bone density changes with age. Based on this model, the FE models of different shapes and densities can be created for population-oriented biomechanical simulation.

II. METHODS
The pipeline for parametric FE model construction is illustrated in Fig. 1. A mesh-based spine anatomy template was first registered to each training subject to obtain the shape correspondence. The statistical shape model was used to learn shape modes from the aligned training set. For material modelling, voxel correspondence of vertebral images was calculated using image registration, and then age-dependent material property regression curve was derived from corresponding voxel. Finally, the shape and density models were combined to generate different FE mesh instances of adjustable shapes and bone densities.

A. TRAINING DATA DESCRIPTION
To cover a proper range of healthy population, we collected retrospective Positron Emission Tomography/Computed Tomography (PET/CT) images from four central hospitals in the northeast, southeast and central areas of China. The PET/CT images were acquired during the past twenty years for health screening purpose. Among these images, 65 subjects (34 females, 31 males) without spine abnormalities (kyphosis, lordosis and scoliosis) were selected, and their CT images were used for model construction. Table. 1 shows the age and gender distribution of the training data set. The CT image acquisition used 120 kV tube voltage and 28-298 mA current. The pixel sizes and slice spacing ranged between 0.59 to 1.37 mm and 1.25 to 3.00 mm, respectively.

B. ETHICS STATEMENT
This study was performed under the ethical approval from Dalian University of Technology Ethics Committees. No patient identification information has been used in this research or presented in this paper.

C. SHAPE MODEL CONSTRUCTION
To endorse the model with the shape variation feature, we adopted the statistical shape model (SSM) method to learn inter-subject spine shape variations from the training CT images. The spine shape was represented as a triangular surface mesh, and the shape variation was represented as the changes of mesh vertex coordinates.
To build the SSM, a prerequisite step was obtaining the surface point correspondence between the training subjects, as shown in Fig. 1. A commercially available highly-detailed spine mesh model [20] was used as the shape template and was non-rigidly registered to all the training images to obtain surface point correspondence. For template registration, vertebral regions were manually segmented from the CT images  using AnatomySketch Software [21], and converted into triangular meshes using the marching cube algorithm [22]. Key anatomical landmarks were manually defined on each vertebra mesh as illustrated in Fig. 2, and then Landmarkguided Robust Point Matching (LRPM) method [23] was applied to register the shape template to each training subject. In this way, the spine shape of each subject was represented as the surface vertex coordinates of the registered template mesh. To ensure accurate template registration, key anatomical landmarks on the vertebra were defined manually as anatomical constraint. To further refine the template registration, we also conducted landmark-constrained intensitybased registration to obtain voxel-level alignment. This was achieved by filling the registered shape into a volumetric label image (moving image), where the voxels inside each vertebral region were labelled with the same value of the corresponding vertebra in the segmented CT label image (fixed image). Then B-spline transform from the moving image to the fixed image was calculated as an optimization problem. To take VOLUME 9, 2021 advantage of key anatomical landmarks, both image intensity similarity and corresponding landmark distance served as optimized objective. The B-spline transform was solved in a multi-resolution scheme with Gaussian pyramid and Adaptive Stochastic Gradient Descent optimizer using Elastix software [24]. The calculated transform was then used to map the initially registered template to each target shape. Therefore, accurate shape correspondence between all training subjects was established.
Once the correspondence was established, Generalized Procrustes Analysis (GPA) [12] was used to align n corresponding spine shapes to eliminate inter-subject rotations, scale and translations. The spinal shape of i th aligned subject can be represented as spatial coordinates of point set . Afterwards, PCA was performed on the aligned spine shapes based on covariance matrix.
Therefore, it was able to calculate eigenvector ϕ s (shape modes of variation) and eigenvalues λ s (variance) of covariance matrix. The resulting eigenvectors were then sorted to define principal modes according to their corresponding eigenvalues in a descending order. Subsequently, a spine instance V generated from the shape space of the training set can be described as where V is the mean shape of the training set, ϕ s is the s th principal shape mode and α s describes the contribution of s th shape mode. Normally, α s ranging between [−2 √ λ s , 2 √ λ s ] results in anatomically plausible shape variations. To determine the how many shape modes should be retained, common practice defining c is to increase the number of modes until the ratio of the accumulated variance to the total variance ( c s=1 λ s )/( n−1 s=1 λ s ) reaches 0.95 [12].

D. MATERIAL MODEL CONSTRUCTION
As explained in the introduction section, CT voxel intensity can be further converted into material property for biomechanical simulations. In this study, We adopted the empirical equations introduced in [18], [25] to calculate material property (BMD and Young's modulus). The conventional SAM method does not characterize the material distribution with respect to independent physiological parameters directly. Since a major physiological factor affecting bone density change is age, we replaced the PCA step of SAM with nonlinear regression regarding age. Details of material model construction are described below.
To model the age-dependent material distribution of training subject, we needed to determine voxel correspondence at first, as shown in Fig. 1. The intensity-based registration was restricted to the each vertebral region to reduce the computation load. Given the spinal segmentation in the previous subsection, a bounding box enclosing each vertebra was automatically calculated, and then the cropped volume of each vertebra was obtained. A random subject image was chosen as reference image, and then all other training subjects were subsequently registered with reference image using landmark-constrained B-spline registration as mentioned above. To reduce the bias of reference selection, we used an iterative step to identify the training subject whose deformation field was the closest to the average deformation as a new reference. The registration process was repeated until convergence, which resulted in voxel correspondence between the finally selected reference and all other training subjects.
Second order polynomial regression was chosen to model age-dependent material distribution after assessing the coefficient of determination of average material distribution per vertebra explained by linear, second order, and third order polynomial regression. High order polynomial regression produced higher coefficient (R 2 > 0.88) of determination than linear regression (R 2 = 0.76). Significant differences (p < 0.05) were observed between linear and high order regression using F-tests, while no differences were found between second and third order regression. The polynomial function coefficients of each voxel were determined using ordinary least squares regression, where f (t) is material property for a voxel and t is the age of subject in years. This procedure was repeated to model material property distribution for each voxel within vertebra. Considering the gender-specific differences of bone density, we conducted separated regression for each gender and obtained two sets of bone density models.

E. CONSTRUCTION OF PARAMETERIZED FE MESH
Based on the SSM and material model, we built an FE model generator with adjustable shape and material properties. The input parameter of the generator is t, d, α, where t is scalar defining age in years; d is a two-element vector controlling resolution of generated tetrahedral mesh, i.e., the maximum size of the surface triangles and the maximum volume of internal tetrahedrons [26]; α = [α i , i = 1, . . . , c] is a c-element vector of shape mode parameter (illustrated in equation 3). The output of the generator was an FE tetrahedral mesh with user-defined shape (controlled by α), bone material property (controlled by t) and tetrahedral mesh resolution (controlled by d).
The FE model of the whole spine was generated in a vertebra-by-vertebra manner. For each vertebra, an initial tetrahedral mesh of the FE model was created from the reference image using the iso2mesh package [26], which took d as the resolution settings. The generated mesh consisted of fournode tetrahedral elements. The material property of tetrahedral node was derived using trilinear interpolation from reference image. Given the shape correspondence between reference shape and SSM obtained from shape model construction, the initial FE model can be morphed into a shape instance generated from SSM. The morphing was calculated by a 3D Thin-Plate-Spline (TPS) spatial transform using the displacements of corresponding surface points.

III. RESULTS
Since our model was developed for population-oriented biomechanical simulation, it was important to know how well the model characterized the shapes and bone material property of different people. The generalization ability of the shape model was evaluated by fitting the model to different subject CT images. We also compared the age-dependent bone density with existing large population statistics.

A. SHAPE MODEL EVALUATION
To visually assess the shape adjustment ability of the SSM, the values of the SSM shape coefficients α i (i = 1, . . . , c) were adjusted between [−2 √ λ i , 2 √ λ i ] and the corresponding spine shape changes were observed. Fig. 3 illustrates the shape variations corresponding to the first for shape coefficients. α 1 corresponds to the global change of the sagittal spine curvature. The mean shape represents a standard healthy ''S'' shape while increasing and decreasing α 1 causes lordosis and kyphosis respectively. α 2 corresponds to the anterior and posterior bending of cervical vertebrae. α 3 is related to the lateral curvature change of the thoracic spine. α 4 adjusts the lumbar vertebrae dislocation and pelvic flatness in the anterior-posterior direction. These changes are learned from the training set and reflect realistic spine morphometric differences in the population. To quantitatively evaluate the generalization ability of SSM, we performed a CT segmentation using SSM. Accurate and automatic segmentation of specific spine shape is important for conducting successful biomechanical simulation. In this test, 17 spine CT images excluded from training set were selected as test samples. Active shape modelling (ASM) algorithm [11] was used to optimize the position, orientation, scale and shape coefficients α of SSM to fit into the vertebral region. Fig. 4 allows us to visually inspect an example of segmentation results of the cervical, thoracic and lumbar vertebrae. Different pseudo colours were used to represent different vertebrae. Accurate model fitting at vertebral level can be visually observed. In order to quantitatively assess the fitting accuracy, each vertebral region of testing CT images were labelled manually as gold standard. Dice coefficient and average surface distance (D surf ) were used as accuracy measurements.
where Dice coefficient reflects the overlapping ratio between the registered model and the ground truth spine region, R R and R S represent the results of ASM segmentation and the gold standard, respectively. ∩ represents the overlapped area of the segmentation and gold standard, The average symmetric surface distance represents the mean Euclidean distances between surface vertices of the segmentation and gold standard, v R and v S represent the surface vertices of the segmentation and the gold standard, respectively, k and j denote the index value of the corresponding vertex. Fig. 5 shows the box plots of the Dice and D surf . For cervical, thoracic and lumbar vertebrae, it can be seen that the segmentation achieved an accuracy of median Dice>0.7, D surf <0.75 mm for cervical vertebrae and Dice>0.8, D surf <0.5 mm for thoracic and lumbar vertebrae, respectively. The cervical vertebrae present relatively lower median value and greater variability of the Dice coefficient than the thoracic and lumbar vertebrae. This is because the cervical vertebra hs smaller size and greater position variability. As the proposed statistical shape model is registered in whole spine level, the registration algorithm mainly focuses on aligning the larger vertebrae with less pose variation, leading to less accurate registration of the cervical vertebrae.

B. MATERIAL MODEL EVALUATION
To evaluate how well the material model reflects agedependent changes, we generate the voxelized L3 images of different ages in the average shape space. Fig. 6(a) illustrates the axial cross section of the voxelized lumbar vertebra images of different ages. The age-dependent intensity change of cancellous bone region is observed in the section images. It also obviously reflects that the female cancellous bone density becomes lower than the male cancellous bone after the age of 50.
To more quantitatively evaluate our model, we also compared our age-related bone density regression curve with the bone mineral density data collected from the literature. Li et al. [27] used Quantitative Computed Tomography (QCT) to measure the mean BMD of L1 and L2 based on a large Chinese population sample set including 4954 adults (2080 males and 2874 females) between the ages of 20 to 80, and the BMD curve regarding different ages were plotted for each gender. Similar to theirs, we calculated the mean polynomial coefficients of the corresponding cancellous bone region using our material model and plotted the derived BMD curve by polynomial regression and empirical equations for comparison. Fig. 6(b) shows a reasonable agreement between the curves of ours and Li et al [27]. Both their and our curves describe a similar trend of age-dependent BMD changes of both genders. The female BMD falls below the male at the age of female's menopause. It is encouraging to see that our results demonstrate the same menopause age (at which the BMD curves of male and female cross each other) as the results of Li et al., meaning that our regression model agrees well with the large population statistics. It can also be observed that regression curves disagree with the large population statistics at the marginal age groups of 20 and 80, because the corresponding training samples at these ages are relatively small.

C. GENERATION AND SIMULATION OF FE INSTANCES
Based on our FE model generator, whole spine tetrahedral models of varying bone density and tetrahedron resolution were created. Fig. 7(a) illustrates a generated whole spine tetrahedral mesh instance. Fig. 7(b) FE demonstrates the zoomed view of lumbar vertebral instances of different ages and mesh resolution settings. It is noted that the number of nodes and elements in the FE mesh can be user-defined and the age-dependent bone density changes can be observed from sagittal cross section.
To demonstrate the ability of our model for populationoriented spine biomechanical simulation, we conducted three FE model of different shape parameters under compression load. Assuming that the whole spine was in a consistent position, uniformly distributed vertical load 800 N was applied on whole spine, and L5 was considered as fixed support. We tested three FE models (191133 nodes and 964843 elements) by tuning shape parameters and keeping vertebral material property at age 30. All elements were considered as linear elastic material. Young's modulus and Poisson's ratio of intervertebral discs (IVD) were set as 10 MPa and 0.4. The contact setting between vertebrae and IVDs was bonded, and contact setting between facet joints was frictionless.  FE analysis was performed using ANSYS V17.2 (ANSYS, Inc). Fig. 8 and 9 indicate the distribution of minimum principal strain and von-Mises stress of three shape modes. It was observed that peak strain concentrated in T4 area, and peak stress concentrated in posterior part of lumbar spine among different shape modes. It was noted that α 1 = −2 √ λ 1 reflecting excessive kyphosis in thoracic spine caused strain and stress increase in whole spine scale, while α 3 = −2 √ λ 3 reflecting slight thoracic scoliosis caused asymmetric concentration of stress at the concave site in the thoracic spine. The current finding demonstrated gross characteristics that peak strain occurred in the middle thoracic spine under compression load, which coincided with previous study [28].

IV. DISCUSSION
In this study, we constructed a fully parameterized whole spine FE model based on statistical shape model and agedependent material model from 65 subjects of different genders and age groups. A wide range of FE model instances can be generated by adjusting the parameters of shape, age and mesh resolution. The novelty of proposed model will empower population-oriented ergonomic product design and musculoskeletal kinetics analysis.
A merit of the model is the ability to emulate realistic intersubject variations of spine shape. As revealed by the validation results, the model can be used to simulate shape instances excluded from training set. The experimental results demonstrate that individual CT data can be used to guide the creation of personalized spine models. As a result, the model can be used to generate different shape instances either by adjusting the SSM shape coefficients or by fitting into individual CT images. This feature offers good flexibility towards population modelling. The first six shape modes of SSM accounted for over 0.95 variation. Compared with previous SSM studies for vertebra in [3], [13], it was noted that our model required less number of modes to describe a given percent of variation. Besides, we also assessed the generalization ability of our model, which produced lower reconstruction error (ASD=0.47 mm) of lumbar spine than [3] (ASD=2.78 mm) but higher reconstruction error (ASD=1.15 mm) of cervical vertebra than [13] (ASD=0.59 mm). This was because our SSM was constructed for whole spine, and the reconstruction process tended to whole scale fitting rather than each vertebra. Unlike studies reported in [9], [29] which offered parametric shape model with independent anatomical measurements, this study took advantage of shape modes naturally produced by SSM, where multiple anatomical measurements contribute to variation simultaneously.
Regarding the age-dependent bone density modelling, the validation results show that our BMD regression curve reasonably agrees with the large population statistics data and our regression curves accurately reflect the age point where the female BMD begins to fall below the male. Despite the prevalence of BMD in clinics, BMD as single factor has some limitations for bone quality assessment [30]. Hence, there is further need to incorporate material properties and anatomical shape into bone quality analysis. Thanks to the voxel correspondence, the material property of each node in the FE model can be interpolated from the voxel grid. Although this study uses tetrahedral mesh as the FE model structure, such an interpolation strategy is also applicable to other FE model structures as long as the node coordinates are provided. This feature offers the potential to extend our model to different FE model structures.
For FE mesh generation, our generator allows the adjustment of mesh resolution. This feature makes our model suitable for the simulations at different scales for further whole spine scale musculoskeletal kinetics analysis and vertebralevel ergonomic study. Although various studies put effort into how anatomical changes in local spine influenced biomechanics, limited studies have been reported for whole spine FE analysis. Therefore, few data were available for comparison with the proposed model. The results of FE simulation under compression load shows shape modes play an important role in affecting distribution of strain and stress on a whole-spine scale. As a proof of concept, the FE simulation demonstrated the value of our model for population-oriented simulation based on adjustable shape and bone material property.
Limitations in this study need to be addressed. First, we used 65 training subjects to cover population of different ages and genders. Although the model has captured realistic shape and density changes, such an amount of training data is still limited. As revealed by Fig. 6(b), The age-dependent BMD curve deviates from the large population statistics at the age group of 20-30 and 70-80, due to insufficient training data at these ages. We will collect more training data to cover the characteristics of population, including healthy people and diseased patients suffering from spinal deformities or osteoporosis. Second, the images used in this study lacked clinical intensity calibration. We used a pseudo-calibration introduced by [31] to address this limitation. Finally, the FE analysis did not consider anatomy such as ligaments, muscles and ribs around spine. We acknowledge that precise FEM simulation is a complex procedure for clinical situation. The simulation results of this paper proves the value of population-oriented FEM with variable shapes and material properties. More accurate biomechanical simulation incorporating surrounding tissues will be added in the future studies.

V. CONCLUSION
In this study, we built a statistical model of spine shape and bone density based on a training set of CT images. The constructed model learned inter-subject variations of spine shape and material property, and provided adjustable parameters for FE mesh generation. Validation results demonstrated that the model was able to accurately emulate different spine shapes, and offer age-dependent bone density for population-oriented simulation. Our future work will focus on increasing the training data set to more accurately characterize the features of healthy and pathological populations.