FEM Simulation of Non-Progressive Growth from Asymmetric Loading and Vicious Cycle Theory: Scoliosis Study Proof of Concept.

Scoliosis affects about 1-3% of the adolescent population, with 80% of cases being idiopathic. There is currently a lack of understanding regarding the biomechanics of scoliosis, current treatment methods can be further improved with a greater understanding of scoliosis growth patterns. The objective of this study is to develop a finite element model that can respond to loads in a similar fashion as current spine biomechanics models and apply it to scoliosis growth. Using CT images of a non-scoliotic individual, a finite element model of the L3-L4 vertebra was created. By applying asymmetric loading in accordance to the 'vicious cycle' theory and through the use of a growth modulation equation it is possible to determine the amount of growth each region of the vertebra will undergo; therefore predict scoliosis growth over a period of time. This study seeks to demonstrate how improved anatomy can expand researchers current knowledge of scoliosis.

INTRODUCTION

Most commonly found in adolescents, scoliosis affects 2 to 3% of the population with over 80% of the cases being idiopathic [1, 2]. While there have only been two known cases of untreated scoliosis directly resulting in death, the condition can lead to back pains and cardiopulmonary complications or failure [1, 2]. White suggested that scoliosis may be an inherited trait, however no discernable pattern has emerged [1, 3].

There are several different theories regarding the cause of scoliosis that have been suggested. Roaf [4], Sevanstik [5] and Pal [6] believed that scoliosis was the result of asymmetric growth of the spine, with Roaf proposing that the anterior components of the spine grew at a faster rate than the posterior and Sevanstik suggesting that asymmetric growth of the thoracic ribcage was the cause. MacEwan [7] believed that scoliosis was caused by deep nerve damage which induces asymmetric muscle loading while Van der Plaats [3], Azegami [8] and Goto [9] believed that spine buckling due to abnormal or delayed growth of spinal ligaments or muscles was the cause. Regardless of the proposed cause, all researchers believed that an upset in the natural ‘balance’ of the spine resulted in scoliosis spine growth.

Another theory behind the initiation of scoliosis was suggested by Stokes [10, 11], Villemure [12, 13] and Lafortune [14]. Termed by Stokes as the ‘vicious cycle’, it was suggested that asymmetric loading on the spine resulted in differences in bone growth rates within an individual vertebra. This would result in asymmetric bone growth which in turn results in increase asymmetric loading [11, 12, 14]. This cycle continues until halted naturally or through external means.

As there has yet to be a confirmed cause for scoliosis, treatment is difficult. If scoliosis curves prove to be progressive, braces such as the Milwaukee Brace or Boston Brace or at times surgical intervention are employed. Recently, the use of finite element analysis (FEA) to perform predictive spine modelling has become an increasingly valuable tool to improve upon treatment of scoliosis. However, despite the success of many researchers in developing a finite element (FE) scoliosis growth model, many of them lack detailed vertebra geometry and appropriate definition of material properties. Models developed by Stokes, Van der Plaats, Villemure and Lafortune all demonstrated how idealized geometry is often employed to model scoliosis growth. While all four researchers were able to simulate scoliosis growth over a period of several years, individual vertebra were often modelled as either a single beam element, a ‘wireframe’ model arranged in an octagonal format or just as a solid cylinder. The intervertebral disk is also often simplified as spring or beam elements with no distinction between the annulus fibrosis or the nucleus pulposus. In other cases, the intervertebral disk is completely ignored, as studies conducted by Stokes show. Models developed by Azegami [8] and Goto [9], while utilizing radiographic images as a geometric basis, suffered from a low element count that resulted in the loss of geometric detail.

