The processing and extraction of patient data described in the following section were executed from asynchronously phantom-based calibrated CT image data (Kaesmacher et al., 2017). We labeled and segmented vertebrae using an automated deep learning-based process for vertebrae segmentation, which is described in detail elsewhere (Sekuboyina et al., 2020). In brief, three artificial neural networks (ANNs) are used to 1) detect the spine, 2) identify and label each vertebra as well as 3) segment each vertebra based on the label. The latter two steps were reviewed by a radiologist and could potentially be corrected. For each vertebra, centroids, as well as segmentation masks for eleven subregions, were generated using a fourth ANN: vertebral body (further divided into the cortex and the trabecular compartment), vertebral arch, spinous process, as well as transverse processes. Before calculating necessary points of interest, the centroids of the first thoracic and last lumbar vertebra were aligned vertically to account for posture differences between supine from CT scans and upright position. Thereafter, these data were used to calculate relevant points for muscle and ligament attachments. Figure 1 shows the overall process.

Subsequently, we defined points of interest for each individual vertebra. Therefore, the algorithm iterated over each vertebra, creating bounding boxes based on its binary segmentations. Corresponding to those bounding boxes, individual subregion segmentations were used to define landmarks for muscle and ligament attachment points by geometrical extreme values as shown in Figure 2A. Thus, depending on the subregion, the most posterior, inferior, superior, or lateral point on the surface was determined by its minimal and maximal coordinate values along the corresponding spatial axis. Based on centroid positions, auxiliary sagittal and horizontal planes were set through the vertebral body to extract its attachment points (Figure 2B). In the horizontal sectional plane, the most lateral points on each side of the vertebral body were extracted. In the sagittal sectional plane, a rectangle is fitted around the subregion of the vertebral body. The corner and center points of the rectangle border were then projected onto the surface by the shortest distance. Using a similar function, attachment points on the vertebral arch were determined via the minimal distance between the anterior border of its sectional plane and posterior points in the vertebral body. The plane was then shifted right and left and the process was repeated.

We assumed the location of the intervertebral joint to be the midpoint of the straight line connecting the central points of the lower and upper endplates of the two vertebrae representing one motion segment (Figure 3). We used a spline interpolation of all vertebral body centroids to define intervertebral joint orientation by calculating the spline derivative at the intersection with the upper endplate of the inferior vertebra.

To account for individual torso weight distribution, a segmentation mask for lung, fat, and muscle/organ tissue was created based on typical CT intensity ranges. Subsequently, the segmented tissue was assigned to the nearest vertebra depending on its minimal anterior–posterior distance. Thus, the torso weight was subdivided into segments for each vertebral level. For each segment, the algorithm calculated its center of mass and total weight corresponding to its individual tissue distribution. We assume to have an average density of 0.25 g/cm³ for the lung, 0.96 g/cm³ for fat, and 1.06 g/cm³ for the remaining soft tissues (Pearsall et al., 1996; Akhavanfar et al., 2018).

The automated modeling of the thoracolumbar spine with generic bodies of the head–neck complex, ribcage, simplified upper extremities, and the pelvic-sacral region is carried out using the multibody simulation software SIMPACK (Dassault systèmes, France). The thoracolumbar spine includes individual information on vertebrae T1-L5, insertion points for muscles and ligaments, spinal alignment, as well as paraspinal soft tissue distribution as described in the previous chapter. Bodies for the head–neck system, ribcage, sacrum, and pelvis are individually scaled according to Winter (2005) and equipped with relevant points for muscle insertion and integrated into the model. Neglecting facet joints and intraabdominal pressure, lumbar intervertebral joints are modeled as actuated spherical joints to ensure necessary stability. The thoracic spine, head–neck, and ribcage are modeled as one rigid body, and segment masses for soft tissue are rigidly fixed to each vertebra according to the calculated centers of mass. Segment masses relative to overall torso mass calculated from the CTs (Winter, 2005) were assigned to the bodies for head–neck and simplified arms. The masses of bony structures were calculated assuming a density of 1.5 g/cm³ (Pearsall et al., 1996; Akhavanfar et al., 2018). Intervertebral discs, as well as ligaments, are modeled as nonlinear, viscoelastic force elements. Occurring moments in the intervertebral discs are characterized by a nonlinear load–deformation relationship (Rupp et al., 2015; White, 2022). Specific data on stiffness and neutral zones of the intervertebral discs were taken from White (2022). The model includes anterior and posterior longitudinal ligament (ALL and PLL), flavum ligament (LF), interspinal ligament (ISL), and supraspinal ligament (SSL). The characteristic force–length curve for the elastic behavior of ligaments (Figure 4) shows a nonlinear toe region in the region of small deformations (A), followed by a linear elastic region before the final failure of the ligament (B). Inspired by the study cited herein(Rupp et al., 2015), their nonlinear force–length characteristics are described in Eq. 1. F e l l lig = 0 , for l lig ≤ l lig,0 K n l l lig − l lig,0 e x p nll , for l lig ≤ l lig,A F A , n + K lin l lig − l lig,A , for l lig > l lig,A (1) with

