Phantom-Based Lumbar Spine Experimental Investigation and Validation of a Multibody Model

At present, experimental and numerical investigations are being exploited hand-in-hand to answer current biomechanical issues effectively. Therefore, the implementation of a joint approach has permitted relevant steps forward in the understanding of various biomechanical aspects, thus permitting relevant improvements in the clinical sphere, such as in surgical or biomedical procedures [1–5], implants realisation [6–12] and prosthesis design [13–17]. On the one hand, experimental analysis provides further insights into materials characterisation at different scales [18–26], objectifies clinical qualitative outcomes [27–35] and lets numerical models be validated [36]. On the other hand, numerical studies allow studying wider test scenarios and inferring physical quantities otherwise tough to figure out because of feasibility and costs reasons [37–50]. Despite human joint kinematic mechanism and intersegmental forces distribution in human joints have revealed to be a fertile topic for this twofold approach, far too little attention has been paid to investigate load-motion response of mul-ti-level segment of human spine. So far, indeed, experimental studies have mainly reported dynamic responses of single intervertebral joints (functional spinal units [FSU]), providing scant attention to spine’s comprehensive biomechanical behaviour [51,52]. Consequently, even the consistency of numerical multi-segment models has been negatively affected since in silico models can only be partially validated by those local results [53–55]. In this framework, the current work intends to provide a joint experimental and numerical contribution to the representation of multi-level lumbar spine biomechanics. In particular, this study set out to assess experimentally the global rotation of a lumbar segment phantom loaded by flexion-extension moments and to validate the corresponding multi-segmental multibody model, discerning the compressive load effects.

This study performed an experimental and numerical joint analysis, providing a comprehensive insight on the lumbar segment of human spine. The research work consisted of two complementary sections: first, a Sawbones (Sawbones Europe AB, Malmö, Sweden) lumbar spine phantom was characterised in flexion-extension motion and second, a multibody model was designed in MSC Adams environment (MSC Software, Hexagon Corporate Services Ltd., UK). The set up used during experimental tests was re-created in the numerical environment in terms of geometry and load characteristics applied; this way, a consistent comparison between experimental and in silico results was achievable.

The experimental tests were performed on a Sawbones spinal phantom replicating not only the lumbar segment but also the adjacent vertebrae T12 and S1 (SKU3430). The phantom includes the intervertebral discs and the main ligaments; authors well distinguished the anterior and posterior longitudinal ligaments, ligamenta flava, intertransverse ligaments, supraspinal and interspinal ligaments. The latter two were created together.

Flexion and extension tests were realised by applying motion-control linear loadings with a linear-torsion test machine (Instron E3000, Instron Corporation, Norwood, MA, USA). Each test was run at a displacement rate of 20 mm/min aiming to reduce the viscous effects given by the materials representing ligaments and intervertebral discs which could generate adverse effects; the maximum linear displacement of the machine actuator was set at 10 mm. These parameters were settled so that they could result suitable to the aim of the work and, with the intention of future studies involving the phantom, to preserve its integrity. Both flexion and extension tests consisted of five different replicas. For each replica, the initial position of the model was not in contact with the test machine, thus resulting slacked. Consequently, since the experimental set up did not introduce an initial pre-stressed condition, the obtained force-displacement curves were re-aligned at 1N threshold (less than the 2% of the maximum load recorded), compensating the initial experimental noise. The experimental set up configuration is shown in Figure 1: the load was transmitted to the upper extremity through spherical and translational joints, suitably designed to apply a load with a constant arm with respect to the global constraint at S1.

The test machine-phantom couplings were specifically designed in CAD environment and 3D printed (StratasysuPrint SE Plus, Stratasys Ltd., Eden Prairie, MN, USA) to hold both ends of the lumbar phantom; moreover, they were prepared to account for the anatomical space orientation of the lumbar spine. Additionally, each test was recorded with a CANON EOS 5D MarkII camera (Canon Inc. Tokyo, Japan) in order to analyse vertebrae’s kinematic patterns along the sagittal plane. To do so, markers were manually positioned on the centre of vertebral bodies and spinous processes and the planar motion tracking was post-processed by the means of GOMCorrelate software.

The multibody model was designed starting from the vertebral CAD geometry made available by Sawbones, which declared a correspondence with the physical model. Passive elements were added to characterise the model with the same anatomical features of the aforementioned phantom. Averaged density of each vertebra, taking care of both cortical and cancellous bone, was suggested by previous literature works [56–58], and the resulting masses were in accordance with previous in vitro studies [58]. Furthermore, facet joints were not only modelled as contact forces between the vertebral bodies, but attractive forces were also added oriented along the surfaces of the facets in order to limit the relative motion, accordingly to the permitted physiological ones.