In addition to geometric limitations, current spine models also suffer from the appropriate definition of material properties. Many models, such as the ones developed by Villemure and Azegami assume linear isotropic cortical and cancellous bone properties. Other studies such as Van der Plaats utilized a single bulk modulus of the vertebral bodies. The intervertebral disks was also commonly simplistically modelled and often assigned a single ‘bulk’ modulus representing an average value of the entire disk. These material property allocations can dramatically alter the physiological response under different loading conditions as previously discussed. As all growth models that utilize Stokes’ “vicious cycle” theory (Villemure, Lafortune, Stokes) require the accurate recording of stress values within vertebral bodies if incorrect material properties are assigned to the model, it would reduce the effectiveness of the model in a clinical setting. Current scoliosis growth models, such as the ones developed by Villemure [12, 13], Stokes [10, 11], Azegami [8] and Goto [9] have been able to successfully model scoliosis growth over a period of several years but none of these models utilize physiologically accurate spine geometry or material properties, which can minimize their clinical impact. However, when focusing on the biomechanical behaviour of intervertebral disks or a functional spinal unit (FSU) to loading such as Yao [15], Wong [16], Fagan [17], Lee [18] and Li [19] have done, these models demonstrated improved vertebra geometry and material property modelling compared to the scoliosis models previously mentioned. This study seeks to bridge the gap between accurate spine component models and functional spinal growth models.

The objective of this study is to create a FE scoliosis growth model which utilizes accurate vertebra geometry and accurate material properties in an effort to demonstrate how increased geometric and material property accuracy can improve upon current scoliosis growth knowledge.

MATERIALS AND METHODOLOGY

Geometry Creation and Material Property Assignment

As increased geometric and material property accuracy would significantly increase the complexity of the FE model, it was determine that in order to model scoliosis growth, only a single functional spinal unit would be modelled. To simplify the analysis of the FSU, the posterior elements were neglected from the model. While in some severe cases of scoliosis, the axial rotation of the spine causes growth deformities in the posterior elements of the spine including this in the analysis would introduce another variable set that could not be clearly identified in x-rays.

CT images of a healthy patient were taken with the images used in accordance with ethics proposal (Ethics Project Number R-2777). Healthy patient images were utilized to ensure that any potential scoliosis growth seen to exclusively to loading conditions as opposed to any potential biased geometric effects. Isolating the L3/L4 lumbar vertebra the CT images were transferred into the Simpleware volumetric creation software (Simpleware Ltd, UK). Creating a three-dimensional image of the relevant CT scans, the faceted volume is then transferred into Pro-Engineer (PTC, MA, USA), a modelling software. As one of the principle theories behind the study of Stokes’ vicious cycle theory is asymmetric loading, it was assumed that each vertebra would be symmetric about the midsagittal plane to reduce the effects of geometric differences skewing the results. The annulus epiphysis was not modelled. Model geometry was based upon actual patient data and to generalize the model, surface imperfections unique to the individual were removed, thus simulating an idealized spine and reducing local stress concentration areas due to geometric effects unique to the individual. A constant thickness cortical shell (0.5mm [9, 20, 21]) was also utilized to model the vertebral body as it would create model generality. Endplates, with a thickness of 0.5mm [22], were fitted to the superior and inferior surfaces of the vertebral; and distal surfaces were smoothed to ensure ease of manipulation. Each vertebral body and corresponding endplates were built separately and assembled such that the L3 vertebra was assembled on top of the L4 vertebral body with a maximum anterior spacing between the two vertebrae of 10.5mm.

To model the different layers within the annulus region, the entire intervertebral disk was modelled in Ansys (Ansys Inc, PA, USA). The perimeter between the inferior L3 endplate and the superior L4 endplate was swept to create the outer disk geometry. The swept geometry was then scaled to create eight different annulus layers as well as the nucleus. The nucleus size occupied approximately 50% of the total disk volume [1].

All FSU components were modelled as a transverse isotropic material using 10-node tetrahedral solid elements to allow the elements to accurately mesh the irregular curved surfaces of the vertebral bodies. The annulus fibrosis was modelled a fibre reinforced matrix. The fibres were modelled using rebar elements with linear isotropic properties and orientated at + 30º from the vertical axis and embedded in a linear isotropic ground substance. In order to account for changes in the lamellae layers as they progress towards the centre of the disk, the fibres had varying elastic modulus with the outer two layers possessing the highest elastic modulus and the third and forth layer possessing a slightly lower elastic modulus and so on. The nucleus was modelled as a nearly incompressible material. The z-axis was defined parallel to the superior/inferior plane, the x-axis in the medial/lateral planet and y-axis as the posterior anterior plane. Table provides a summary of the material properties used in the model. A total of 524,916 elements were used in the finite element model. Fig. () displays the FEA model utilized.

