In-silico model development and validation of the L5-S1 spinal unit

Abstract

The L5-S1 segment of the spine is highly susceptible to injury, frequently causing low back pain. The segment has gained a lot of scientific interest, leading to many experimental works that can be found describing its biomechanical characteristics. But, there is a lack of work focusing on its computational studies, which can significantly aid its further studies. In the current study, a subject-specific single-segment finite element model of the L5-S1 unit was developed from a T2-mapped MRI scan. This study is mainly intended to probe the requirements for modelling the annulus of the disc and also attempts to understand the role of ligaments exclusive to the L5-S1 spinal unit to establish its validated finite element model. The annulus was represented by two different forms of hyperelastic material models (isotropic and anisotropic) for which the constants were determined from experimental data found in the literature. Their ability to impart the required characteristic was tested for the finite element model to mimic the experimental responses during sagittal and lateral moment loads. A comparison of results with the two material models is also discussed for other valuable parameters like contact pressure at the facets, maximum von-Mises stresses in the vertebrae, ligament strains, and midplane Tresca shear stresses of the annulus. The anisotropic Gasser-Ogden-Holzapfel (GOH) model was observed to deliver a response that consistently showed good compliance with the experimental response and hence, it is recommended for the computational studies of this segment.

Keywords:

• L5-S1 spinal unit
• finite element method
• Gasser-Ogden-Holzapfel
• Ogden
• T2-mapped MRI

1. Introduction

The L5-S1 spinal unit has a lot of clinical significance, as the common location for the initiation of low back pain is usually the annulus of the intervertebral disc corresponding to this segment and the L4-L5 segment (Schwarzer et al., 1995). Some of the innate features associated with the vertebrae, ligaments, and the disc in this segment along with its location in the human body contribute to significant differences in its biomechanical characteristics compared to the rest of the sections of the spine (Bogduk, 2005). Apart from the occurrence of diseases like herniation or disc prolapse, disc degeneration, disc disruptions, and stenosis which are frequently observed in the L3-S1 region of the spine (Adams et al., 2015; Bogduk, 2005; Graham, 2018; How et al., 2015; Sehgal, 2000) the L5-S1 segment is susceptible to other severe diseases like spondylolysis and spondylolisthesis (Deng et al., 2017) increasing the chances of leading to a traumatic condition called sciatica (Maheshwaran S et al., 1995).

In-silico studies of this segment can be useful to increase the understanding of the underlying causes of the problems that occur in this region. Besides being a cost-effective approach, it provides us the capability to perform tests beyond the physical constraints which may occur in-vivo or in-vitro experiments. It can also prove efficient in testing the efficiency of surgical procedures or other treatment modalities without experimenting on living specimens. Most of the single segment finite element studies have been performed over the L4-L5 and above spinal units (Fasser et al., 2022; Manickam & Roy, 2022; Sengul et al., 2021; Vinyas et al., 2020; Zhu et al., 2022). The amount of exploration on the L5-S1 segment in-silico models and its further improvement are considerably inadequate both in single-segment as well as multi-segment models (Vinyas et al., 2020). The limited experimental studies on this segment to support the in-silico studies are majorly responsible for this situation.

Nevertheless, to our knowledge the immediate predecessor of the current study is the work carried out by Charriere et al. (2003) which is recorded as the first attempt to develop and validate the in-silico L5-S1 single segment model. The work considered geometrical approximations for modelling the disc and explicit modelling of ground substance and fibres of the annulus. The hyperelastic law that was used for the ground substance of the annulus was incapable of capturing its behavior in both the tensile and compressive regime. Also, the validation was carried out at a higher load range where the non-linear behavior of the spine structure is not so evident as in lower load (White & Panjabi, 1990). The current study is highly motivated and focused towards the development of single segment L5-S1 in-silico model and evaluate the capability of two different hyperelastic material models to represent the annulus of the intervertebral disc which is hypothesized to have a better potential than the previously used approach.

2. Material and method

2.1. Hardware and software used

The source data was obtained from GE HDxt Signa 1.5 T MRI scanning machine. The modelling and analysis were performed on a Lenovo Desktop PC with Intel® Xeon® W-2155 CPU at 3.3 Ghz processor, 10 cores, 20 logical processors, 32 GB RAM, and Nvidia graphics of 2GB. MIMICS v 15.0 (Materialise Inc.). Hypermesh v14.0 (Altair Engineering Inc., USA) was used for finite element mesh generation and enhancement of mesh quality. Further preprocessing and simulations were performed on ABAQUS 2017 software. Origin 2017 (OriginLab Corp., USA) was used for data extraction from graphs of experimental studies and for plotting graphs. MATLAB R2021 was used for the determination of hyperelastic material constants.

