The effects of tether pretension within vertebral body tethering on the biomechanics of the spine: a Finite Element analysis

This study investigates the biomechanics of the spine after insertion of vertebral body tethering (VBT) with different cord pretensions. For that purpose, a Finite Element model of the average thoracolumbar spine was stepwise calibrated and validated. The VBT instrumentation was inserted in the left side of the L1-L2 segment with different cord pretensions. As a second test, the L1-L2 segment was submitted to an external pure moment of 6 Nm in left and right lateral bending. The range of motion (ROM) for the spine with VBT was determined with respect to its initial post VBT position. Pretension forces of 100 N and 300 N resulted in a change of scoliotic angle of 2.7° and 5.3° to the left side of the spine, respectively. The ROM of the native spine was 4.5° in right lateral bending and reduced to 1.8° and 1.4° for the cases of the spine with a cord pretension of 100 N and 300 N, respectively. In left lateral bending, the absolute ROM of the native spine was 4.6°. For the cases of a cord pretension of 100 N and 300 N, the spine bent 1.9° and 0.8° to the left side from its initial post VBT position, respectively.


INTRODUCTION
Adolescent idiopathic scoliosis (AIS) is a complex 3D deformity of the spine that worsens during skeletal growth.The spine develops a side-to-side curvature, usually in an elongated "S" or "C" shape in the frontal plane instead of growing straight (Šarčević and Tepavčević, 2019).Anterior vertebral body tethering (VBT) is a fusionless treatment option that is considered a new promising method for the management of selective cases of AIS (Pehlivanoglu et al., 2020).It is a growth modulating technique that derives benefit from the Heuter-Volkman principle which states that increased mechanical compression load in comparison with normal values slows growth and distraction accelerates it (Stokes et al., 1996).During surgery, this effect is achieved by tensioning the polyethylene-terephthalate tether (cord) that is connected to screws inserted in the lateral convex side of the vertebral bodies.One of the main advantages of VBT over spinal fusion is the ability to preserve the mobility of the spine while maintaining growth potential.
Multiple clinical reports have shown that VBT is safe and effective in correcting deformity of the spine (Crawford and Lenke, 2010;Miyanji et al., 2020;Samdani et al., 2015).Recently, Miyanji et al. (2020) reported the cases of 57 consecutive patients who underwent VBT mostly at the thoracic region (96.5%).The mean preoperative major curve of 51° was significantly improved to a mean of 24.6° at the first postoperative visit, and 16.3° at one year (Miyanji et al., 2020).Biomechanical studies on human spines confirmed that VBT can preserve motion of the spine in flexion-extension and axial rotation while reducing motion in lateral bending (Lavelle et al., 2016;Nicolini et al., 2021b).
Despite the promising results, overcorrection and under-correction after VBT are common causes of revision surgery (Buyuk et al., 2021).The intraoperative tensioning of the cord is one key component to achieve appropriate curve correction but its effects are poorly investigated.Lavelle et al. (2016) evaluated the effects of pretension of the tether at the most superior and most inferior vertebrae as opposed to all vertebrae.The thoracic spines were tested by applying only a pretension of 100 N on the cord as also performed in the biomechanical study of Nicolini et al. (2021b) with thoracolumbar spines.However, the effects of further tensioning were not investigated although a tensioning of 200 N and 300 N are common in clinical scenarios.Therefore, it remains unknown how much tension to place on the scoliotic spinal segments within the VBT construct although this is a key engineering parameter to control correction of the deformity.
Numerical data that could contribute to the clinical and biomechanical studies as well as to the improvement of the VBT technique remains sparse.Finite Elements (FE) models of the spine have made important contributions in the field of biomechanics but accurate and clinically relevant modeling of the human spine remains challenging (Dreischarf et al., 2014).Few FE models (Ezquerro et al., 2011;Heuer et al., 2007b;Schmidt et al., 2007) are calibrated considering experimental data from stepwise resection studies (Li et al., 2017;Zhang et al., 2011).The experimental procedure is conducted by testing the spine in different loading directions and gradually resecting the passive spinal structures to obtain their individual contribution to the biomechanics of the spine.Posteriorly, each spinal structure is sequentially added to the FE models and calibrated by varying and choosing the material properties that best fit the experimental results (Ezquerro et al., 2011;Nicolini et al., 2021a;Schmidt et al., 2006).This method allows for more precise modeling than obtaining material properties from the literature data and accounts for preferably physiological behavior of the spine (Damm et al., 2020;Schmidt et al., 2007).Moreover, most FE models do not take into consideration the regional variance of geometrical and material properties of the intervertebral disc (IVD) (Cassidy et al., 1989;Guerin and Elliott, 2006;Holzapfel et al., 2005) although they are important characteristics for the accuracy when predicting internal stresses such as the intradiscal pressure.
Another important step in the development of FE models of the spine is the construction of the geometry since it strongly influences the response of numerical results (Niemeyer et al., 2012).One possible approach to this problem is to choose a geometry obtained via image segmentation of an individual (e.g.computed tomography) and calibrate its material properties from the aforementioned experimental data.The geometry is derived from a real spine and consequently has feasible anthropometric proportions.However, it may carry deviations from an anthropometric mean that may induce numerical results with biases.This may result in difficulties in obtaining a good agreement with the average/median values regarding the biomechanics of the spine reported in the literature.An alternative approach is to construct or controllably modify a geometry to statistically represent the average values of anthropometric data.In this case, the model may carry some geometric relationship that is not verified in a human being since it was "synthesized".On the other hand, with the statistically representative construction of the anthropometric data, accompanied by an observation of the anatomical proportions and proper assembly, it is possible to minimize the risk of incurring the abovementioned problem.Furthermore, this allows for estimating and reporting average values, which is the preferred method of data publication for most clinical and biomechanical studies.
The main objective of this study is to investigate the biomechanical behavior (change of scoliotic angle and range of motion) of the spine after insertion of VBT with different cord pretensions.For that purpose, a FE model statistically representative of the average spine was developed, stepwise calibrated, and validated.The study is of clinical interest Latin American Journal of Solids and Structures, 2022, 19(3), e442 3/17 since the cord pretension is an important factor to avoid overcorrection and under-correction of scoliotic spines treated with VBT which are common causes of revision surgery.It also fills a gap in the literature since there is no study that numerically evaluated the effects of different cord tension within VBT on thoracolumbar spines.