Validation of the FSU geometic creation procedure and material property assignment was achieved by following a well defined procedure for FEM [24].

Growth Modulation Equation

The governing growth modulation equation is based on that utilized by Stokes and Villemure:

is defined as the strain rate of the vertebral body, σ is the compressive stress seen by the vertebral body due to mechanical loading, β is the growth sensitivity factor defined as 1.5MPa-1 and subscript m denotes the baseline value seen for healthy, non-scoliotic spines. The baseline strain rate of the vertebral body is calculated from data made available by Stokes [11] and summarized in Table for the three years of evaluation in this study. This equation was applied to only the vertebral bodies themselves. The intervertebral disk does not undergo any growth as a result to loading conditions.

To model scoliosis growth, the strain rates are converted into thermal loads as it can both expand and contract a volume much like bone growth and resorption [9]. The thermal growth equation is defined as.

and

where, ε is thermal strain, ΔTis the thermal load, Δt is the time period, is the strain rate determined from the growth equation above and α is the thermal expansion coefficient (27.5*10-6 mm/oC from [25]). This approach could be applied to all directions, but for this study growth will only be applied along the length of the vertebral body. For this particular study, the time period Δt is defined as one year, resulting in ε to be equal to

In order to implement the growth model, five steps are required to be iterated for each successive ‘year’ of study:

Application of physiological symmetric loading conditions as seen in healthy, non-scoliotic spines;

Application of asymmetric loading conditions as seen in scoliotic spines;

Applying the acquired vertical stress data of each element achieved in steps 1 and 2 into equation 1 as σ and σ respectively. The resulting values are then used in equations 2 and 3 to determine the temperature needed to achieve the required growth measured as strain. The baseline biological growth used in this step varies as outlined in Table ;

Application of the temperature loads on the original FE model; and,

Geometry update to the ‘grown’ FE model.

Loading Conditions

In order to accurately characterize forces applied to the L3 vertebra from gravitational forces, as well as stabilizing muscle loads, and to reduce load application complexity, it was determined that an equivalent load of 614.3N [13, 14], equalling to 0.52MPa of pressure on the superior L3 endplate would be applied. This value is valid for an adolescent with a mass of approximately 45kg at the L3 vertebra level. For asymmetric loading conditions, a minimum of 614.3N is applied to the endplate and a maximum of 637.65N or 0.54MPa of pressure. This maximum value represents the load applied on the healthy L4 vertebra, calculated by Villemure and Nachemson [26] by placing a pressure transducer within the intervertebral disk at the region of interest. The effects of the preceding vertebral body mass i.e. the mass of the L3 vertebra on the L4 vertebra was assumed negligible compared to the applied load. No moment loads were applied as compressive loads would constitute a majority of loading circumstances.

FEA was conducted through the use of Ansys (Ansys Inc, PA, USA) to determine the compressive vertical stresses within each of the vertebral bodies for both load cases. The stress data for each of the elements within the L3 and L4 vertebral bodies is then recorded and passed through a custom Matlab (Mathworks Inc., MA, USA) script that solved the equation (1) through (3) for the strain rate and thermal load for each of the elements in question. The thermal load is then applied to the FE model and the vertebral bodies ‘grown’. For all thermal loading cases, endplate surfaces adjacent to the intervertebral disk are fixed to ensure that the vertebral bodies grew outwards from the disk. Once the newly deformed geometry is saved, the same pressure loads are applied and the process is repeated for years two and three.

DISCUSSION

Scoliosis Growth

As demonstrated in Fig. (), under asymmetric loading the spine begins to display wedging effects, an effect witnessed by studies performed by Villemure and Lafortune. In addition, Fig. () clearly displays how this type of growth behavior is not limited to the sagittal plane and that the asymmetric loading is in fact causing three dimensional deformations. As can be seen, anterior regions of the vertebra appear to grow slightly less than the posterior. This type of three-dimensional growth behaviour is consistent with reported scoliosis growth [1].