2.2. Source of geometry

CT scan is the most commonly used source for finite element model creation related to the spine. As CT does not provide any data about the soft tissues, it leads to the limitation of modelling the intervertebral disc on the basis of generalized anthropometric data which might not be valid for every subject. As the current study is focused on the intervertebral disc, using a T2-mapped MRI as the source was found to be more appropriate to accurately model the disc. The specified MRI scan provided the anthropometric data of the vertebra as well as the annulus and nucleus of the intervertebral disc. The parameters of MRI were altered (Vinyas et al., 2022) to make the data compatible with the image processing software, MIMICS v 15.0. The ligaments were modelled using anatomical points as suggested in the literature.

2.3. Mesh creation

The cortical shell and the cancellous core of both the vertebra were assigned the 4-noded tetrahedral elements. The thickness of the cortical shell was considered to be 1 mm (Sadegh Naserkhaki et al., 2016) which was implemented by keeping the edge length of the element as 1 mm. The elements for the cancellous core were set in such a way that their size was the same as the elements of the cortical shell near the interface between them and it gradually increased on approaching the center of the vertebra. Ligaments were represented with 3D Truss elements. The attachment points for the ligaments were chosen by the description of the anatomical markers as suggested in the literature. The articular facet between the two vertebrae was separated by a gap of 0.5 mm to avoid the initial overclosure of elements which adversely affects the analysis. The normal contact property was defined as hard contact. Figure 1 shows the in-silico model of the L5-S1 segment.

Figure 1. Different views of L5-S1 in-silico model, (a) Sagittal view, (b) Posteriolateral view, (c) Anteriolateral view, Transparent vertebrae, and opaque IVD, and ligaments in-silico model of L5-S1 segment (d) Posteriolateral view and (e)Anteriolateral view.

Multiple methods have been proposed and implemented by various researchers in the past to represent the annulus of the disc in finite element models (Vinyas et al., 2020). They can be broadly classified into two main categories, which are: (i) Type-A models, with annulus ground substance and the fibers defined by exclusive elements, and (ii) Type-B models, with ground substance and the fibers defined in an element as a continuum. Most of the studies use the type-A models, where the annulus fibres are modeled with truss elements (Goel et al., 1995; Little et al., 2007; Polikeit et al., 2003; Ramakrishna et al., 2018a; Rao & Dumas, 1991; Shirazi-Adl et al., 1984; Weisse et al., 2012; Xie et al., 2017; Zhang et al., 2020) or surface rebars (Łodygowski et al., 2005; Natarajan, 2006) reinforcing the solid elements. A few studies incorporated the type-B models (Guan et al., 2006; Momeni Shahraki et al., 2015; Moramarco et al., 2010) where the material model includes the effects of both the entities. The current study uses the type-B modelling approach for the annulus of the disc. The annulus and the nucleus of the disc were meshed using 8-noded hexahedral elements as this type of element is found to be more stable, less influenced by mesh size, and provide better accuracy in predictions. Mesh sensitivity was checked for the disc based on the maximum displacement, maximum Tresca stresses and intervertebral disc pressure. The convergence was achieved for 8 layers of elements in both radial and axial directions, since there was less than 5% variation (Jones & Wilcox, 2008) in the values for the selected parameters on further refinements.

2.4. Material models for parts other than the annulus

The vertebrae are almost rigid structures compared to the rest of the parts in the unit. Thus, the cortical bone and cancellous bone were assigned isotropic linear elastic properties. The ligaments play a crucial role in regulating the biomechanical movement of the in-silico model. They actively resist to deformation in a highly non-linear nature only to a tensile load and slacken during compression. This non-linear behavior of the ligaments was implemented using hypoelastic material model. The nucleus was considered an incompressible fluid with a pressure of 1 MPa. The details of all the properties are listed in Table 1.

Table 1. Material properties and element types used of the parts involved in the in-silico model

