## Abstract

Low back pain is a major cause of disability and requires the development of new devices to treat pathologies and improve prognosis following surgery. Understanding the effects of new devices on the biomechanics of the spine is crucial in the development of new effective and functional devices. The aim of this study was to develop a preliminary parametric, scalable and anatomically accurate finite-element model of the lumbar spine allowing for the evaluation of the performance of spinal devices. The principal anatomical surfaces of the lumbar spine were first identified, and then accurately fitted from a previous model supplied by S14 Implants (Bordeaux, France). Finally, the reconstructed model was defined according to 17 parameters which are used to scale the model according to patient dimensions. The developed model, available as a toolbox named the lumbar model generator, enables generating a population of models using subject-specific dimensions obtained from data scans or averaged dimensions evaluated from the correlation analysis. This toolbox allows patient-specific assessment, taking into account individual morphological variation. The models have applications in the design process of new devices, evaluating the biomechanics of the spine and helping clinicians when deciding on treatment strategies.

## 1. Introduction

Low back pain (LBP) is one of the main musculoskeletal disorders with major disability effects on the population worldwide [1]. In England, LBP is one of the first causes of activity limitation, sick leave and hospitalization. The related economic effects burden governments, individuals and the society more widely [2] and in the USA, from 1998 to 2008, the healthcare costs have increased from $4.3 to $33.9 billion [3]. Typical non-surgical treatments are muscle relaxation and anti-inflammatory medications, steroid injections, physical therapy and spinal manipulations [4]. Surgical intervention is the last resort, with the most common intervention being spinal fusion [5]. Spinal fusion is associated with prolonged recuperation time, loss of mobility at the fused level and has been shown to increase stress at the adjacent unfused levels [6,7] potentially resulting in degeneration and pseudarthrosis [8]. With these limitations in mind, there is scope to develop new devices to improve the efficacy of the treatment and surgery, and improved methods to assess such new assistive devices.

To meet design and regulatory requirements, a medical device is subjected to a review between each phase of the product design [9]. Each of these reviews introduces additional costs to the device development delaying its final release. New technologies, such as three-dimensional (3D) printing, have assisted with this iterative design process, allowing for the rapid production of cost-effective prototypes. However, the design process of a medical device has to be evaluated in relation with the bodies and tissues with which it interfaces. These interactions can be simulated using finite-element (FE) models and can be combined with iterative optimization of the design topology and mechanical properties.

Lumbar spine models available in the literature are mostly subject-specific models based on magnetic resonance (MR) or computed tomography (CT) imaging [10,11], or models based on averaged approximated dimensions which use relatively simplified geometries [12–16]. While models based on medical images precisely represent the subject-specific geometry, the process to generate these models is both time consuming and expensive. Moreover, to provide wider understanding, beyond a single individual, subject-specific models require a number of models to be solved, for statistical power [17]. However, idealized models based on average dimensions often lack the anatomical detail that is necessary to be of clinical value, with their geometries typically too simplistic. Further, several studies have highlighted the importance of anatomically representative geometry in simulations of the spine [18,19] due to its effects on the intradiscal pressure, the range of motion and facet joint contact forces. Clearly, this is of relevance to any spinal devices developed which use such geometry for the spine.

Recently, the Food and Drug Administration (FDA) and the Medical Device Innovation Consortium (MDIC) focused on the necessity of improving the regulatory system delivering new devices more quickly through the use of computer modelling and simulations [20]. Accordingly, parametric and anatomically accurate models are needed to implement the range of combinations of geometrical features necessary to evaluate the huge variety of clinical cases that can be addressed using the designed device. This methodology would speed up the acceptance of new devices, reducing the risk of failure of the device. Although parametric models have great potential in evaluating a wide range of cases, in order to incorporate them as part of the regulatory workflow, considerable work is needed to verify and validate the results obtained against experimental studies.

The aim of this study was to develop an automated technique to obtain a population of anatomically representative models, which can be used to evaluate the effects of spinal implants on distinct anatomical features of the spine. In this study, a software toolbox named the lumbar model generator (LMG) was developed and implemented using Matlab. Using a parametrized baseline model, the LMG can be used to create a population of geometric models of the lumbar spine (from the L1 to L5 including the intervertebral discs (IVDs)), whose surfaces and solid regions are meshed allowing for direct use in FE models. The parametric model generated by LMG has an anatomically accurate geometry, as evaluated through a comparison with the male dataset of the Visible Human Project (VHP) [21], described in §4.

