Musculoskeletal multi-body models (MSK) of the human spine allow for investigations, such as soft tissue characterisation, which would present ethical and practical challenges in vivo (Wang et al., 2021). Additionally, MSK models allow for investigating situations which could present a risk to human subjects (e.g., large loads on the cervical spine) (Silvestros et al., 2019). The study by Silvestros et al. showed a proof-of-concept application to investigate injury mechanisms of the cervical spine using an MSK model (Silvestros et al., 2019). Beaucage-Gauvreau et al. developed a model to allow for the non-invasive investigation of lower back loads during lifting activities (Beaucage-Gauvreau et al., 2019). Other studies have investigated ways of improving surgical treatments by simulating different instrumentation strategies for the same patient (La Barbera et al., 2021). Further, MSK models allow for sensitivity analyses, such as the effect of intervertebral joint (IVJ) location and soft tissue properties on the compressive loads that the vertebrae experience (Senteler et al., 2016; Byrne et al., 2020). They, therefore, have the potential to reduce the risk of spinal injuries and improve spinal surgery outcomes (Beaucage-Gauvreau et al., 2019; Gould et al., 2021; Dall’Ara et al., 2022), and so reduce the cost to individuals and society.

The passive soft tissues between adjacent vertebrae (i.e., the intervertebral joint and ligaments), commonly referred to as the IVJ, have different mechanical properties which contribute to the functional behaviour of the spine and are crucial to maintaining spinal stability (Widmer et al., 2020). MSK models of the spine often simplify the IVJ to three rotational degrees of freedom (DoF) with a fixed centre of rotation (Petit et al., 2004; Han et al., 2012), although more recent studies have included translational DoF (Meng et al., 2015; Ignasiak et al., 2016; Arshad et al., 2017). The mechanical properties of the IVJs are not included in all models and if present are simplified to a lumped parameter spring-damper model, commonly referred to as a bushing force (Silvestros et al., 2019; Wang et al., 2019; 2021; Alizadeh et al., 2020; Gould et al., 2021). Parameters for the bushing force are taken from the literature (Han et al., 2012; Christophy et al., 2013; Ignasiak et al., 2016; Senteler et al., 2016; Arshad et al., 2017; Zhang et al., 2020; Remus et al., 2021; Lerchl et al., 2022) and applied in each DoF (Silvestros et al., 2019; Wang et al., 2019; 2021; Alizadeh et al., 2020; Gould et al., 2021). Despite these simplifications, representing the IVJ with a bushing force is complex due to the interaction between joint pose and stiffness which influences the predicted joint loads and muscle forces (Byrne et al., 2020). Furthermore, although rarely accounted for, the IVJ properties vary strongly between specimens and spine levels in all DoF (Panjabi et al., 1976; Andersson and Schultz, 1979; Miller et al., 1986; McGlashen et al., 1987; Izambert et al., 2003; Petit et al., 2004; Renner et al., 2007; Garges et al., 2008; Doulgeris et al., 2014; Schmidt et al., 2015; Newell et al., 2019; Palanca et al., 2020) and influence predicted muscle forces, intervertebral disc loads, and range of motion (Arshad et al., 2017; Wang et al., 2019; Byrne et al., 2020). Moreover, several studies have highlighted the importance of determining subject- or specimen-specific and condition-specific stiffnesses which can result in more accurate kinematic predictions (Ghezelbash et al., 2016; Alizadeh et al., 2020; Wang et al., 2020; 2021). Therefore, accurate subject-specific characterisation of the IVJ is necessary to prove the reliability of models.

Personalisation of the IVJ stiffness has been done using hybrid models (Remus et al., 2021) and by using finite element models which are integrated into MSK models (Byrne et al., 2020). Another method to estimate subject-specific properties has been to scale the stiffness based on anthropometric data (Ghezelbash et al., 2016). The meta-analysis of studies of the rotational behaviour of human cadaveric spine segments by Zhang et al. has established a regression model of the moment-rotation behaviour of the IVJ (Zhang et al., 2020) which could be used to assign stiffnesses to the IVJ in MSK models. However, few studies have directly optimised the stiffness within MSK models. Wang et al. developed a generalised stiffness model from literature data of the IVJ which was incorporated into an MSK model (Wang et al., 2020), using this model they went a step further, optimising a subject-specific stiffness based on in vivo motion (Wang et al., 2021). One study has used a genetic algorithm with MSK simulations to optimise the IVJ stiffness of the model for ex vivo porcine cervical spine specimens (Silvestros et al., 2019). These studies considered the level dependency of the IVJ stiffness at the individual IVJ levels (Silvestros et al., 2019; Wang et al., 2021). To simulate spinal surgeries Petit et al. developed a method using functional bending tests to optimise patient-specific stiffnesses, however, this approach simplified the level dependency by dividing the thoracolumbar spine into three regions (Petit et al., 2004).