Individual lengths for ligament segments are measured directly in the model in neutral position, which is what we considered upright standing. Values for initial ligament lengths were calculated considering values for pre-strain from the studies cited herein (Nachemson and Evans, 1968; Aspden, 1992; Robertson et al., 2013). The main difference to the mechanical law provided by Rupp et al. is that we take relative strain values for l_{lig, A} and l_{lig, B} instead of absolute elongations. Therefore, we guarantee uniform preloads within one ligamentous structure for the neutral position. Parameters for strains ϵ _A,ϵ _B, and Forces F _A and F _B at points A and B, are taken from the study referred herein (Chazal et al., 1985).

Nine muscle groups of the lower back and abdomen are included in the model as 103 point-to-point actuators: rectus abdominis (RA), internal obliques (IO), external obliques (EO), psoas major (PM), quadratus lumborum (QL), multifidus (MF), longissimus thoracis pars lumborum (LL), iliocostalis lumborum (IL) and the interspinales lumborum (IS) (Figure 5). Globally acting muscles RA, IO, and EO, as well as those muscle fascicles of LL, QL, and IL attached to the ribcage, are simplified each to one actuator per side. Muscles acting locally on the lumbar spine are modeled in detail based on attachment points taken from a cadaver study (Bayoglu et al., 2017a). Muscle fascicles attached to the same subregion were combined and physiological cross-sectional areas (PCSA) based on Christophy et al. (2012) were assigned to respective fascicles.

The MBS model calculates necessary joint moments M to hold the imposed static positions, which are transferred to Matlab via a SIMULINK model, where muscle forces are calculated using a static optimization approach (Gagnon et al., 2001). In order to increase the chances of finding a global optimum, we used the globalsearch solver (Global Optimization Toolbox, Matlab 2020b) to solve the following optimization problem: M i n i m i z e C o s t F u n k t i o n = ∑ i = 1 n F i P C S A i 3 (2) subject to equality constraints c e q = M x , 1 … M z , i = 0 (3) and bound constraints 0 ≤ F i ≤ σ m a x P C S A i (4)

Focusing on vertebral loading in the sagittal and frontal plane, equality constraints consider respective moments (x for frontal, z for sagittal) occurring in each lumbar intervertebral joint. Only active forces are taken into account neglecting the passive elastic behavior of muscular tissues. To guarantee compliance with the equilibrium conditions for all load cases, maximal muscle stress (MMS) was first set to 0.6 MPa, (Arjmand et al., 2009), and then to 1 MPa (Bruno et al., 2015; Beaucage-Gauvreau et al., 2019; Favier et al., 2021).

For model validation, we built models based on two nondegenerative spine datasets (Figure 6). We selected datasets of two young patients (1 M, 1 F) with anthropometric data as comparable as possible to the subjects in comparative studies (Wilke et al., 2001; Takahashi et al., 2006; Rohlmann et al., 2008).

Anthropometric data, such as body height and weight, was calculated in reference to individual spine height and torso weight from CT data according to the study cited herein Winter (2005). To characterize individual spinal alignment, we measured kyphosis and lordosis angles for T1-T12, and L1-S1 in the sagittal plane, as well as Cobb angles for C7-T12, in the frontal plane and these are summarized in Table 1.

We evaluated predicted muscle forces and spinal loading for various activities. The load cases were selected based on previous studies on the measurement of intradiscal pressure and muscle activation (Wilke et al., 2001; Rohlmann et al., 2008). Table 2 summarizes all investigated load cases.