The models can be reconstructed through the definition of 17 parameters. The parameter set is either determined directly from subject-specific measurements, or can be estimated from correlation analyses based on subject age and height (described in §3). Thus, geometric models are fully parametrized and scalable, so a range of anatomical geometries can be easily generated and replicated. In this paper, the capabilities of the LMG toolbox are described, which includes (i) the methodology for developing the geometry, (ii) the correlation analyses implemented to evaluate the anatomical dimensions, and (iii) the process to obtain the meshed solid model ready to use in FE software.

The innovation of the LMG includes the generation of an accurate geometry of the lumbar spine, based on a set of parameters, and the automatic pre-processing of the solid meshed model. This model is compatible with FE simulations (here briefly introduced but to be included in further publications). Previous studies developed automatic models for the lumbar [18,22] and cervical spine [23,24], based on highly simplified geometries or reconstruction from scans. Comparing these studies with the current study, this toolbox generates models with an accurate geometry based on the definition of only a few parameters. The material properties are easily controlled and implemented in the model to reproduce a model for an individual subject. The FE pre-processing is then performed with the automatic definition of boundary conditions and contact properties and then simulation can be directly run from Matlab, using FEBio (FEBio Software Suite). As far as the authors know, the LMG is the first toolbox which allows the accomplishment of the entire workflow (described in figure 1) from the generation of a geometrical model, the pre-processing of the FE model and then obtaining the solution of the analysis.

## 2. The lumbar model generator toolbox

### 2.1. General features of the toolbox

A LMG software toolbox that can generate the geometry for biomechanical models of the lumbar spine was developed and implemented in Matlab (Matlab^{®}, R2017a, 9.2.0.538062, The MathWorks Inc., Natick, MA, USA). The workflow of the toolbox is shown in figure 1. The toolbox generates a complete lumbar spine model, including the five vertebral bodies (L1–L5) and the four IVDs interposed between them. The main geometric features of the vertebrae and IVD follow recommendations from previous studies [13,25] reported as linear and angular parameters used in the generated model (figure 2). Models are parametrized such that they can be generated using two alternative techniques: (i) based on subject-specific dimensions (which can be directly inputted by the user) directly measured from subject-specific data (e.g. image data and scans) or (ii) using average dimensions derived from correlation analysis based on subject height and age (§3).

Once the geometrical data have been generated, the triangulated surface geometries can be exported to a stereolithography (STL) file, which is compatible with computer-aided design and FE analysis software packages. Further, the surface model can be meshed with solid elements and exported to FE software. Using FEBio to run the FE analysis, the meshed model can be pre-processed. Defining materials, boundary conditions and contact properties through Matlab, the simulations can be requested as an output of the LMG.

### 2.2. Geometrical model

#### 2.2.1. Vertebrae model

An STL file of a lumbar spine (50th percentile) was supplied by an industrial partner (S14 Implants, Pessac, France) used in a previous study [26]. This model was used as a template to reproduce the geometry of the lumbar vertebrae and the IVDs. The two key assumptions were (i) the geometry of the healthy spine is reproduced and (ii) the lateral symmetry, across the mid-sagittal plane, applies to both the vertebrae and IVD [15].

Each vertebra was divided into four regions: the vertebral body, the pedicles, the transverse processes and the spinous process (figure 1*b*). The surfaces characterizing each region were identified and the best fitting polynomial curves were selected to be used to build an accurate and scalable model of the vertebrae. The solid model was parametrized, identifying the dimensions reported in literature for each region, to obtain a scalable model. The parameters identified (listed in figure 2) can be independently scaled according to the dimensions obtained from the subject-specific scans or from the correlation analysis (described further in §3). The following vertebral dimensions have been implemented:

— width of the upper and inferior endplates (EPWu, EPWi);

— depth of the upper and inferior endplates (EPDu, EPDi);

— pedicle height and width (PDH, PDW);

— spinal canal depth and width (SCD, SCW);

— width of the upper and inferior transverse process (TPWu, TPWi);