Part name	Element type	Young’s Modulus (MPa)	Poisson’s ratio
Cortical (S Naserkhaki et al., 2018)	C3D4	12,000	0.3
Cancellous (S Naserkhaki et al., 2018)	C3D4	200	0.3
Nucleus (Guan et al., 2006)	C3D8H	1	0.499
Annulus	C3D8H	Material property
Annulus in SMiso		Ogden hyperelastic material model constants ${\mu }_1$= 0.0575 MPa and ${\alpha }_1$=37.9622
Annulus in SManiso		GOH hyperelastic material model constants C1 = 0.4323 MPa, $\mathit{\kappa }=0$, k1 = 2.1638 MPa, k2 = 200.588, D = 0
Ligaments (Momeni Shahraki et al., 2015)	3D Truss element	Hypoelastic constitutive law (Tension only), Elastic Modulus in MPa (Strain %)
Anterior Longitudinal Ligament (ALL)		7.8 ($\mathit{\epsilon }<12\mathrm{\%})$, 20 ($\mathit{\epsilon }>12\mathrm{\%})$
Posterior Longitudinal Ligament (PLL)		10 ($\mathit{\epsilon }<10\mathrm{\%})$, 20 ($\mathit{\epsilon }>10\mathrm{\%})$
Ligamentum Flava (LF)		15 ($\mathit{\epsilon }<6\mathrm{\%})$, 19.5 ($\mathit{\epsilon }>6\mathrm{\%})$
Interspinous Ligament (ISL)		10 ($\mathit{\epsilon }<15\mathrm{\%})$, 11.6 ($\mathit{\epsilon }>15\mathrm{\%})$
Supraspinous Ligament (SSL)		8 ($\mathit{\epsilon }<20\mathrm{\%})$, 15 ($\mathit{\epsilon }>20\mathrm{\%})$
Capsular Ligament (FC)		7.5 ($\mathit{\epsilon }<25\mathrm{\%})$, 30($\mathit{\epsilon }>25\mathrm{\%})$

2.5. Material model for the annulus

The annulus can be perceived as a structure similar to composites. It consists of a gel-like ground substance that binds the stack of fibers together which are oriented at ± 30° to the endplates (Bogduk, 2005). The fiber behavior to load is analogous to that of ligaments. Thus, the ground substance becomes responsible for load sharing during compression and the fibers take charge during tensile loads. Experimental studies on the annulus by researchers like Wagner & Lotz, (2004) and O’Connell et al., (2012) have presented its non-linear nature and also developed constitutive laws to define the behavior. In reality, though the annulus is viscoelastic, it is ignored in the current study as the physiological movements of the spine considered in this study like flexion/extension and lateral bending are static in nature. Earlier studies have proven that the hyperelastic material models can be used to represent the disc for such loads (Eberlein et al., 2004, 2001; Rohlmann et al., 2006)

The type-B model requires the treatment of the features of both the fibers and the ground substance as a continuum in an element. So far, there exists only one study (Guan et al., 2006) that proposes the idea of using an isotropic hyperelastic material model of type-B. The idea presented by Guan et al. was novel and shown to be effective in their assessment of range of motion (ROM) in sagittal and lateral movements for the in-silico model for L4-S1 multi-segment. This method was not found to be used in any other studies thereafter and hence received the interest to be used in the current study. The type-B anisotropic material model was essentially the GOH model and its manifestations (Eberlein et al., 2001; Gasser et al., 2006) that have been incorporated in some of the in-silico studies (Momeni Shahraki et al., 2015; Moramarco et al., 2010). Moramarco et al. (2010) implemented the material model using the anisotropic hyperelastic constants characterized by Eberlein et al. (2001) on L1-L2 intervertebral disc, in their study of the L1-S1 multi-segment in-silico model. In their study, particularly the L5-S1 segment of the in-silico model was observed to produce the highest deviation from the experimental response among all the other segments. Though the factor of material model usage may not be completely responsible for this result, the possibility of its contribution to this anomaly cannot be ignored. The study by Eberlein et al. (2001) used two separate uniaxial tensile test data pertaining to circumferential and longitudinal samples to obtain the constants for their material model. Whereas in-vitro loading conditions, the disc would be subjected to bi-axial load. The effectiveness of using the biaxial data for obtaining the constants for the material model was also proven in the study by Momeni Shahraki et al., in their study of the L4-L5 segment.

In this study, the simplified approach of using the isotropic hyperelastic model proposed by (Guan et al., 2006) was investigated and compared with the anisotropic hyperelastic model of Gasser-Ogden-Holzapfel (Eberlein et al., 2001; Gasser et al., 2006) to assess their predictive capabilities. From here on the in-silico segment model with implementation of isotropic and anisotropic hyperelastic material model for the disc will be referred to as SM_iso and SM_aniso respectively.

2.6. Isotropic hyperelastic material model

Neither the description of the isotropic hyperelastic material model nor the reason for its choice was presented in the study by Guan et al. (2006). Thus, the implementation of this approach required the evaluation of the best possible material model to represent the annulus behavior. Three of the available material models in ABAQUS software which are Mooney Rivlin (Mooney, 1940; Rivlin, 1948), Yeoh (Yeoh, 1993), and Ogden (Ogden, 1972) were chosen for assessment. The strain energy density function (SEDF) corresponding to the above three material models is described as follows: -

(1) $W_{\mathit{Yeoh}}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+C_2{\left(I_1\left(\mathit{C}\right)-3\right)}^2+C_3{\left(I_1\left(\mathit{C}\right)-3\right)}^3+\frac{H}{2}\left(J^2-1\right)$(1)