Ligaments: The model distinguished the anterior (ALL) and posterior (PLL) longitudinal ligaments, ligamenta flava (LF right-left), supraspinal (SSL) and interspinal (ISL) ligaments and intertransverse ligaments (ITL right-left) for each FSU. The insertion points were positioned based on anthropometric data, and a particular focus was given to those spanning the entire lumbar segment, in order to respect the intrinsic spine’s curvature. In vitro studies showed that spine’s ligaments know a pre-strained condition when the spine is at upright neutral position and characterised by stress–strain curve with two main distinguishable regions: (1) non-linear segment (toe region) and (2) linear segment. Due to these aspects, it is widely spread in multibody modelling literature to reproduce ligaments’ biomechanical behaviour through tension-only force constrained to follow the line of sight between their attachments [57, 59–63]. Accordingly, all the ligaments were described as a pre-tensioned spring element in parallel to a damper. The rest length $l_0$ was inversely computed from the distance between the defined attachment points and the pre-strains $\varepsilon_0$. Table 1 reports the assigned initial strain of ligaments, $\varepsilon_0$, and the strain at the toe-linear regions transaction, $2 \varepsilon_r$; all the values fall within the physiological ranges provided by in vitro studies [62, 64-66]. Starting from the suggested mechanical characterisation given by Putame et al. [2], we developed the following law to describe ligaments' stress-strain (given $\varepsilon=\frac{l-l_0}{l_0}$):

where $k_n$ corresponds to a stiffness per unit strain $[\mathrm{N}]$, calculated from data provided by Pintar et al. [63], $c$ is the damping constant, $l$ and $v_{\text {rel }}$ are the distance and relative velocity of ligament's attachment points, respectively. The choice to go to adimensioned stiffnesses guarantees the removal of possible biases due to different initial lengths between in vitro and in silico data. Finally, little adjustments were performed in order to make each taut length fall in its corresponding toe region range after the dynamic transitory.

Intervertebral discs: The implementation of the 6DOF of the intervertebral discs was obtained by the means of bushing force elements, typically described as non-coupled stiffness and damping matrices. Rotational and translational stiffnesses were initially extracted from reference values [58, 67–69], and then slightly adjusted with respect to the experimental behaviour revealed. Intervertebral discs were supposed at rest in the original spine stance of Sawbones CAD geometry. The orientation of discs’ reference system was set for each FSUin order to comply withits local geometry. Therefore, the longitudinal axis was directed along the conjunction line of the centroids of the adjacent endplates, while the anterior–posterior axis parallel to the line joining the endplate centroid and the spinous process midpoint of the upper vertebra. The subsequent third axis was directly obtained by the previous ones to get a right-handed reference system.

Pre-load: A compressive pre-load aims to simulate the physiological compressive load that the lumbar spine encounters supporting the weight of the upper part of the body. The pre-load was designed in forms of a follower load of 800 N, also satisfying the absence of induced transmission of moment or shear forces to the vertebrae. In the multibody environment, this kind of pre-load was created by compressing each bushing with a force acting along each disc’s longitudinal axis, approximately passing through each FSU’s instantaneous center of rotation [70–73].

Contacts: The only spots where contacts could take place where in correspondence of adjacent vertebrae’s facet joints. Thus, deformable contacts were settled to be consistent with the particular cartilage present therein. The following relationship was adopted:

Where $K$ is the contact stiffness, $\delta$ and $\dot{\delta}$ are the penetration depth displacement and velocity, respectively, $e$ is a non-linear scaling factor, and $C$ is a sigmoid damping function based on $\delta$. Values were extracted by previous studies [53]. Due to the small entities of the displacement applied, contacts appeared during extension only at the end of the motion.

Two multibody models, differing only for the presence of the compressive follower-load, were compared with the experimental results in a twofold way, dynamic and kinematic. On the one hand, flexion-extension moment-angle behaviours of the in silico models were compared to the one resulted experimentally. Experimentally, the moment was obtained from the sagittal vertical force measured by the load-cell of the testing machine multiplied by its constant arm; the angle measured, named $\vartheta_{S 1, T 12}$, was formed by the vertical axis and the line joining the center of mass of S1 (obtained from the multibody and add as a marker on the S1 constraint) and the T12 vertebral body center. The angular displacement was post-processed thanks to the motion tracking tool.

Vertebrae's range of motions (ROM) were considered in a double form: first, we measured the variation of the angle described by the vertical axis and the line joining the centres of mass of adjacent vertebral bodies $\left(\mathrm{CM}_{\mathrm{VB}}\right), \Delta \beta_{\text {rel}}$; second, we measured the variation of the horizontal angles of the lines joining the center of the spinous processes and their corresponding vertebral $\mathrm{CM}_{\mathrm{VB}}, \Delta \varphi_h$. The $\Delta \beta_{\text {rel}}$ provides details about the relative displacements between vertebrae, while the $\Delta \varphi_h$ gives evidence of the variation of orientation of each vertebra.