— spinous process length (SPL);

— vertebral posterior body height (VBHp);

— pedicle sagittal inclination (PDIs);

— pedicle transverse inclination (PDIt);

— intervertebral disc width and depth (IVDw, IVDd);

— intervertebral disc height (IVDh);

— lumbar curvature (

*α*)

#### 2.2.2. Intervertebral disc model

The key variable to describe the disc geometry is the IVD height [18,19], thus only the superior and inferior surfaces have been reconstructed with the best fitting polynomial curves. The height is not constant throughout all the geometry, but it is characterized by different heights in the posterior and anterior aspects [27,28]. However, the average height for the IVD has been used for the LMG, consistent with the literature [29], and the anterior and posterior aspect were linearly scaled. Moreover, the perimeter of each endplate was simplified to be kidney-bean shaped [25,26,27].

The IVD was exported as composed of two different parts, the annulus fibrosus (AF) and the nucleus pulpous (NP). The volumetric percentage (VP) of the two bodies has been reported in the literature [30–32] and a proportion of 44% NP and 56% AF, corresponding to a healthy spine, was used to guide the model [30]. However, these values can be altered to simulate different pathologies. Disc degeneration affects the mechanical behaviour of the IVD, with changes in the composition of the AF and NP, and through structural changes of the disc. Thus, by implementing the volumetric fraction and the height and width of the disc as variables, it is feasible to use the model developed to simulate the mechanical behaviour of degenerated discs [33].

### 2.3. Three-dimensional orientation

The IVDs and vertebral bodies generated from the lumbar generator (§2.1) were arranged in 3D space (figure 3). Four types of lordosis can be identified in the population [34] and in this paper, the 3D orientation typical of the fourth type, in which the five lumbar vertebrae are in lordotic configuration, is implemented. The lumbar sector can be approximated as an arc of circumference, where the angle *α* (figure 2) is the angle between the lines connecting posteriorly the upper surface of the L1 and the lower surface of the L5. In this study, the value of 43.39° was used, based on available literature [35]. The centroid of each vertebral body and IVD bodies were evaluated and used to distribute the bodies in 3D space. The vertebrae have been rotated to follow the lumbar angle (*α*), aligned over the lumbar curvature and spaced using the IVD height at each level. The lumbar curvature is a further input parameter (*α*) which can be varied by the user if required for instance, to simulate pathological conditions.

## 3. Correlation analysis and evaluation of dimensions

The LMG includes a function which enables the automated generation of full vertebrae and discs based on the stature and age of a subject. This follows previous studies which have correlated these two parameters to vertebral dimensions [13,14,25,29,36–38]. According to Jason & Taylor [38], the combined length of the cervical, thoracic and lumbar spine (C–T–L) can be correlated with the stature of a subject:
3.1where *a* and *b* are the correlation coefficients, listed in table 1, between the stature (*h*) of a person and the length of the cervical–thoracic–lumbar segments () of the spine. Then, the posterior height for each vertebral body of the lumbar spine can be evaluated as a percentage of the total length, (table 2) [37].

Likewise, the IVD height is correlated with the age of the subject:
3.2where *h*_{IVD} is the height of the IVD, *n*_{age} the age of a person, and *c* and *d* the correlation coefficients listed in table 3 [39].

The dimensions of the vertebra, described in figure 2 have been correlated in previous studies to the vertebral posterior height (*V*_{BHp}) [13,14] based on datasets originally published by Panjabi [25]. Owing to the axial symmetry hypotheses and to the more complete description of the anatomical features, the correlation analysis used by Breglia [14] has been implemented in the LMG:
3.3where *V*_{BHp} is the vertebral body posterior height, *V*_{par} identify the parameters identified in figure 2 *f* and *g* are correlation coefficients shown in table 4.

The entire sets of equations based on correlation analyses from the literature, described above, were used to generate a full lumbar spine model. Once the age and the height of a patient were defined, these values were automatically evaluated in the script and sent as input parameters ready for use in the LMG (§2), which generates the geometrical model in the order of seconds.

## 4. Accuracy