(2) $W_{\mathit{Mooney}-\mathit{Rivlin}}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+C_2\left(I_2\left(\mathit{C}\right)-3\right)+\frac{H}{2}{\left(\mathit{J}-1\right)}^2$(2)

(3) $W_{\mathit{Ogden}}=\frac{2{\mu }_1}{{{\alpha }_1}^2}\left({{\lambda }_1}^{{\alpha }_1}+{{\lambda }_2}^{{\alpha }_1}+{{\lambda }_3}^{{\alpha }_1}-3\right)+\frac{1}{D_1}{\left(\mathit{J}-1\right)}^2$(3)

Where, ${\mu }_1$,${\alpha }_1$ and C_i’s, are constants, H and D₁ are indeterminate Lagrange multipliers which are determined from the boundary conditions as a hydrostatic pressure, I₁(C) and I₂(C) are the first and second invariants of Right Cauchy Green Tensor C, J is the Jacobian, λ₁, λ₂ and λ₃ are the stretches in principle directions. As the hyperelastic material models considered in Eq. (1(1) $W_{\mathit{Yeoh}}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+C_2{\left(I_1\left(\mathit{C}\right)-3\right)}^2+C_3{\left(I_1\left(\mathit{C}\right)-3\right)}^3+\frac{H}{2}\left(J^2-1\right)$(1) –3(3) $W_{\mathit{Ogden}}=\frac{2{\mu }_1}{{{\alpha }_1}^2}\left({{\lambda }_1}^{{\alpha }_1}+{{\lambda }_2}^{{\alpha }_1}+{{\lambda }_3}^{{\alpha }_1}-3\right)+\frac{1}{D_1}{\left(\mathit{J}-1\right)}^2$(3) ) are isotropic, for a uniaxial loading condition the stretches in the directions other than the loading direction are equal. Also, as most of the hyperelastic materials can be considered as incompressible, it leads us to the following conditions—

(4) ${\lambda }_1{\lambda }_2{\lambda }_3=1$(4)

(5) ${\lambda }_2=\frac{1}{\sqrt{{\lambda }_1}}$(5)

The Cauchy stress, ${\sigma }_i=\frac{\mathrm{\partial W}}{\mathrm{\partial }{\lambda }_i}{\lambda }_i+\mathit{H}$(6)

The 1^st Piola Kirchoff stress, $P_i=$ $\frac{{\sigma }_i}{\mathit{\lambda }}$(7)

Solving equations (1)(1) $W_{\mathit{Yeoh}}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+C_2{\left(I_1\left(\mathit{C}\right)-3\right)}^2+C_3{\left(I_1\left(\mathit{C}\right)-3\right)}^3+\frac{H}{2}\left(J^2-1\right)$(1) , (2(2) $W_{\mathit{Mooney}-\mathit{Rivlin}}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+C_2\left(I_2\left(\mathit{C}\right)-3\right)+\frac{H}{2}{\left(\mathit{J}-1\right)}^2$(2) ) and (3(3) $W_{\mathit{Ogden}}=\frac{2{\mu }_1}{{{\alpha }_1}^2}\left({{\lambda }_1}^{{\alpha }_1}+{{\lambda }_2}^{{\alpha }_1}+{{\lambda }_3}^{{\alpha }_1}-3\right)+\frac{1}{D_1}{\left(\mathit{J}-1\right)}^2$(3) ) using equations (4)(4) ${\lambda }_1{\lambda }_2{\lambda }_3=1$(4) to (7), the 1^st Piola Kirchoff stress corresponding to different strain functions was evaluated to be -

(8) $P_{\mathit{i}\left(\mathit{Yeoh}\right)}=2\left({\lambda }_1^1-{\lambda }_1^{-2}\right)\left[C_1+2C_2\left(I_1-3\right)+3C_3{\left(I_1-3\right)}^2\right]$(8)

(9) $P_{\mathit{i}\left(\mathit{Mooney}-\mathit{Rivlin}\right)}=2\left[\left({\lambda }_1^1-{\lambda }_1^{-2}\right)C_1+\left(1-{\lambda }_1^{-3}\right)C_2\right]$(9)

(10) $P_{\mathit{i}\left(\mathit{Ogden}\right)}=2\frac{{\mu }_1}{{{\alpha }_1}^2}\left[{\lambda }_1^{{\alpha }_1}-{\lambda }_1^{-\frac{{\alpha }_1}{2}}\right]$(10)