As only a FSU was modelled in this study, a Cobb angle calculation is not possible. However, a study conducted by Lafortune analyzed the wedge angle of each vertebra of a scoliosis spine. The wedge angle was measured by finding the angle between a horizontal plane and the superior/inferior endplate surfaces. This measurement was achieved through the use of ImageJ (National Institutes of Health, US). Horizontal lines were first drawn on an image screen shot and then a slanted line that ran parallel to the endplate surfaces. Table displays a summary of the wedge angles calculated while Fig. () compares the wedge angle of this study with ones conducted by Lafortune and Villemure.

As Fig. () shows, the wedge angle of the L3/L4 vertebra in Lafortune’s study translates to about 0.85º after a period of one year. No additional data was provided after a period of one year. Villemure’s model indicates a maximum wedge angle of the entire scoliotic spine after one year is approximately 1.5º and after two years, a maximum wedge angle of 9º. This sudden increase in wedge angle is most likely due to the accumulation of wedging from multiple vertebral bodies as opposed to a single vertebra. In addition, Villemure gave no mention of where the maximum wedge angle was located and most likely reflects the accumulation of multiple vertebra wedge angles as opposed to a single vertebral body. If the cumulative wedge angles for this particular study are taken into account, then this study would see a maximum wedge angle of approximately 1.014º, a value similar to the value seen by Villemure after one year of growth. Additional areas that could cause loading differences can be attributed to slight variations in loading conditions and differences in material properties, vertebra location, geometry and growth sensitivity factor.

Non Progressive Scoliosis Growth

While scoliosis growth is clearly modelled by this spine model after three years, the simulation also suggests that given the current loading conditions, this particular scoliotic spine would be non-progressive. This hypothesis is most evident when viewing the temperature distributions of each vertebral body (Fig. ). The different coloured regions indicate differences in thermal loading with yellow regions having a higher applied thermal load (and thus increased growth) than the green regions. As can be seen from image (a), after one year of growth it is evident that the left portion of the spine will be undergoing greater growth than the right due to the predominant yellow area. When viewing the thermal loads applied in year three however, it is evident that the thermal distribution is becoming more uniform throughout both L3 and L4 vertebras with the green thermal regions progressively getting larger. This indicates that a greater portion of the vertebral body is undergoing slower growth. Should this pattern continue, it can be suggested that at a certain wedge angle the thermal load and therefore the growth rate throughout the entire vertebra will be uniform, thus indicating that the asymmetric load will result in non-progressive growth.

This hypothesis of a non-progressive scoliosis spine is supported by theoretical strength of materials calculations [27] of a curved beam. Referring to Fig. (), it can be seen that for uniform growth to occur the stress distribution along the entire cross section of the curved beam must be equal and constant throughout the entire length of the beam. Consequently, the moment due to the horizontal force (designated as Px) must counteract the force due to the asymmetric load moment (Ma) and the moment due to the vertical force (Py). Recall that the resultant load P is always applied perpendicular to the top surface of the cylinder/vertebra.

As illustrated in Fig. (), as the beam curve or wedge angle (θ1) progressively increases, so do the value of Px and its moment arm while the value of Py and its moment arm decreases. Consequently so long as the loading remains constant, there will be a wedge angle in which the summation of moments due to the applied loads and the moment due to asymmetric loading equals zero at a certain cross sectional area of the cylinder. At this region, only stress due to pure axial compression is present and therefore, the presence of uniform body growth.

To test this hypothesis, the stress values across the entire length of the cylinder (set as the same height as a single vertebral body), seen in Fig. () was calculated using Microsoft Excel©. It was determined that at a wedge angle of 0.511º, the bending stresses at the base (i.e. the fixed end) of the cylinder equaled zero. As this does not necessarily mean that the stresses throughout the vertebral body are also zero, stress calculations must be made along the entire length of the curved beam. These calculations were conducted using identical loading conditions to that of the FE model and assuming that the body is composed of two co-axial cylinders: a smaller cylinder encompassed by a larger one. The applied loading was adjusted to reflect the load seen by only the cortical bone so that the theoretical stress values could be compared to that of the FE model.