MATERIALS AND METHODS
The FE model of the thoracolumbar spine was developed in Abaqus® (Dassault Systèmes, Waltham, MA, USA) as shown in Figure 1.

Geometry
The geometry of the vertebrae and sacrum (T10-S1) was derived from CT images of a man who had no history of spinal pathology.It was modified to represent an average spine of a young man following a manual iterative procedure.Firstly, the original geometry of each vertebra and sacrum was modified by translating and rotating the control points that define the shape of the surfaces (Non-Uniform Rational B-Spline -NURBS) using the software Rhinoceros® (Robert McNeel & Associates, Seattle, USA).The dimensions considered during this process (Figure 2) were the superior and inferior vertebral body width (SVBW and IVBW), superior and inferior vertebral body length (SVBL and IVBL), anterior and posterior vertebral body height (AVBH and PVBH), lateral vertebral body height (LVBH), transverse process length (TPL), vertebral canal superior width and length (VCSW and VCSL), superior and inferior interfacet width (SFFW and IFFW), interfacet height (FFH), superior and inferior facet length (SFL and IFL), superior and inferior facet width (SFW and IFW), superior and inferior longitudinal facet angle (SLFA and ILFA), and superior and inferior transverse facet angle (STFA and ITFA) as well as the morphology of the vertebral endplates (Masharawi et al., 2010(Masharawi et al., , 2005(Masharawi et al., , 2004;;Masharawi and Salame, 2011;van der Houwen et al., 2010).Secondly, the vertebrae were positioned in the 3D space considering the intervertebral disc (IVD) height and lordosis angle (Been et al., 2011(Been et al., , 2010;;Boulay et al., 2006;Cheng et al., 1998;Eijkelkamp et al., 2007;Labelle et al., 2005;Nissan and Issachar, 1986;Vialle et al., 2005;White and Panjabi, 1990).The geometries of the IVDs were created by sketching their lateral profile and using the morphology and positions of the vertebrae.The nucleus pulposus was dimensioned with 44% of the total cross-section area of the IVD (White and Panjabi, 1990) with its center located approximately at 3.5 mm posterior from the center of the IVD (Noailly et al., 2007;Weisse et al., 2012).
The following equation was used to compare the dimension of the model (  ) with the anthropometric data (  ) with respect to the anthropometric standard deviation (): (1) Figure 2 Dimensions considered during the modification of the geometry of the vertebrae.

Mesh
Tetrahedral elements were used to represent the vertebrae, as they can conform to complex geometries (Table 1).The cortical bone of the vertebral body was modeled using shell elements with a wall thickness of 0.29 mm (Ritzel et al., 1997).The IVD mesh was created using hexahedral elements.Only 1% of the FEs representing the vertebrae and IVD had an aspect ratio greater than 3, which is acceptable according to Burkhart et al. (2013).Tie constraint was used to guarantee continuity between the superior and inferior surfaces of the IVDs and their adjacent vertebrae since their meshes are different.
A mesh quality assessment was performed to evaluate the optimal mesh density that provides accurate results while keeping simulation time tractable.For that purpose, the L1-L2 IVD was tested using the traditional pure moment load protocol with a magnitude moment of 7.5 Nm in flexion.The range of motion (ROM, angular displacement of L1 vertebra with respect to L2 vertebra) and pressure at the center of the IVD were evaluated for different mesh sizes using quadratic shapes functions.A mesh sensitive analysis was not performed with the vertebrae since they have a minor effect on the movements of the spine compared to the soft tissue.Additionally, linear interpolation was used for the FEs representing the vertebrae.

Constitutive models
The nucleus pulposus was modeled as a hyperelastic Mooney-Rivlin material using the strain energy density function Latin American Journal of Solids and Structures, 2022, 19(3), e442  5/17 where  10 ,  01 and   are the material constants defining matrix stiffness and compressibility.The elastic volume ratio is represented by   while  and  are the first and second invariant of the right Cauchy-Green deformation tensor, respectively.The subscript (. )  indicates the belonging of the material parameters to the nucleus pulposus.
The annulus fibrosus was modelled using the hyperelastic Holzapfel-Gasser-Ogden model (Holzapfel et al., 2000).The total strain energy density function is composed of the strain energy density function of the isotropic matrix represented by the Neo-Hookean model and the strain energy density function of the fibers.It has the form with where  10 and  are material constants of the annulus ground substance related to its stiffness and compressibility.The material parameters  1 and  2 describes the exponential stress-strain relation for the collagen fibers, and  represents the level of dispersion in the fiber directions.The constitutive model considers that the fibers provide resistance only to traction.Additionally,  4() are the pseudo-invariants of the distortion part of the right Cauchy-Green deformation tensor and unit vectors in the fiber directions.The number of fiber families  was considered equal to two.Moreover, the orientation of the fibers at the anterior aspect of the IVD with respect to its transverse plane is defined by the parameter .
The annulus fibrosus was divided into five radial and five circumferential sections (Figure 3) to simulate the variance of properties in the annulus found in histological anatomy studies and experimental studies (Ebara et al., 1996;Eberlein et al., 2001;Holzapfel et al., 2005;Zhu et al., 2008).The coefficients  1 and  2 defining the material properties of the fibres were varied linearly in the circumferential and radial direction.The coefficients  1 ,  2 and   were defined to be responsible for creating a gradient of  1 ,  2 and  values along the perimeter of the annulus fibrosus, respectively.Similarly,  1 and  2 were set to control the change of  1 and  2 values along the radial direction of the annulus fibrosus, respectively.For these reasons, the subscripts (. )  and (. )  indicate the parameter is responsible for creating a gradient of the material properties values in the circumferential and radial direction of the annulus fibrosus, respectively.Operationally, the properties of each layer of the IVD are calculated based on the values of the outer layer of the most anterior aspect of the IVD (region A1, Figure 3).For instance, a  1 equal to 0.1 means that the  1 value of region B1 and  1 value of region C1 are 10% and 20% higher than the  1 value of region A1, respectively.The spinal ligaments were modeled as bar elements of nonlinear stiffness connecting the adjacent vertebrae (Figure 4).The gap at each pair of facet joints was modeled using soft contact interactions of the type nonlinear penalization.The model does not consider any tangential resistance at the surface of the facet joints since their coefficient of friction is negligible (Cramer and Darby, 2014).The input data was calculated using the equation where  is the contact pressure and ℎ is the overclosure while ,  and ℎ 0 are material constants.

Material properties and calibration
The cortical and cancellous bones of the vertebral body were modeled as an isotropic elastic material with a Young Modulus of 12,000 MPa and 100 MPa, respectively (Guan et al., 2006;Kulduk et al., 2015;C.-L. Liu et al., 2011).The posterior elements of the vertebrae were simplified by considering only one type of bone with intermediate stiffness considering a Young modulus of 3,500 MPa (Ayturka and Puttlitz, 2011;Guan et al., 2006;Shirazi-Adl et al., 1986).
The mechanical properties of the soft tissues and the contact of the facet joints were calibrated considering the inverse order of experimental resection studies (Coombs et al., 2013;Damm et al., 2020;Heuer et al., 2007a;Jaramillo et al., 2016).The spinal structures were sequentially added to the FE models and calibrated by varying and choosing the set of material properties that best fit the experimental results (Ezquerro et al., 2011;Nicolini et al., 2021a;Schmidt et al., 2006).
The experimental tests were performed as described in our previous studies (Beckmann et al., 2020(Beckmann et al., , 2018;;Nicolini et al., 2021b).The spine specimens had the vertebra, intervertebral discs, and the seven aforementioned ligaments preserved.They were tested under a pure moment load and a rotational speed of 1°/s in flexion-extension, lateral bending and axial rotation with the custom-built machine (Nicolini et al., 2021b).The load was recorded by torque sensors that provided a maximum error of 1% relative to the target value while the kinematics of the vertebrae were measured by a electromagnetic tracking system that provides deviations lower than 0.1 mm and 0.1° (Beckmann et al., 2018).The data was recorded at a frequency of 100 Hz (Brandes et al., 2021).Some of the L1-L2 spine specimens were tested in different resection states to obtain the individual contribution of the IVD, the facet joints and the ligaments to the kinematics of the spine.
For the calibration procedure, the experimental data was prepared in such a way that the FE model could: 1) replicate the median experimental moment-ROM curves of twelve L1-L2 intact segments in the intact state (some already published by Nicolini et al. (2021b), and 2) reproduce the proportional change of ROM after removal of a spinal structure as observed in our resection study with five L1-L2 spines in the resected state.For this purpose, the experimental median ROM of twelve L1-L2 intact segments (  ) were scaled using the equation: Latin American Journal of Solids and Structures, 2022, 19(3), e442 7/17 where   represents the predicted ROM of the spine with defects while  is a coefficient defined as The scale factor  represents a proportional increase in ROM caused by the resection of a spinal structure.It was calculated by dividing the ROM obtained after the resection of a spinal structure (  ) by the ROM of the native spine (  ).This process was performed for each of the five spine segments and the median was taken.Moreover, this rule was applied for each loading direction and datapoint of the moment-ROM curve individually (Nicolini et al., 2021a).
The R-squared function comparing the numerical values with the experimental ones was implemented as a cost function to be minimized.The calibration of the material properties of the IVD were calibrated considering motion within the three anatomical planes while the material properties of the contact of the facet joints were calibrated considering axial rotation.The following design vectors defining the mechanical properties of the IVD and the facet joints were set to be optimized via an in-house developed Genetic Algorithm.
An optimization algorithm was used to unequivocally determine the force-displacement curves that represent the material properties of the ligaments and replicate the experimental moment-ROM data of the spine in flexion-extension.During all simulations, the degrees of freedom of the lower part of the L2 vertebra were constraint and the L1 vertebra was bent similar to the experimental flexibility tests.For the sake of shortness, this article presents only the calibration results of the L1-L2 segment.

Verification and validation
After the calibration of the material properties of the L1-L2 IVD under a pure moment, the tests of Holzapfel et al. (2005) were numerically replicated to verify the regional variation of the tensile properties within the IVD.They performed tensile tests with single lamellar annulus specimens whereby the lamellar fiber direction was aligned with the load axis.The specimens were obtained from the ventro-lateral external (VLe), the ventrolateral internal (VLi), the dorsal external (De), and the dorsal internal (Di) of the L1-L2 IVD.Moreover, the forces at the contact of the facet joints were verified with the experimental values of Niosi et al. (2008) by testing the L1-L2 segment under a pure moment of 7.5 Nm in different loading directions.
For validation purposes, the L1-L2 segment was submitted under multiple bending moments in flexion-extension or lateral bending combined with axial rotation.The numerical ROM values were compared to those obtained experimental by Wilmanns et al. (2021).

Modeling and simulation of the spine with VBT
The geometry of a VBT system was obtained from a manufacturer (Globus Medical Inc, Pennsylvania, USA) and tested using the model of L1-L2 segment (Figure 5).The screw was designed with an inner diameter of 4.75 mm and an outer diameter of 5.5 mm.Oriented by an experienced surgeon, the screws were inserted close to the foramen with bicortical purchase at the left lateral side of the L1 and L2 vertebrae and parallel to the frontal plane.The anchors had an inner and outer diameter of 5.5 and 6.5 mm, respectively.The screws and anchors were modeled as an isotropic material made of titanium with a Young's Modulus of 110 GPa and a Poisson's ratio of 0.3 (Lo et al., 2011;Niinomi et al., 2016).The mesh was created using tetrahedral elements with linear shape functions.Tie constraint was used to attach each anchor and screw to the vertebra.The cord of 4 mm diameter was modeled using a bar element with a Young's Modulus 1,500 MPa in tension and negligible resistance to compression (Kulduk et al., 2015;C. L. Liu et al., 2011).The first test evaluated the changes in the coronal plane angles (scoliotic angles) due to a cord tension of up to 300 N within VBT.In the second test, the spine was simulated under an external pure moment up to 6 Nm in left and right lateral bending (Lavelle et al., 2016;Nicolini et al., 2021b) considering an initial cord pretension of 0 N, 50N, 100 N, 200 N and 300 N. Pretension was applied by controlling the initial stretch of the bar element representing the cord.The ROM was determined with respect to the spine at its respective initial post VBT position (no external load independent of cord pretension).

Geometry
The final geometric model of the vertebrae presented an average of the absolute  values of 18% and a maximum absolute Z value of 74% considering the dimensions SVBW, IVBW, SVBL, IVBL, AVBH, PVBH, LVBH, TPL, VCSW, VCSL, SFFW, IFFW, FFH, SFL, IFL, SFW, IFW, SLFA, ILFA, STFA and ITFA (Table S1 of supplementary material).The morphology of the vertebral endplates was within 10% of the SD of the anthropometric data.The average of the absolute Z values was 29% for the wedging angle of the vertebrae and IVDs, 53% for sagittal dimensions including lumbar lordosis, sacral slope, pelvic tilt, pelvic incidence and thoracic kyphosis, and 20% for the anterior and posterior height of the IVDs.The comparison of the dimensions of the geometric model and the experimental data are presented in Table S1, Table S2, Table S3, and Table S4 of supplementary material.The lordosis angle of the model represents a young male with a type 3 lordosis (Roussouly et al., 2005).

Mesh
The mesh-sensitive analysis showed an asymptotic behavior of the ROM and intradiscal pressure with the refinement of the mesh (Figure 6).An intermediate mesh of 27991 nodes with quadratic shape functions was selected to model the IVD which had a difference of 0.6% of ROM and 4.5% for the intradiscal pressure when compared to the results of the finest tested mesh (68592 nodes).
Figure 6 Effect of mesh refinement on numerical ROM and intradiscal pressure at the center of the disc.

Calibration and validation
After calibration, the model of L1-L2 IVD was capable of reproducing the experimental moment-ROM curves presenting an average R-squared of 0.9 considering flexion, extension, lateral bending, and axial rotation (Figure 7).The material parameters are shown in Table S5 of supplementary material.The numerical stress-stretch ratio curves of single lamellae obtained from four anatomical regions of the annulus also presented a similar trend compared to the experimental results (Figure 8).All numerical data points were within SD of the experimental data and the average Rsquared was 0.95.The tensile moduli for external lamellae were approximately two times higher than those for internal lamellae.The tensile moduli of the ventral layers were only 10% higher compared to the dorsal layers.Both results agree with the values found by Holzapfel et al. (Holzapfel et al., 2005).
The fiber orientation angle with respect to the transverse plane gradually increases from 30° at the anterior aspect to 42° at the posterior aspect of the modeled IVD.This is consistent with the histological values for the anterior and posterior part of the IVD obtained by Holzapfel et al. (2005) (25.7 ± 5.8° and 49.3 ± 8.0°) and Cassidy et al. (1989) (28° and 45°).The total criss-cross angle of the fibers at the outer layer of the model was 60°, which is within standard deviation compared to the experimental results of 51.9 ± 8.4° obtained by Guerin et al. (Guerin and Elliott, 2006).The model of the L1-L2 segment agreed with the experimental results for different resection stages in flexion, extension, lateral bending, and axial rotation (Figure 9).The average R-squared for all loading directions and resection stages was 0.85.For the native spine, the mentioned variable was equal to 0.83.The model predicted a force at the facet joints of 0, 5, 4, and 79 N whereas Niosi et al. (2008) measured values of 2 ± 5, 13 ± 14, 11 ± 11, and 56 ± 17 N for flexion, extension, lateral bending, and axial rotation, respectively.The comparison of the numerical and experimental data of the spine under combined loading presented an R-squared of 0.98, 0.94, and 0.90 for flexion, extension, and lateral bending combined with axial rotation, respectively.

Effects of VBT
The scoliotic angle of the spine gradually increased with increasing of the cord pretension (Figure 10).Pretension forces of 50 N, 100 N, 200 N, and 300 N resulted in a change of scoliotic angle of 1.5°, 2.7°, 4.3°, and 5.3° to the left side of the spine, respectively.
Herein, the ROM is presented as the angular displacement of the L1-L2 segment from its neutral position or initial post VBT position to the final position caused by an external moment of ±6 Nm.The ROM of the native spine was 4.5° in right lateral bending and reduced to 1.9°, 1.9°, 1.8°, 1.6°, and 1.4° for the cases of the spine with a cord pretension of 0 N, 50 N, 100 N, 200 N and 300 N, respectively (Figure 11).In left lateral bending, the absolute ROM of the native spine and the instrumented spine with no pretension was 4.6°.For the cases of a cord pretension of 50 N, 100 N, 200 N, and 300 N, the spine bent 3.1°, 1.9°, 1.0°, and 0.8° to the left side from its initial post VBT position, respectively.For the case of 50 N and 100 N, the moment-ROM curve in left lateral bending presented an inflection close to the point in which the cord tension reduced to zero due to the applied external moment.

DISCUSSION
VBT is an interesting technique that has been approved in 2019 by the Food and Drug Administration (FDA) for US patients with AIS.As an emerging surgical technique, there is a paucity of studies on VBT especially regarding the effects of its engineering parameters that can substantially affect the biomechanics of the spine.This paper presented the first Latin American Journal of Solids and Structures, 2022, 19(3), e442 12/17 numerical study that evaluated the change of scoliotic angle and range of motion of the spine after insertion of VBT with different cord pretensions.It was motivated by the common clinical scenario related to overcorrection and under-correction of the spine treated with VBT that leads to revision surgery and the fact that the cord pretension is a key parameter to achieve appropriate scoliotic curve correction.For that purpose, this paper also presented the development of a FE model of the spine that is statistically representative of the average spine, stepwise calibrated, and validated.The proposed method successfully generated a geometric model with dimensions close to the anthropometric data and within SDs.The iterative process is justified by the difficulties of guaranteeing geometric values close to the experimental data, smooth surfaces, and a reasonable anatomical assembly while dealing with data from different studies and dimensions that are connected.This includes the alignment of the vertebral bodies, the vertebral canals, and the facet joints of adjacent vertebrae.Based on the obtained results, we consider that the model is statistically representative of the average of the human thoracolumbar spine.
The mesh-sensitive analysis showed an asymptotic behavior of the ROM and intradiscal pressure with the refinement.Based on these results, the selected mesh of approximately 28000 nodes with quadratic shape functions was considered dense enough for practical purposes and could capture important characteristics such as the morphology of the vertebral endplates.It also keeps the simulations numerically tractable since the clock time for the simulations of one spinal functional unit under pure moment within the anatomical planes took less than 5 min using a computer with an i7 processor and 8 GB of RAM.
The results of the developed FE model agreed with experimental results regarding the 1) ROM of the spine under pure moment in different loading directions, 2) stress-stretch ratio curves from tensile tests of single lamellae of four anatomical regions of the annulus, 3) angle of the fibers for different positions within the IVD, 4) forces at facet joints, 5) ROM of the spine under combined loading.We attribute this achievement due to several aspects.First, the geometric model was generated considering the average anthropometric data.This probably reduced the divergence between the numerical values and the experimental median/average values regarding the biomechanics of the spine presented in the literature.Second, a proper mesh size was chosen based on a mesh sensitive analysis.Third, each spinal structure of the FE model was calibrated considering their individual contribution to the motion of the spine using data from cadaveric stepwise resection studies.This allowed to properly dose the stiffness within the spinal structures.For instance, a satisfactory agreement of the numerical forces at facet joints was possible because the soft tissue was calibrated.A too soft IVD, for example, would increase the stress at the facet joints in axial rotation leading to a divergence of the facet joints forces.Furthermore, a powerful Genetic Algorithm was used to tune 12 parameters representing the material properties of the IVDs including characteristics that are not present in most published FE models such as the regional variance of fibers angles and material properties of the IVD.The algorithm is also able to deal with challenging aspects that can be present during the determination of the material properties of the IVD such as multiple local optima and non-smooth objective functions (Nicolini et al., 2021a).The Genetic Algorithm was also used to determine the parameters that define the contact of the facet joints.Moreover, the force-displacement curves representing the material properties of the ligaments were unequivocally determined.By varying the material properties of the soft tissue and the contact of the facet joints and choosing the set of material properties that produce the best fit between the numerical and experimental results gave a better agreement with the experimental data rather than directly applying the material properties from the literature (Ezquerro et al., 2011;Nicolini et al., 2021a;Schmidt et al., 2006).
The numerical results confirmed that cord pretension is a key factor to change the scoliotic angle of the spine.A pretension force of 100 N and 300 N resulted in a change of scoliotic angle of 2.7° and 5.3° to the left side of the spine, respectively.This value is associated with the flexibility/stiffness behavior of the spine and may be predicted based on the moment-ROM curve of the native spine under pure moment as an approximation.The cord was attached 26.2 mm from the middle sagittal plane of the spine.Therefore, a pretension of 100 N produces a moment of 2.6 Nm which also provides an angular displacement of 2.7° according to the moment-ROM of the native spine (Figure 11).In the cadaveric study of Lavelle et al. [24], a value of 9.9±5.5°was found for the case of a T4-T12 continuous sequentially tensioned with 100 N, which means 1.24° per spinal segment when assuming a uniform distribution along the spine.The seventeen patients studied by Newton et al. (2018) presented an average thoracic curve magnitude of 52 ± 10° preoperatively and 31 ± 10° immediately postoperatively.Considering that there was an average of 5.8 segments tethered per patient and assuming uniform distribution, the change of scoliotic angle was 3.6° per segment.However, Newton et al. (2018) did not report the value of the tension applied to the cord during surgery.
Primary bending stability in the opposite direction to insertion of VBT is an important requirement for VBT (Nicolini et al., 2021b).In right lateral bending, the stiffness generated from cord tension resulted in a reduction of ROM of the spine.The instrumentation with no cord pretension and with a pretension of 300 N were capable of reducing 59% and 69% of the ROM of the native spine, respectively.The relatively small change between the case of no pretension and the case with a cord pretension of 300 N suggests that the cord pretension is capable of changing the neutral position of the spine but does not add much stiffness to the spine in right lateral bending.The moment-ROM curve of the native spine in left lateral bending was not affected after inserting a cord without pretension since it was modeled with no resistance to compression.However, the implant inserted with pretension added stiffness to the spinal system.The ROM was reduced in left lateral bending mainly because the soft tissue with nonlinear stiffness was tensioned at the right side of the spine after insertion of VBT.
This study also presents some limitations.The effects of VBT were investigated using a FE model representing the average of adult spines, while VBT is usually used in adolescents with scoliosis.This limitation could potentially be overcome by scaling the geometric model and creating scoliosis.However, it is not fully understood the biomechanical effects after such modification and it is not guaranteed that it would provide a geometry that represents an adolescent.Moreover, there are no stepwise resection studies on adolescent scoliotic spines which could provide experimental data for the calibration of the FE model.We opted to use the model of a non-scoliotic spine to allow comparison with the biomechanical experiments with cadaveric spines.Furthermore, the bone parts were modelled as isotropic material (Guan et al., 2006;Kulduk et al., 2015;C.-L. Liu et al., 2011), whereas they are supposed to be anisotropic materials.This simplification would have a minor effect on the results of our study since most of the motion of the spine is defined by the properties of the soft tissue due to the differences in stiffness compared to the bone.However, a more sophisticated modelling might be necessary if one intends to perform an analysis involving the stresses acting in the vertebrae.