Stiffnesses have been personalised in flexion/extension and lateral bending under flexion/extension motion and lateral bending motion using motion capture data from in vivo experiments, personalised stiffnesses in axial rotation were not reported (Wang et al., 2021). Given the direction-dependent nature of the IVJ stiffness (Miller et al., 1986; McGlashen et al., 1987; Renner et al., 2007) further work is needed to predict stiffnesses under axial rotation. Furthermore, optimising the stiffnesses in the rotational DoFs required an estimation of the joint kinematics using an optimisation algorithm as the motion capture data came from in vivo experiments (Wang et al., 2021). Wang et al. suggested that with more accurate motion tracking data, stiffnesses could be optimised without calculating the joint kinematics through optimisation, this is possible with data from ex vivo experiments (Silvestros et al., 2019) but has not been done for rotational DoFs.

While previous studies have shown optimisation of the subject or specimen-specific stiffnesses improves the accuracy of the spinal MSK models (Silvestros et al., 2019; Wang et al., 2021), the aims of the previous studies were not to investigate the inter-specimen variability or the variation between spine levels of optimised stiffnesses. A study focusing on the inter-specimen and spine level variability which optimises stiffnesses within MSK models is lacking and would evidence the need for specimen-specific optimisation which accounts for IVJ level dependency and loading conditions. Moreover, a study using highly accurate motion tracking in combination with MSK models to simultaneously optimise the stiffnesses of the lumbar spine in all rotational DoF for multiple loading conditions would strengthen and build upon the existing studies.

Therefore, the current study aimed to investigate the variability of optimised stiffnesses across different specimens and different spinal levels in all rotational DoF through numerical simulations using accurate motion capture data collected in vitro from human specimens to allow for highly accurate tracking. This also included investigating these aspects under different loading.

This study reanalysed a subset of the imaging, load, and motion capture data from the results of the in vitro biomechanical study by Volkheimer et al. (Volkheimer et al., 2018). More specifically, the experimental study (Volkheimer et al., 2018) provided data for six T12-sacrum spine segments from human cadavers (Table 1). The specimens were acquired from the Science Care (United States) donation program. In brief, the specimens were cleaned of soft tissue leaving the vertebrae, intervertebral discs, all ligaments, and facet joints intact. The sacrum was completely constrained and pure moments were applied to T12 using a spine tester (Wilke et al., 1994). Three loading conditions were individually applied, flexion/extension, lateral bending, and axial rotation. Each loading condition was applied as a pure moment by means of a gimbal and stepper motor (Wilke et al., 1994). At all joints there were 6 DoF, therefore although a pure moment was applied, the nature of spinal biomechanics meant that the joints experienced coupled moments. 3.5 loading cycles were applied at 1.0°/s in flexion/extension and lateral bending and 0.5°/s in axial rotation with a peak moment of 7.5 Nm. The loading rates were selected to reduce the influence of creep and inertia on the results and the load magnitude was based on recommendations within literature (Wilke et al., 1998). The load cycle for evaluation was selected as the one with the smallest range of the coupled moments. The first cycle was excluded a priori as the initial recorded load was not always 0N. The moment in the loaded DoF and the resulting coupled moments in the non-loaded DoF were measured using a six-component load cell mounted above the specimens (FT 1500/40, Schunk GmbH, Lauffen/Neckar, Germany) (Figure 1A). Three reflective markers were attached to the anterior surface of each vertebra (Figure 1B) (Volkheimer et al., 2018). The vertebral motion was measured with a motion tracking system (Vicon MX13+, Vicon Motion Systems Ltd., Oxford, UK) including six infrared cameras. From this, the rotation in each DoF at each joint was calculated, the measured rotations were in agreement with other in vitro studies (Panjabi et al., 1994). These rotations were then provided as inputs for the current computational study. Additionally, the specimens were imaged twice:

• CT scans with a pixel size of 0.39 × 0.39 mm and slice increment of 0.5 mm (thickness 1.0 mm) were obtained with a Brilliance 64, Philips CT device using a voltage of 120 kVp and a tube current of 356 mA,

First, a summary of the optimised stiffnesses and kinematic errors of the three types of models (literature, optimised, and cross-validation) considering the loading direction is presented. Then the overall inter-specimen variability for both the uniform and level-dependent models is discussed. Next, the variation between spine levels is discussed, first addressing the uniform models and then the level-dependent models. Finally, a direct comparison is made between the results of the optimised uniform models and the level-dependent models.

This study aimed to investigate the inter-subject variation and the difference between spinal levels under different loading conditions. Six specimen-specific models were constructed and employed to simulate left lateral bending, left axial rotation, and flexion experiments. The IVJ stiffness was optimised with an optimisation algorithm which minimised the predicted motion error. The optimisation was performed considering (i) no level dependency, and (ii) level dependency implemented as a fixed ratio between the IVJ levels.