Experimental tests reveal an optimal reproducibility both in flexion and extension motion- control loads: in flexion, the momentum calculated at the maximum angular displacement (2°) shows a standard deviation of 0.15 N*m less than the 4.4% of the corresponding mean value. Concerning the extension moment, the deviation from one replica to the other is even less (0.05 N*m, 2.1 %). Figure 2 outlines, concurrently with the experimental results, the moment-angle profile of the two in silico models: the one without preload (nFL) and the other subjected to 800 N of follower load (FL). Interestingly, the rachis phantom shows a stiffer response in flexion than in extension along the small angular displacements investigated. Concerning the numerical curves, the follower load has a marked impact as the angle increases. In both cases, the application of the follower-load enables the model to describe better the experimental pattern.

At 2° of flexion, the nFL curve diverges from the experimental mean value by a deviation of 0.8 N*m (23%). Conversely, the FL curve overestimates the experimental results with a lower error, maximum at 2° of bending: 0.31 N*m (9%).

In extension, the FL curve displays a satisfactory moment-angle profile with respect to experimental data: the model load curve fell within the experimental range for most part of the plot, and the maximum separation from it, amounts to -0.16 N*m (-6.9 %). The model without pre-load assumes more deviation from all the other cases as the maximum extension moment differs of -0.71 N*m (- 30.2%).

Finally, the kinematic analysis demonstrates good consistency between the in silico models and the lumbar phantom. Figure 3 and Figure 4 illustrate the $\operatorname{ROM}\left(\Delta \beta_{\text {rel}}\right.$ and $\left.\Delta \varphi_h\right)$ calculated as the difference of the angles at rest configuration and at $2^{\circ}$ bended. Results show a cranio-caudally angular displacement decreasing pattern: upper vertebrae display greater mobility than the distal ones. The high reproducibility of the experimental tests is confirmed by the low variability of kinematic results, almost negligible. Furthermore, extension motion-control loading induces wider variation of angles for upper vertebrae: T12-L1 and L1-L2 $\Delta \beta_{\text {rel}}$ are almost $1^{\circ}$ greater than in flexion. This cannot be said for the lower part of the lumbar segment, where L3-L4 behaviour is almost the same in the two circumstances and L4-L5 demonstrates a more pronounced mobility in flexion ($0.5^{\circ}$ greater). Therefore, the gradient of $\Delta \beta_{\text {rel}}$ ranges along the spine segment reveals being higher in extension than in flexion. The chart in Figure 4, reporting the results of $\Delta \varphi_h$, illustrates interesting aspects: the direction of the motion-control loading seems not varying the orientation variances, except for the extremes of the lumbar segment. Indeed, L1 registers an inclination greater of $0.6^{\circ}$ in flexion than in extension, while L5's pitch angle decreases of almost $0.5^{\circ}$ in extension.

Finally, observing the in silico models, it is possible to assess that the follower loaded multibody model reveals a very strong kinematic accuracy since in almost all cases it provides angular values included in the experimental range. On the contrary, the multibody model not pre-compressed generally shows larger estimates.

The purpose of the current study was to determine the comprehensive behaviour of the lumbar spine phantom when subjected to flexion-extension motion-control load. Moreover, two multibody models were recreated starting from its geometry to understand whether the introduction of a follower load provided greater adherence to the experimental results.

In our study, the load resulted applied to the extremity of the phantom is consistent with the ones applied during in vitro tests [70, 74]. Anyway, our study makes use of a linear motion-control load to execute the torque actions on the phantom, whereas in literature, the adoption of a pure-moment loading is spread [75]: the resulted stiffer behaviour in flexion than in extension is not fully in accordance with in vitro state of art [70, 76–78]. Authors supposed that the reason of this discrepancy could be led by the way of load application, which could influence the behaviour of the spine (motion-control vs load-control). Moreover, the motion was not applied directly to the extremity of the phantom but through a constant arm. Finally, another aspect to keep into consideration is that in [77] the bending stiffness of a lumbar phantom resulted greater in flexion than in extension as highlighted in this study.

Finally, the designing of the numerical multibody models corroborates the effects of the follower load depending on the direction of the bending, already registered in in vitro and in silico studies. Anteriorly, the follower load induces an increasing bending ROM under for a given load, while posteriorly, it is the opposite.

Besides that, the introduction of the follower load makes the multibody model fit to the experimental results, both in terms of moment-angle profile and kinematic analysis and allows the numerical multibody model to be considered validated with respect to the Sawbones phantom. Therefore, starting from the validated numerical model both spinal pathological conditions and post-orthopaedic surgical outcomes can be consistently simulated. This way, future studies will be pursued to objectify their effects and compare them with the corresponding physiological results. To the best knowledge of the authors, this aspect has not been exploited yet in spinal biomechanics literature [79, 80].

None of the authors have any conflicts of interest to disclose.