Additional calculations determined that at a wedge angle of 1.53º, the right asymmetric growth occurring in the top half of the body will be completely offset by the left asymmetric growth occurring in the bottom portion of the body resulting in zero net asymmetric growth. This switch in asymmetric growth direction occurs slightly below the midpoint of the cylinder as expected since flexural loading will be greater further from the point of application. This is illustrated in Fig. ().

While it could be argued that the lack of posterior elements represents a regression in FEM when compared to models developed by Stokes and Villemure, it is important to recall that none of the previous scoliosis models modelled both cortical and cancellous bones and none were capable of simulating the unique properties of the intervertebral disk. In addition, this model utilized a more accurate physical representation of the vertebral body than the cylinders and wireframe models often employed in scoliosis models.

LIMITATIONS AND FUTURE WORK

As with all FEM, there are several limitations of the modelled developed in this study, the lack of posterior elements being the most obvious. Given that the vertebral bodies are the primary load bearing structures within the spine and that the loading conditions factored in muscle stabilization forces, it is not believed that the presence of posterior elements would change the vertebra body stresses due to the loading conditions, although it may have an effect when vertebra wedge angles increase significantly. It could also be argued that the L3/L4 vertebra does not provide an accurate representation of common scoliotic curves which occur primarily in the thoracic region. While a valid point of contention, it must be noted that the purpose of this study was to demonstrate the effects of accurate geometry and material property assignment would have on scoliosis modelling.

Another model limitation is the modelling of entire vertebral body growth as opposed to just the growth plates. While vertebral body growth does chiefly occur at the growth plates, there is a lack of consensus as to the thickness of the growth plates themselves making it difficult to model. Nevertheless, an analysis of the superior and inferior surfaces of the vertebral bodies in Fig. () clearly demonstrate that even if growth were to only take place at these surfaces the wedging effect of the vertebral body would still take place.

Future studies will work towards the inclusion of posterior elements and the modelling of the entire vertebral column.

CONCLUSIONS

With three years of growth, the FSU clearly demonstrated the presence of vertebra wedging and three dimensional asymmetric growth, characteristic of scoliosis. In addition, the progression of this particular model also suggested that under the current geometric, material property and loading conditions, the scoliotic curve would be non-progressive.

This particular model demonstrated that the use of physiologically accurate material properties and spine geometry in scoliosis finite element analysis was capable of not only modelling scoliosis growth, but also demonstrating non-progressive scoliosis growth. This information can provide valuable insight into the type of loading conditions that would cause progressive and non-progressive scoliosis.

It is important to note that while this particular FEA model is not clinically applicable in its current form, the results shown demonstrates the potential importance of improving scoliosis FE models by implementing accurate geometry and more physiologically similar material properties. Results seen from a single FSU demonstrates that it possesses great potential in predicting scoliosis growth behaviour as well as aiding researchers and clinicians in determining the best course of treatment for those afflicted with this condition and eventually the cause of this condition.

Table 1

Material Properties of the Functional Spinal Unit Components Utilized in this Study

	Elastic Modulus (MPa)			Poisson’s Ratio			Shear Modulus (MPa)
	Exx	Eyy	Ezz	νxy	νxz	νyz	Gxy	Gxz	Gyz
Cortical Bone [23]	11,300	11,300	22,000	0.484	0.203	0.203	3800	5400	5400
Cancellous Bone [23]	140	140	200	0.45	0.315	0.315	48.3	48.3	48.3
Annulus Fibers [21]	550/485/420/360			0.30			-
Ground Substance [21]	4.2			0.45			-
Nucleus Pulposus [19]	4			0.4999			-
Cartilage Endplates [23]	23.8			0.40			-

Table 2

Vertebra Baseline Growth Rates [11]

Age	Spinal Column Height (mm)	Spine Growth Velocity (mm/year)	Strain/year
8	205.9	16.5	0.080311
9	220.9	15.4	0.069927
10	235.1	14.2	0.060333

Table 3

Functional Spinal Unit Wedge Angles Achieved for Each Year of Growth

Year of Growth	Superior Wedge Angles	Inferior Wedge Angles
0	0	0
1	0.635	0.379
2	0.653	0.396
3	0.675	0.407