This study built on previous studies that have sought to determine subject-specific stiffness properties in two degrees of freedom through optimisation (Silvestros et al., 2019; Wang et al., 2021) by simultaneously optimising the stiffness in all three rotational degrees of freedom simultaneously. Rather than using data from in vivo experiments which require optimisations to calculate the joint kinematics, this study used ex vivo dataset of the lumbar spine which provides more accurate motion tracking data and removing the need for an optimisation to calculate the joint kinematics (Wang et al., 2021). Previous studies have effectively used ratios to distribute the motion across the IVJ levels (Bruno et al., 2015; Burkhart et al., 2020; Wang et al., 2021), the present study attempted to distribute the stiffness properties across the IVJ levels with a similar approach which to the best of the Authors’ knowledge has not previously been attempted.

Within the general population, a large variation of intervertebral joint stiffness is expected between individuals. Optimising stiffnesses using the tracking error reflected this, as a wide range of optimised stiffnesses were calculated (Figure 5). These stiffnesses fell within the range of experimentally reported stiffnesses (Andersson and Schultz, 1979; McGlashen et al., 1987; Garges et al., 2008; Palanca et al., 2020). This implies that the proposed optimisation approach was effective at identifying specimen-specific stiffnesses and able to capture the inter-specimen differences.

For the models without level dependency, the optimisation tended to average out the errors across the spine (observed when generic stiffness values were implemented) and the resulting stiffness was representative of an overall spinal stiffness for each specimen in each loading condition. Thus, similar kinematic errors were observed across all spine levels. Despite assuming no level dependency, substantial improvements were seen between the literature and optimised stiffnesses for each specimen (Figures 6–8). Similarly, when using level-dependent stiffnesses the optimiser was attempting to average out the errors across the levels; however, the distribution of the stiffnesses across the joint levels was defined by a fixed ratio. The level-dependent models with optimised stiffnesses were associated with markedly smaller tracking errors (Figures 6–8) compared to level-dependent models with literature stiffnesses. This supports the findings of the study by Wang et al. (Wang et al., 2021), that specimen-specific stiffness (level-dependent stiffnesses or uniform) is needed to accurately predict the kinematics. Furthermore, the presence of significant differences in the accuracy of the cross-validation models supports the need for specimen or subject-specific properties. Additionally, the variation of the kinematic accuracy between optimised stiffness models indicates factors other than the stiffness should also be considered, such as joint pose (Senteler et al., 2018; Byrne et al., 2020).

The optimisation of the uniform stiffnesses across spinal levels (Models_{Optim_Unif}) resulted in similar errors at all spinal levels for some specimens while for other specimens there were high errors (>2°) at some levels and small errors (<0.1°), although statistically significant differences were not found. This could suggest the extent of the level dependency varies between specimens.

The extent of the motion is largely dependent on the IVJ stiffness, which has been reported to vary between spine levels (Pintar et al., 1992; Renner et al., 2007). The initial stiffnesses were based on literature data for the L3L4 IVJ level, yet the errors at the L3L4 IVJ level were comparable to the errors at the other levels when looking at the generic stiffnesses (Models_{Lit_Unif}). This suggests that using a generic but level-dependent stiffness may not offer any improvement over a generic uniform stiffness in the prediction accuracy. In the present study, introducing a generic but level-dependent stiffness resulted in higher errors than when a uniform stiffness was used (Figure 9). Therefore, the introduction of level dependency via a fixed scaling factor is not suitable. However, this could also be due to the specific scaling factors used which were more suitable for some joint levels (resulting in very low prediction errors, <0.5°) and less so for other joint levels, typically the extreme joints (resulting in comparatively larger errors). The scaling factors were calculated from data presented in the meta-analysis by Zhang et al. (Zhang et al., 2020) (maximum rotations without a compressive preload), and the moments reported were higher than the moments used in the experiments the current study simulated. Therefore, if a ratio is used to describe the variation of the stiffness between levels, it needs to be calibrated for the loading conditions.

Considering the influence of the loading conditions, the predicted stiffnesses varied greatly for the different loading conditions. The predicted kinematic errors were markedly higher for the level-dependent models with literature stiffness under a lateral bending load, this implies the scaling ratio introduced to define the level dependency was more inappropriate for lateral bending than for axial rotation and flexion. This could indicate that the level dependency changes with the loading condition.

For all loading conditions, the stiffnesses were optimised simultaneously in all three rotational DoF for each loading condition. This is necessary as the spine exhibits coupled behaviour (Wu et al., 2014). The study by Meng et al. found that either coupled or uncoupled stiffnesses could represent the spinal properties however coupled and uncoupled stiffnesses should not be used interchangeably (Meng et al., 2015). In the present study, the spinal stiffnesses were modelled as uncoupled, resulting in a limited accuracy of the models in the DoF in which the load was not applied. Although the stiffnesses were optimised in these DoF their value did not change much relative to the initially set values, which were based on the literature (Table 3, Supplementary Data SB). This could be due to the motion being smaller in the unloaded directions, thus reducing the sensitivity of the cost function to these directions. Future research should address this limitation by either introducing coupling terms or by identifying a cost function that results in a better representation of the spinal stiffness in the DoF in which the load is not applied.