Since the subjects from in vivo studies were all men (Table 2), we only used the model based on the dataset of the male for validation. Prior to the simulation of dedicated load cases, we used an optimization routine to determine the optimal position for upright standing. Assuming the optimal position to be energy efficient, joint angles of the lumbosacral spine (T12—Sacrum) were optimized, subject to minimization of occurring joint moments in the sagittal plane. We adopted determined joint angles as the starting posture for a neutral stance in the further course.

Applied flexion was assumed to be 40% sacral rotation and 60% lumbar flexion (Liu et al., 2018). According to the studies referred herein (Wong et al., 2006; Christophy et al., 2012), lumbar flexion was distributed 25.5% for L1/L2, 23.1% for L2/L3, 20.4% for L3/L4, 18.5% for L4/L5, and 12.5% for L5/S1.

Ligament modeling was evaluated based on their stress states during each load case with a focus on whether they were within a physiological range or whether failure could already be expected. Calculated lumbar loads were compared to respective vertebral load measurements and predicted muscle forces were evaluated in the context of measured EMG signals from experimental studies.

To determine the potential of our pipeline regarding the systematic investigation of individual spinal loads, we compared generated results for both patients considering their individual spinal alignment. For this, compensation angles for an upright position, as well as estimated muscle forces and spinal loads, were compared, put into context with the identified curvature of each subject in the frontal and sagittal plane, and analyzed for plausibility.

For initial simulations with LF pre-strain taken from the study referred herein (Nachemson and Evans, 1968), the occurring LF forces exceeded physiological maximum forces even for low-intensity flexions ( < 10°). We reduced pre-strain from 10% to 5%, which still lies within the standard deviation of experimentally determined values (Nachemson and Evans, 1968) (Figure 7). Mean normalized forces were within a physiological range for all ligaments during investigated loading tasks. However, despite adjusted pre-strain, the average LF force reaches close to 100% with standard deviations of up to 37% during 30°flexion. The remaining ligament showed low (ALL, PLL, SSL) to moderate (PLL, ISL) loading, with a less than 2% and 50%, respectively.

The estimated muscle forces of the erector spinae correlated closely (r = 0.95) with measured EMG-signals from Takahashi et al. (2006) (Figure 8). Initial simulations with a detailed muscular architecture according to Bayoglu et al. (2017a), under consideration of physiological MMS 60 N/cm², could not satisfy equilibrium conditions for all models, even for moderately intense activities, such as 30° flexion. Even increasing the MMS to 1 MPa was not sufficient to reliably satisfy the equilibrium conditions for all cases, though this mainly affected high-intensity cases, such as extensive flexion with additional weight. Therefore, we adapted muscle properties according to a validated musculoskeletal model from literature (Christophy et al., 2012). However, Christophy’s model includes no muscle fascicles for the IS, nor those fascicles of the MF attached to the thoracic spine, which we added and in order to guarantee compliance with the equilibrium conditions, equipped with PCSAs of 1 cm², which is rather at the higher end of the range of measured values for MF fascicles (Bayoglu et al., 2017a).

Overall, the estimated compression loads on intervertebral joints for various static loading tasks showed a good correlation with reported spinal load measurements (r = 0,98) (Figure 9). However, our model slightly underpredicted compression forces normalized to upright standing compared to measured forces from the study cited herein (Takahashi et al., 2006) by up to 16%. Normalized compression forces for upright standing with 20 kg weight held 55 cm from the Sacrum were overestimated by 33% compared to data from the respective comparative study (Wilke et al., 2001) (Figure 10).

Our results demonstrated that our patient-specific models are well suited to investigate interindividual differences. In correlation with individual spinal characteristics, distinct differences in muscle activity and vertebral loading between investigated subjects could be identified. Thus, subject two showed generally higher muscle activity than subject one (Figure 11). The compensation angles for the upright standing position were in the same range. In both cases, the balancing movement started to a large extent from a sacral rotation. However, subject one showed a rather extensive compensation in the lower mid-region of the lumbar spine, whereas subject two compensated within a smaller range in the upper lumbar region (Figure 12). Subject one showed a smaller overall compensatory flexion with -3,9° than subject two (-4,5°). The estimated shear forces showed considerable differences, especially regarding anterior–posterior loading of up to 123 N anterior and 84 N posterior for subjects one and two, respectively. Please note that this relates well to the differences in lumbar lordosis between both subjects. In terms of lateral shear forces, the differences in the thoracolumbar transition are particularly noticeable and considerably higher for the subject with thoracic scoliosis.