The experimental data of uniaxial tests on the annulus specimen from (Wagner & Lotz, 2004) was extracted using Origin software. The method of least square curve fitting was used to identify the most suitable material model among the above-mentioned material models for representing the annulus material which could confirm with the experimental data. An in-house MATLAB code was developed and used for this purpose. Figure 2 shows the result of the analysis. The Mooney-Rivlin material model showed a high level of deviation from the required data. The curve corresponding to the Yeoh material model was close to the experimental data in most of the loading regime, but it showed a major deviation at the end of compressive load. It is clearly evident from the graph that the Ogden material model provides the best fit to the experimental data in both the loading regimes. The ability of Ogden material model to represent asymmetric tension-compression behavior is also supported in an article (Moerman et al., 2016) and thus, it was selected to represent the annulus material. The values of the constant that were found are listed in Table 1.

Figure 2. Comparison of three isotropic hyperelastic material models to fit to the experimental results.

2.7. Anisotropic hyperelastic material model

The possibility of using the Gasser-Ogden-Holzapfel (GOH) material model for representing the annulus has been explored by researchers in the past as well as in the current study. It had been developed with the intention of representing fiber distribution in a matrix medium for arterial walls which was further extended to other applications with structural similarity like the annulus of intervertebral disc.

The SEDF of GOH material model is formulated as (Gasser et al., 2006)

(11) $\mathit{W}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+\frac{k_1}{2k_2}\sum _{\mathit{i}=1}^n\left\{\exp \left(k_2{\overline{E_i}}^2\right)-1\right\}+\frac{1}{D}\left(\frac{J^2-1}{2}-ln\mathit{J}\right)$(11)

Where,

(12) ${\bar{E}}_{\mathit{i}}=\mathit{\kappa }\left(I_1\left(\mathit{C}\right)-3\right)+\left(1-3\mathit{\kappa }\right)\left(I_{4\mathit{i}}-1\right)$(12)

where $⟨\cdot ⟩$ are Macaulay brackets, $C_1,k_1,k_2,\mathit{\kappa }$ are material constants, D is an indeterminate pressure. $\mathit{\kappa }\in \left[0,1/3\right]$ characterizes the degree of anisotropy of the structure such that, κ = 0 represents a perfect alignment of fibers whereas κ = 1/3 signifies an isotropic distribution of fibers. As the fibers in the annulus are considered to have an organized alignment of ± 30°, κ was considered as “0” which reduces the expression for ${\bar{E}}_{\mathit{i}}$ to (I_4i-1). The anisotropic invariants $I_{4\mathit{i}}$ are defined as

(13) $I_{4\mathit{i}}=\left(a_{0\mathit{i}}\otimes a_{0\mathit{i}}\right):\mathit{C}$(13)

to describe the direction-specific mechanical properties of collagen fibers, where ${\mathit{a}}_{0\mathit{i}}$ are directional unit vectors representing the orientation of the i^th family of collagen fibers ($\mathit{i}=1,2,3,\dots ,\mathit{N}$) and (:) is a symbol of contraction. $I_{4\mathit{i}}$ is equal to the square of stretches along the direction of the vector $a_{0\mathit{i}}$. The third term in Eq. (11)(11) $\mathit{W}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+\frac{k_1}{2k_2}\sum _{\mathit{i}=1}^n\left\{\exp \left(k_2{\overline{E_i}}^2\right)-1\right\}+\frac{1}{D}\left(\frac{J^2-1}{2}-ln\mathit{J}\right)$(11) accounts for incompressibility.

Assuming two family of fibers ($\mathit{N}=2$) embedded symmetrically with an angle $\mathit{\gamma }$ between the circumferential direction and the mean orientation $a_{0\mathit{i}}$ of the fiber families. The fourth invariants are equal because of the symmetry, and can be expressed as

(14) $I_{4\left(1\right)}=I_{4\left(2\right)}=I_4={\lambda }_1^2{sin}^2\mathit{\gamma }+{\lambda }_2^2{cos}^2\mathit{\gamma }$(14)

Then, the SEDF of Eq. (11)(11) $\mathit{W}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+\frac{k_1}{2k_2}\sum _{\mathit{i}=1}^n\left\{\exp \left(k_2{\overline{E_i}}^2\right)-1\right\}+\frac{1}{D}\left(\frac{J^2-1}{2}-ln\mathit{J}\right)$(11) can be rewritten as

(15) $\mathit{W}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+\frac{k_1}{2k_2}\left\{\exp \left(k_2{\left({\mathit{I}}_{4\left(1\right)}-1\right)}^2\right)-1\right\}+\frac{k_1}{2k_2}\left\{\exp \left(k_2{\left({\mathit{I}}_{4\left(2\right)}-1\right)}^2\right)-1\right\}+\frac{1}{D}\left(\frac{J^2-1}{2}-ln\mathit{J}\right)$(15)

Since the fourth invariants are equal, the above expression can be simplified as