The accuracy of the model generated through the LMG was assessed through comparison with the male dataset of the VHP [21]. The male dataset was segmented from the axial cryosection photographs of a 38 year old man of height 180.34 cm. The geometry of the lumbar vertebrae and IVDs were segmented from the visible human male dataset into two-dimensional (2D) axial images by professional clinicians [21]. These segmented 2D axial images were then processed in Matlab and were imported into Mimics v. 17.0–v. 19.0 (Materialise, Leuven, Belgium) for reconstruction of three-dimensional (3D) geometry of the lumbar vertebrae and IVDs. Two models were generated, importing the dimensions with the two procedures described above (§1):

The vertebral dimensions were generated following correlation analyses, importing only the age (38 years old) and height (180.34 cm) of the VHP dataset.

The vertebral dimensions were directly measured from a 3D reconstruction of the VHP (tables 5 and 6), and implemented directly into the LMG (i.e. following case 2 in §2.1).

The accuracy of the generated vertebrae and IVD models were evaluated against the VHP model using Cloud Compare (v. 2.9, GPL software, retrieved from http://www.cloudcompare.org/). The Iterative Closest Point algorithm was used to register the STL models and assess the RMS (root mean square) error values. Moreover, the accuracy of the 3D orientation of the generated model, using the VHP dimensions, and the VH model has been examined, evaluating the RMS errors.

## 5. Finite-element model pre-processing

### 5.1. Solid tetrahedral meshing

The geometrical model produced by the LMG is in the form of a point cloud and this section outlines the steps required to import the model into FEA software. In this section, the meshing procedures are described. Once a model has been created and meshed, the LMG allows the user to choose to work either with commercial or open-source software in the model implementation. Alternatively, the model can be prepared using the LMG toolbox and solved directly with FEBio.

The pre-processing of the geometry developed (defined in §§2 and 3) for FEA was performed in two separate steps for the vertebrae and the IVD.

#### 5.1.1. Vertebral bodies

The point cloud of each vertebra was meshed in Matlab, using the GIBBON toolbox [40] and Tetgen [41]. In order to simulate the cortical shell in the vertebral bodies, during the meshing procedure, the internal cancellous core was selected and different element indexes have been assigned. The thickness of the cortical shell can be defined by the user and any value required can be assigned, along with different material properties (figure 4).

#### 5.1.2. Intervertebral discs

The high complexity of the IVD microstructure, where collagen fibres are embedded in the ground substance, required the definition of a structured mesh. A custom algorithm was developed in Matlab, using the Gibbon toolbox, to mesh the AF and the NP. The parameters to define the mesh size are the perimeter points of the AF (pp), the volumetric percentage (VP) of the IVD, the number of layers (nl) which will be implemented in the AF and the number of elements in the axial direction (nz) (figure 5). Initially, a two-dimensional mesh structure was defined, which combined the concentric alignment of the elements and an internal rectangular grid, subsequently smoothed with an elliptical perimeter to improve the mesh quality. The positions of the perimeter points are concentrically scaled and replicated nl times to reproduce the concentric alignment of the collagen fibres. The mesh element size depends on the VP, on the number of layers and on the number of points taken on the perimeter (figure 5). The elements in the AF were oriented following the perimeter and arranged in concentric layers, in order to mimic the layers of fibres of collagen and assign their material properties so as to represent fibre-orientation.

#### 5.1.3. Solid mesh quality

In order to check the mesh quality in the generation of several models, the mesh quality was evaluated through a sensitivity analysis. In the case of the vertebral bodies, two cases were evaluated:

different mesh dimensions were considered (1.2–2.2 mm);

three models were generated based on different person heights (1.75 m, 1.80 m and 1.82 m).

The mesh quality criteria evaluated from TetGen included the aspect ratio, the face angles and the dihedral angles (the definitions are reported in [41]). The evaluation of the mesh quality does not have a unique definition, but it depends on the application [41,42]. Consistent with previous studies, the mesh quality was judged satisfactory when at least 95% of the elements had an aspect ratio less than 4, and face angle less than 100° and dihedral angles were less than 130° [42–45]. A mesh convergence test was performed in FEBio for the L1 vertebra at different mesh sizes (1.4 mm, 1.6 mm, 1.8 mm, 2.0 mm, 2.2 mm), applying 1 MPa pressure over the superior surface. The difference in the stress between the model with the finest mesh and the others was evaluated and the results led to the selection of 1.6 mm with an error of less than 5% (as compared to subsequent mesh refinement). The frequency distribution of the quality criteria of the L1 vertebrae was then compared either with the distribution of the other vertebrae (L2–L5) or the distribution of the models based on different body dimensions, through a quantile–quantile plot (QQ-plot).

The IVD mesh was evaluated through a sensitivity study on varying the parameters which influence the mesh size pp (64, 72, 80, 88, 96), nl (8, 9, 10, 11, 12) and nz (6, 8, 10, 12, 14), described above. The influence of the VP has not been reported in this study because it is not a parameter relevant to the structural properties of the IVD. A pressure of 0.5 MPa has been applied to the upper surface with the lower surface fully constrained, further, neo-Hookean material properties were assigned to both the AF and NP (*E*_{AF} = 5 MPa, *ν*_{AF} = 0.3; *E*_{NP} = 3 MPa, *ν*_{NP} = 0.3, where *E*_{AF} and *E*_{NP} are the Young's modulus and *v*_{AF} and *v*_{NP}, the Poisson's ratio for the AF and NP [11] automatically converted into Lamé coefficients by FEBio).

## 6. Model evaluation

### 6.1. Comparison between the lumbar model generator generated model and the Visible Human Project model

A comparison between the VHP spine and the generated models are shown in figure 6, and the RMS values are reported in table 7. The quantitative analysis highlighted the areas of the model which differ most from the VHP model. In particular, the largest differences have been identified in the posterior vertebral structures, where the superior and the inferior articular facets, lamina and pedicles of the generated geometry are less detailed. For all the vertebrae, the maximum RMS error for the models generated with the correlation analysis was higher than for the models generated using the VHP dimensions. In both cases, the L5 vertebrae geometry showed the highest RMS values. Lower RMS values are shown when the dimensions are directly measured on the subject-specific model and then imported into the algorithm. In fact, using the correlation analysis method based on previous studies [13,14], the model is affected by the grade of correlation between the variables. High and moderate errors have been shown in the dimensions, which had low (*R*^{2} < 0.5: PDW, SCW, SCD, TPW, PDIt, PDIs) and moderate correlation coefficients (*R*^{2} between 0.5 and 0.8: SPL, PDH), provided by Breglia [14]. Figure 7 shows the accuracy of the 3D orientation between the generated model and the VH model. Owing to the supine position of the cadaveric specimen during the image acquisition, the curvature of the lumbar curve is reduced as compared to an upright position. Nevertheless, using CloudCompare a best fit registration has been performed, obtaining a mean RMS value of 2.73 mm.

### 6.2. Mesh quality evaluation

Assessment of mesh convergence demonstrated a suitable compromise between a mesh size (which preserves geometric precision) and variation in predicted results. In the case of vertebral bodies, the differences in stress from the finer model are less than 5%, where the fine model has a mesh size of 1.2 mm. The mesh convergence study on the IVD found that the differences in the stress values were less than 6% varying the nl, nz and pp parameters (figure 8).

A sensitivity analysis was performed to evaluate potential differences in the mesh quality throughout the vertebrae (L1–L5) or on varying the dimensions selecting different body heights of a patient. The mesh quality obtained in both approaches was checked through a QQ-plot where the vertebra L1 with element size of 1.6 mm, generated as average model based on a person of 1.75 m height, was taken as reference. The QQ-plots in figure 9 demonstrate that there were no differences between the distributions on varying either the vertebrae considered (figure 9*a*) and the dimensions of the specimens (figure 9*b*). The sensitivity analysis on varying the geometrical dimensions, in both the IVD and vertebrae bodies, showed that there were no effects on the mesh quality.

The mesh quality reported that more than 95% of elements had an aspect ratio of less than 4, more than 98% of elements had a face angle less than 100° and more than 95% of elements had a dihedral angle less than 130°, consistent with other studies [39]. Quantitative results of the distribution of the dihedral angle are shown in figure 10, where the two histograms show the distribution of the maximum and minimum dihedral on the total number of elements.

## 7. Discussion

The LMG is a toolbox which allows the generation and simulation of FE models for the lumbar segment. It is possible to obtain averaged models, imported using as input only the height, age and gender of a patient, where the dimensions are evaluated from a subroutine based on the correlation analysis described above. Another option is to input the data evaluated from subject-specific scans, or alternatively opting for a hybrid use of the tool, changing the dimensions desired from the averaged model. In the case of the former, it is not necessary to have clinical data available and a population of models can be generated based on the gender, age and height of patients. The output can then be pre-processed in commercial software able to read STL file or input file for commercial software can be requested (e.g. Abaqus, Hyperworks). The main novelty is the possibility to obtain from the LMG, using as input a set of dimensions: (i) a geometrical model, based on dimensions obtained from subject-specific scans or on correlation analysis; (ii) obtain a solid meshed model; (iii) pre-process the model directly in Matlab and (iv) run the simulations using FEBio.

The LMG could be seen as having two end-user applications: (i) clinical/industrial applications where a device is assessed using FEA on a range of models; (ii) engineering science applications where multiple variables are assessed to determine the sensitivity of whole spine mechanics to these variables. Pathological conditions can be evaluated implementing models with different lumbar curvatures to compare to the averaged lumbar curvature, or changing the volumetric ratio between the NP and AF or varying their material properties [46].

Considering that the LMG was meant to be an averaged model, which cannot take into account the subject-specific anatomical variations, the models generated through the LMG toolbox are anatomically realistic. Moreover, subject-specific models are affected by errors due to the segmentation procedure. In fact, the accuracy of the geometry obtained from scans depends on the type of scan used (CT or MR), their respective resolutions (inter-slice resolution ranging between 0.6 and 2.5 mm), the technique used for the segmentation and the operator expertise [47,48]. Hence, the RMS errors obtained in this study are comparable to the range of accuracy of subject-specific models and is thus considered acceptable. The highest differences between the VPH model and the generated one have been highlighted in the accuracy test and they are principally localized in the posterior area of the vertebrae. The effect of the geometry of the spinous processes on the biomechanics of the spine has been evaluated [18,19]; however, no direct influence has been identified. The caveat being the evaluation of interspinous devices, for which this model may have great limitations. The limitations due to the simplified geometry of spinous processes will be investigated further in future work through FE analysis. The RMS error values showed good agreement with the subject-specific model, even though greater difference has been highlighted at the L5 level. This is justified by assessing previous studies [13,49,50], which did not consider the L5 vertebrae in morphometric studies due to the high variability of its geometry. Moreover, the higher RMS values are related to the lumbar vertebrae obtained from the correlation analysis, which are affected by uncertainties. In fact, the correlation analyses performed by Breglia [14] and Kunkel [13] are both based on the Panjabi dataset [25], which evaluated 12 specimens and larger studies are required to obtain improved statistics. A brief preliminary correlation analysis was evaluated with the data collected by Wolf *et al.* [51] and Alam *et al.* [52] and compared with the previous correlation analyses. Evaluating these correlations on other morphological studies, conducted by Alam *et al.* (*n* = 55) and Wolf *et al.* (*n* = 49), it is possible to obtain higher values for the *R*^{2} values for regression analysis in certain cases. Unfortunately, the measurements between the different studies could not be merged because those studies did not consider the same anatomical parameters, gender or ethnic group, which would affect the whole correlation analysis [52,53]. Further studies are required to obtain more complete datasets of morphological measurements and to evaluate new correlation analyses. This preliminary model is to be developed more, the fitting of the posterior part of the vertebrae, the connection between transverse process, lamina and facet joints shapes, shall be improved and further parameters (i.e. thickness and height of the spinous processes) will be added to better describe their geometry.

Furthermore, the results on the mesh quality and its reproducibility on different geometries showed that the automated generation leads to valid meshed models. Using a fixed mesh size, according to the mesh convergence test, it showed that the mesh quality of the vertebrae was not sensitive to the variation of the vertebral dimensions. Hence, the FE model can be directly set up to run the simulations without other time-consuming actions.

Several models have been developed based on subject-specific datasets [54–57]. Subject-specific models have the advantage to reproduce the anatomy of the patient accurately, with the possibility to include soft and hard tissue with high resolution. However, they reproduce the anatomy of an individual subject and reconstructing these individual models is time expensive. Further, they require extensive pre-processing. This final point is actually a barrier to clinical implementation because of the need to have someone dedicated to develop and solve the models. In reply to these issues, the LMG offers anatomically representative models, scalable to average human dimensions, which do not require extensive segmentation processes. Moreover, the pre-processing is accelerated further because the models can be directly meshed and the material properties, boundary and loading conditions can be imported in order to obtain a ready-to-solve FE model.

In the last decades several models have been created using average dimensions [15,16,18,58] obtaining models to get quantitative analysis over the biomechanics of the lumbar spine. However, all of them have used approximations in the anatomy of the vertebrae. Campbell *et al*. [59] developed an automatic tool to reconstruct models from data scans without user intervention and with low computational cost. This enables highly detailed models from subject-specific data. However, the limitation of this is that it requires having data from clinical studies. As far as the authors know, our current study is the first parametric model which generates a full FE model and includes the potential for anatomical variability and the flexibility to input the dimensions of an individual subject. The advantages of this tool in the short term are that the models could be used to perform sensitivity analysis, varying the dimensions most subjected to change in the population. These models could be used in the design process of new devices, or to develop custom made implants and assess their performance, as recently recommended by the FDA and MDIC [20]. The biomechanics of the spine has been identified as being sensitive to parameters such as the lumbar curvature, the vertebral body height, IVD height and the width of transverse processes [60,61]. Therefore, the evaluation of the functionalities of new devices and how they influence the biomechanics of the spine, in several anatomical configurations, will lead to the design optimization and customization of new implants. The intention is to implement the lumbar model to assess the validity of devices such as BDyn, a posterior stabilization device produced by S14 Implants and described elsewhere [62], to assess the performance of the device and its effect on the spine. Further optimization of posterior stabilization devices is of value because, unlike fusion, such devices retain motion at the spinal segment of interest.

The current study has focused on the concept of developing an anatomically accurate but automated model. Future releases of the toolbox are planned to include the definition of the 1D nonlinear spring elements to simulate the action of ligaments [54,63], which will be placed in the corresponding anatomical location. The facet cartilage will be included, defined as a 1D element in between the facet joints and the definition of the facet contact properties will be automatically implemented. The toolbox has been released on Zenodo [64] and the development version is freely available on GitHub. The future release will be provided with a complete GUI (Graphic User Interface) with more complete pre-processing features which allows implementing further loading conditions and exporting the file ready to solve if using other solvers. In future studies, FEA will be evaluated to assess the influence of the anatomical features (i.e. the VP ratio and the vertebral dimensions) on the biomechanics of the lumbar spine and how it is affected using devices such as BDyn and disc replacement devices.

## 8. Conclusion

The LMG toolbox has been developed with the intent of helping in the design optimization of spinal devices such as posterior stabilization devices as well as the development of custom devices. The developed toolbox enables an automated workflow which is user independent and fully compatible with open-source software (Octave, FEBio, Calculix). It constitutes a tool that can be used in clinical studies to improve the decision-making process to select the best intervention. The clinicians, supported by experts, using the GUI would be able to simulate and understand the effect on the biomechanics of the specific patient taking into account the anatomical variation. In fact, it will enable the use of average dimensions or importing the dimensions measured from data scans or evaluating a hybrid model where the dimensions of a desired structure can be altered, then evaluating the specific case to treat. In the short term, this toolbox will be released using an online platform with a user-friendly GUI, which will allow the quick evaluation of the biomechanics of specific cases. It can then be used to aid clinician's practice when assessing the biomechanics of the spine of their patients, and it could lead to improve the decision-making process to select the best intervention.

## Data accessibility

The LMG toolbox is freely available on GitHub [64].

## Authors' contributions

C.E.L. conceived the presented idea and developed the toolbox, including the geometrical development of the model and validation, participating in the study design and drafting the manuscript. K.M.M. participated in developing the custom meshing algorithm used; P.V.S.L., K.M.T. and D.R. participated in the study design related to the VHP model, D.M.E. and D.E.T.S. participated in the study design and supervised the project. All the authors contributed in revising the manuscript.

## Competing interests

There are no conflicts of interest to declare.

## Funding

This study was supported by the European Commission under the 7th Framework Programme (grant agreement no. 604935).

- Received November 3, 2017.
- Accepted December 1, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.