Establishment of the finite element model
The cervical finite element model was established based on a previously validated model built by Deng et al [27]. This model was constructed using CT data from a healthy Chinese male volunteer (age 30 years, height 170 cm, body weight 68 kg). All patient gave informed consent and Chinese Ethics Committee of Registering Clinical Trials approved the study. The model incorporated various components such
as vertebrae, intervertebral discs, facet joints, and ligaments. Specifically, the vertebral model encompassed cortical bone, cancellous bone, end plates, and end plate cartilage. The average thicknesses of the cortical bone, end plate, and end plate cartilage were 1.5 mm, 0.5 mm, and 1.0 mm, respectively. The intervertebral disc was composed of nucleus pulposus, annulus ground substance and reinforced fibrous membrane. The annulus fibrosus matrix was composed of five layers of hexahedral elements, while the reinforced fibrous membrane consisted of five pairs of concentric shell grids. The fiber angles were set to ± 25° ~ 45°, which gradually changed from the outside to the inside annulus fibrosus along the radial direction. Spring elements were used to simulate the ligaments, which limit segmental range of motion. The model included the anterior longitudinal ligament, posterior longitudinal ligament, ligament flavum, alar ligament, supraspinous ligament, interspinous ligament, intertransverse ligament, transverse atlantooccipital ligament, posterior atlantooccipital membrane, posterior atlantooccipital ligament and articular capsule ligament. Facet joints in the model included articular cartilage and articular capsule ligament. The articular cartilage attached to the facet was simulated by a layer of hexahedron elements, and the articular capsule ligament was simulated by spring elements. The intact model (C2-T1) included 254484 nodes and 345367 elements, as shown in Fig. 1.
All methods were conducted in accordance to relevant guidelines and regulations.
Constitutive model and parameters
The constitutive model and parameters utilized in this study for simulation are listed in Table 1 [28, 29] and Table 2. [30,31,32,33] The bone was modeled by isotropic plastic material. The strength of the end plate was defined as 1/3 of the cortical bone. The cartilage and vertebra were connected with the same node and modeled using elastic material. The nucleus pulposus and annulus fibrosus matrix were represented using viscoelastic and elastic materials, respectively.
Calibration and validation
All six degrees of freedom of all nodes of the inferior end plate of T1 were constrained as boundary conditions. A reference point was established 2 mm above the top of the odontoid process, which was located at the rotation center of the cervical spine. The reference point was coupled with all element nodes of the C2 odontoid process. According to the right-hand law, ± 0.5 Nm, ± 1 Nm and ± 2 Nm torques were applied to the reference point to simulate six movements, including cervical flexion, extension, left and right rotation and left and right bending. Figure 2 illustrates a comparison of the ROM in the current FE model with similar studies conducted in vitro and using computational methods, specifically focusing on flexion and extension movements. [34,35,36]
Establishment of surgical models
We modified the C5/6 segments of the validated intact model to construct PECFD models. In the PECFD model, a portion of the right facet joint in the C5/6 segment, along with sections of the articular capsule ligament, posterior longitudinal ligament, and ligamentum flavum, were removed. Then, ICEM CFD was used to remesh the C5/6 segment and construct the PECFD model. The post-operative skeleton hard tissue model and corresponding schematics are shown in Fig. 3.
Models of PECFD with three different extents of disc resection area were established (Fig. 4). [37,38] The first PECFD model considered the case in which the annulus fibrosus was bit off but the nucleus pulposus was not involved (Model 1, M1). In the second case, the annulus fibrosus was bit off, and the margin of the nucleus pulposus was involved (Model 2, M2). The third case considered the condition in which part of the nucleus pulposus (less than 50%) was bit off (Model 3, M3).
CMT simulation
One coupling node was placed 1 mm directly above the odontoid of C2, and the upper surface of C2 and the coupling point were kinematically coupled. The lower end of T1 was fixed. As shown in Fig. 5, a load of 55 N was applied on the coupling point to simulate head weight, and then a moment of 1 Nm was applied on the coupling point to simulate cervical flexion, lateral bending, and axial rotation. As shown in Fig. 5a, another two coupling nodes were placed at the anterior 70 mm of C4 and the 5 mm behind the spinous process of C5. This is to simulate the force acting on the chin in the massage technique, which propagates through the head to the surface of C2. These two nodes were used as the points for applying forces produced by CMT. This study considered both right and left sides of manipulation directions. In the first case, a force of 68 N to the right was applied on the C4 coupling point, and a force of 32 N to the left was applied on the C5 coupling point to simulate the CMT in the right direction and was referred to as load1. This is to simulate the effect of surface force on the skin behind the neck in the manipulation. As shown in Fig. 5b, the directions of the above two forces were reversed to simulate the CMT in the left direction, as referred to as load2.
This study measured and evaluated indexes of cervical stability and maximum stress concentration, including the cervical maximum displacement (MD), cervical rotational angle (RA), maximum stress of annulus fibrosus matrix, intervertebral disc pressure (IDP) and stress of facet ligament (FL). Displacement of the maximum displacement node of the entire cervical spine referred to the distance from the tip of the C2 cervical odontoid process to the original position. The rotation angle was projected to each coordinate axis plane to obtain the angle of extension, rotation and lateral bending when the C2 vertebra moved to the maximum position.