(16) $\mathit{W}=C_1\left(I_1\left(\mathit{C}\right)-3\right)+\frac{k_1}{k_2}\left\{\exp \left(k_2{\left(I_4-1\right)}^2\right)-1\right\}+\frac{1}{D}\left(\frac{J^2-1}{2}-ln\mathit{J}\right)$(16)

Where, $I_1\left(\mathit{C}\right)=\mathrm{tr}\left(\mathit{C}\right)={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2$

Therefore, the principal stress components can be obtained using equation 6 as

(17) ${\sigma }_1=2C_1{\lambda }_1^2+4k_1\exp \mathit{Q}\left[{\lambda }_1^2{sin}^2\mathit{\gamma }+{\lambda }_2^2{cos}^2\mathit{\gamma }-1\right]\left[{\lambda }_1^2{sin}^2\mathit{\gamma }\right]+\mathit{H}$(17)

(18) ${\sigma }_2=2C_1{\lambda }_2^2+4k_1\exp \mathit{Q}\left[{\lambda }_1^2{sin}^2\mathit{\gamma }+{\lambda }_2^2{cos}^2\mathit{\gamma }-1\right]\left[{\lambda }_2^2{cos}^2\mathit{\gamma }\right]+\mathit{H}$(18)

(19) ${\sigma }_3=2C_1{\lambda }_3^2+\mathit{H}$(19)

Considering biaxial loading condition for the reasons explained before. If σ₁ and σ₂ are the principal stresses in the loading plane, it yields to ${\mathit{\sigma }}_3=0$ which becomes the boundary condition to find the indeterminate pressure, H. These stresses are expressed as

(20) ${\sigma }_1=2C_1\left({\lambda }_1^2-{\lambda }_3^2\right)+4k_1\exp \mathit{Q}\left[{\lambda }_1^2{sin}^2\mathit{\gamma }+{\lambda }_2^2{cos}^2\mathit{\gamma }-1\right]\left[{\lambda }_1^2{sin}^2\mathit{\gamma }\right]$(20)

(21) ${\sigma }_2=2C_1\left({\lambda }_2^2-{\lambda }_3^2\right)+4k_1\exp \mathit{Q}\left[{\lambda }_1^2{sin}^2\mathit{\gamma }+{\lambda }_2^2{cos}^2\mathit{\gamma }-1\right]\left[{\lambda }_2^2{cos}^2\mathit{\gamma }\right]$(21)

Where, $\mathit{Q}=k_2{\left({\lambda }_1^2{sin}^2\mathit{\gamma }+{\lambda }_2^2{cos}^2\mathit{\gamma }-1\right)}^2$

An inconsistency was observed in the values adopted for the constants $C_1,k_1,k_2$ associated with the GOH material model in different studies (Eberlein et al., 2001; Moramarco et al., 2010; Shahraki, 2014) which demands for the need of evaluating these constants. The evaluation of constants was performed by an in-house code developed in MATLAB. The experimental data for biaxial stress state was obtained from O’Connell et al., (2012) using Origin software. The curves corresponding to 1:1 strain ratio (equibiaxial condition) were chosen for this purpose. The factor corresponding to shear stress was assumed to negligible. The R-square error for the curve-fit was computed as 0.81. The values of the constants for which maximum compatibility was achieved between the experimental data and the material model are enlisted in Table 1. Figure 3 shows the orientation of fibres as vectors in the annulus.

Figure 3. Fibre orientation in the Annulus shown as vectors (a) θ = +30°, (b) θ = −30°.

2.8. Boundary conditions

The sacrum was fixed with ENCASTRE on its centroid, which was connected to the nodes of the sacrum by multi-point constraints. As the intention of the study is to evaluate the material models to display the non-linear behavior of the disc to a better degree of compliance, a low moment load of 4 Nm was used, since the disc is proven to show high non-linearity in its behavior for loads with small magnitude (White & Panjabi, 1990). Pure moment of 4 Nm was applied on the mid-reference point of the upper surface of L5 vertebrae which was kinematically constrained to the nodes of that surface. This moment load was aligned to mid-sagittal and mid-coronal planes to produce the desired movements in the respective planes.

3. Results and discussions

The current study aims to compare the parameter associated with the validation of the L5-S1 spinal unit model with the use of isotropic and anisotropic hyperelastic material model for the annulus of the intervertebral disc. The range of motion for sagittal and lateral movements was used for this purpose from published data (Guan et al., 2007) as it is found to be the only available experimental data for the range of load considered in this study. Apart from it, the observations related to the contact pressure at the facets, the maximum von Mises stress on the vertebra, the midplane Tresca shear stresses of the annulus and strain in the ligaments for the two material models in different motion scenarios are discussed later.

Almost all of the in-silico studies in the literature which involved the L5-S1 segment with few exceptions (Charriere et al., 2003; Ramakrishna et al., 2018b) carried out before, that involves the L5-S1 segment, have been observed to include transverse ligaments as one among the different other ligaments, which is an anomaly with respect to clinical facts. The L5-S1 segment does not have the transverse process, whereas it has the iliolumbar and lumbosacral ligaments which are not found elsewhere in the spine. Besides this, the material properties of these ligaments are unavailable in any literature. An attempt was made to consider these ligaments in the developed in-silico model with the properties of anterior longitudinal ligaments. It was observed that there was no appreciable change in the response obtained in the in-silico model through this method. Also, the experimental data considered for the comparison has been obtained after the excision of these ligaments in their study (Guan et al., 2007). Hence, the in-silico model developed for this study excludes both the iliolumbar and lumbosacral ligaments.

3.1. ROM curves

Figure 4 shows the response curves for sagittal and lateral movements of the L5-S1 segment model. The measure of relative error between the curves corresponding to the two models with the experimental response was used to assess the efficiency of the model. In flexion, the response of SM_iso shows a slightly better performance with a lower relative error of 15% compared to 26% for SM_aniso in flexion, but during extension the value of relative error reached undesirably high value of 60% for SM_iso compared to 33% for SM_aniso. The response of SM_aniso was stiffer compared to SM_iso throughout the entire range of load. Though there was a slight slip of the SM_aniso response from the required range in the beginning of extension moment, the overall response curve was inside the acceptable domain of values. The response from SM_aniso showed excellent agreement with the experimental response curves throughout the entire range of moment in lateral direction when compared with SM_iso. The relative error of SM_aniso for right and left lateral moments were 3% and 5% respectively when compared to the computed values of 33% and 36% for SM_iso. Unusually high deviation from experimental response curve was observed in the initial stages of application of moment in both the left and right lateral bending motion for SM_iso. This deviation reduces and the curve lies inside the acceptable range of values on moving towards the end of the load.

Figure 4. Range of motion plots for, (a) Flexion-Extension movement and (b) Left and Right Lateral movements.

3.2. Facet contact pressure

It is clearly evident from Figure 5 that except for the lateral moment to the left, there is significant difference in the facet contact pressure between the SM_iso and SM_aniso models for the remaining types of moments. 40% rise in contact pressure was observed for flexion in SM_iso in comparison to the contact pressure for SM_aniso. The contact pressure for SM_aniso in the extension moment was 18% higher than the contact pressure for SM_iso. Similar contact pressure was obtained for both the models in lateral moment to the left. No contact pressure was produced in the lateral bending movement to the right for the SM_aniso compared to 27 MPa of contact pressure observed for SM_iso.

Figure 5. Comparison of facet contact pressure (MPa) for all the considered movements between the two models (FLEX- Flexion, EXT- Extension, LLAT—Left lateral bending, RLAT—Right lateral bending).

3.3. Maximum von mises stress on the vertebrae

The maximum stress in flexion and extension moments according to the von-Mises criterion was predicted on the right pedicle of the L5 vertebra for both SM_iso and SM_aniso. The left pedicle was found to be the next location of high stress. The difference in values of stresses between the two models in flexion and extension were found to be 7% and 12% respectively. The location of maximum stress in lateral bending was observed to creep in the posterolateral region near to the endplate of the L5 vertebra which was in contact with the disc. Except for the case where SM_aniso was subjected to lateral bending to the right, the remaining cases with lateral bending moments experienced the peak stress towards the posterior and to the same direction as that of the moment. As opposed to this observation, in the case of SM_aniso loaded with lateral bending moment to the right, the peak stress was observed on the left posterior corner of the endplate of L5 vertebra. A significantly high difference was found in values of stresses between the two models in left and right lateral moments. The percentage change of values for the two models in the left lateral movement was 54% and that for the right lateral moment was 58%. Figure 6 summarizes the above data.

Figure 6. Comparison of maximum von-Mises stress (MPa) on the vertebrae for all the considered movements between the two material models (FLEX- Flexion, EXT- Extension, LLAT—Left lateral bending, RLAT—Right lateral bending).

3.4. Strain in the ligaments

The ligaments were modeled to be responsive to tensile load. It is evident in Figure 7, that the peak strains were experienced by capsular ligament in all the movements. In majority of the cases the ligaments in the SM_iso model showed higher strains compared to the ligaments in SM_aniso which is the result of higher flexibility of the SM_iso model than the SM_aniso model. The highest strain of 72% was observed for the capsular ligament in the left side of SM_iso model followed by 66 % strain in the capsular ligaments on the right side as well as in the interspinous ligaments of the same model during flexion moment. The rest of the ligaments had low levels of strain ranging from 10 to 15% and there was no contribution of anterior longitudinal ligaments in flexion, as expected. The same trend is also displayed for SM_aniso. For the extension moment, the strain was observed only in the capsular ligaments and the anterior longitudinal ligaments with a dominance in the amount of strain on capsular ligaments. The maximum strain in capsular ligament was observed as 62% on the right side of the SM_iso model whereas it was 10% in the anterior longitudinal ligaments.

Figure 7. Comparison of ligament strain for all the considered movements between the two models (FLEX- Flexion, EXT- Extension, LLAT—Left lateral bending, RLAT—Right lateral bending).

3.5. Midplane shear stresses of the annulus based on Tresca criterion

Shear stresses have been considered in this study as they are identified as the cause for initiation of annular disruptions (Costi et al., 2007). From Figure 8 it can be observed that the SM_aniso model indicated higher midplane shear stress in the annulus compared to the SM_iso in all the movements. The higher stress values that are observed to occur in SM_aniso are attributed to the cumulative effect of ground substance and fibres of the annulus incorporated in GOH model. The percentage of difference in values between the two models for flexion, extension, left lateral bending and right lateral bending were 146%, 27%, 36% and 160% respectively.

Figure 8. Comparison of maximum shear stresses in midplane of annulus for all the considered movements between the two models.

The extension motion was found to be of least risk as compared to all the other motions that were considered in the study. This observation was in agreement to the study by Moramarco et al., (2010) where shear stress concentration was used to assess the damage initiating location. The disc in the L5-S1 region was found to be the most susceptible to damage. Also, it was observed in their study that flexion had the most damaging effect among the different types of motion. In the current study the lateral movements showed the possibility to cause even more damage in the annulus compared to flexion. This may be attributed to various factors like the inclusion of transverse ligaments or the higher load of 10 Nm or the simplified geometry of the disc or even the material model considered in their study.

The location of peak stress in the two models had a similarity for all the cases of movement. It was observed where the disc was expected to experience tension. The high stress in flexion was found to be at the right posterolateral region of the annulus. During extension, the right anterolateral portion of the annulus experienced the high stress. In the case of left lateral bending, the peak stress occurred at the right lateral region and vice versa was observed for right lateral bending in both the models.

4. Conclusion

The in-silico model of L5-S1 segment that was developed in this study includes a more geometrically accurate model of the disc. It was observed that the ligaments which are exclusively found in this region (iliolumbar and lumbosacral ligaments) did not affect the ROM of the in-silico model in this study. Therefore, those ligaments were omitted from the model. The validation was carried out with a range of motion at a load range of low magnitude which is more critical. The differences in validation capabilities and the responses of other valuable parameters with the use of GOH anisotropic hyperelastic material model and Ogden isotropic hyperelastic material model are presented. The SM_aniso model proves to be a stiffer model compared to SM_iso as it is clearly visible in the ROM curves and the midplane Tresca shear stresses of the annulus. The SM_aniso model is found to remain in the corridor of acceptance with more consistency compared to the SM_iso model. The results indicate that it is preferable to use the GOH model for the annulus as it is in better agreement with the in-vitro results. This model can be an effective substitute for the in-vitro specimens to perform further studies related to the ailments and surgical interventions which may occur in the L5-S1 segment.

Consent for publication

All authors consent to the publication of this manuscript.

Ethics approval and consent to participate

An institutional ethical clearance has been obtained from Kasturba Medical College and Kasturba Hospital’s institutional ethical committee. [IEC Project No.IEC:931/2018].

Author contribution statement

Dr. Raviraja Adhikari is responsible for the ideation and guidance. Dr. Shyamasunder Bhat N helped to procure the source data and provided the necessary guidance. Dr. Subraya Krishna Bhat developed the MATLAB code to test the mathematical models. Mr. Vinyas developed the finite element model, performed simulations, analyzed the results and wrote the main manuscript. All the authors have reviewed the manuscript.

Acknowledgements

The author would like to express their heartfelt gratitude to Manipal Academy of Higher Education, Manipal for providing the required facilities to perform this study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

  Vinyas

Vinyas is working as Assistant Professor-Senior Scale in the Department of Mechanical and Industrial Engineering, Manipal Institute of Technology, Manipal Academy of Higher Education, Manipal, Karnataka, India. He holds B.E. (Mechanical Engineering) and M.Tech (Computer Aided Mechanical Design and Analysis) degrees. He has more than eight years of teaching and research experience. His area of interest includes Biomechanics of the spine, Finite element analysis, Fatigue, and Continuum damage mechanics. The author has published articles related to Finite element study of spine and Continuum damage mechanics.