CONCLUSION
This numerical study evaluated the change of scoliotic angle and range of motion of the spine after insertion of VBT with different cord pretensions.We conclude that the cord pretension substantially influences the resulting scoliotic angle while the stiffness of the tether provides primary stability to the spine for lateral movement in the opposite direction of the implant.The increase of cord pretension leads to a reduction of ROM of the spine in lateral bending.Due to the very satisfactory agreement between the numerical and experimental results (geometry, ROM, stresses in the IVD, facet joint contact forces, etc.), we conclude that the model is calibrated, validated and has the potential be used as a tool to assist on design, testing and evaluation of spinal implants.

Figure 1
Figure 1 Perspective and lateral view of the intact T10-L3 spine.

Figure 3
Figure 3Intervertebral disc divided into five radial and five circumferential sections to account for the regional difference in material properties.The figure also indicates the subregion of C4.

Figure 5
Figure 5 Perspective view of the L1-L2 segment instrumented with vertebral body tethering.

Figure 7
Figure 7 Comparison of the numerical and experimental moment-ROM curves of the L1-L2 under pure moment in different loading directions.For the first subplot, positive and negative values represent flexion and extension, respectively.The average R-squared was 0.9 considering all loading directions.

Figure 8
Figure 8 Comparison of the numerical and experimental engineering stress-stretch ratio curves of single lamellae obtained from four anatomical regions of the annulus fibrosus: the ventro-lateral external (VLe), the ventro-lateral internal (VLi), the dorsal external (De), and the dorsal internal (Di).The average R-squared for all loading directions was 0.95.

Figure 9
Figure 9 Comparison of the numerical and experimental results of the L1-L2 segment under pure moment in different loading directions and resection states.The error bars indicate the difference between the numerical and experimental results.

Figure 10
Figure 10 Change of scoliotic angle as a function of the tension applied on the cord within anterior vertebral body tethering.A tensioning of 300 N resulted in a change of scoliotic angle of 5.3° as shown in the subfigure on the right side.

Figure 11
Figure 11 Moment-ROM curves for the spine with different cord pretension in left and right lateral bending.Positive and negative values represent right and left lateral bending, respectively.
The effects of tether pretension within vertebral body tethering on the biomechanics of the spine: a Finite Element analysis Luis Fernando Nicolini et al.Latin American Journal ofSolids and Structures, 2022, 19(3), e44213/17

Table 1
Element types used for each structure in the FE model.