Development of a Detailed Human Spine Model with Haptic Interface

The human spine is an important structure, providing strength and support as well as permitting the body to bend, stretch, rotate and lean. It is also a vulnerable part of our skeleton open to disease and injuries like whiplash, low back pain, scoliosis etc. Whiplash is a frequent consequence of rear-end automobile accidents and has been a significant public health problem for many years with injury to one or more ligaments, intervertebral discs, facet joints or muscles of the neck. Low back pain is a very common disease and strongly associated with degeneration of intervertebral discs (Luoma et al., 2000) seen in people with sedentary jobs spending hours sitting in relatively fixed positions, with the lower back forced away from its natural lordotic curvature. This causes health risks of the lumbar spine, especially for the three lower vertebrae L3-L5. 80% of people in the United States will have lower back pain at some point in their life (Vallfors, 1985). Scoliosis is a less common but more complicated spinal disorder, being a congenital 3D deformity of the spine and trunk affecting 1.5-3% of the population. In severe cases, surgical correction is required to stabilize the scoliosis curvature. Studies into the treatment of these spinal diseases have played an important role in modern medicine. Many biomechanical models have been proposed to study dynamic behavior and biomechanics of the human spine, to develop new implants and new surgical strategies for treating these spinal diseases.


Introduction
The human spine is an important structure, providing strength and support as well as permitting the body to bend, stretch, rotate and lean. It is also a vulnerable part of our skeleton open to disease and injuries like whiplash, low back pain, scoliosis etc. Whiplash is a frequent consequence of rear-end automobile accidents and has been a significant public health problem for many years with injury to one or more ligaments, intervertebral discs, facet joints or muscles of the neck. Low back pain is a very common disease and strongly associated with degeneration of intervertebral discs (Luoma et al., 2000) seen in people with sedentary jobs spending hours sitting in relatively fixed positions, with the lower back forced away from its natural lordotic curvature. This causes health risks of the lumbar spine, especially for the three lower vertebrae L3-L5. 80% of people in the United States will have lower back pain at some point in their life (Vallfors, 1985). Scoliosis is a less common but more complicated spinal disorder, being a congenital 3D deformity of the spine and trunk affecting 1.5-3% of the population. In severe cases, surgical correction is required to stabilize the scoliosis curvature. Studies into the treatment of these spinal diseases have played an important role in modern medicine. Many biomechanical models have been proposed to study dynamic behavior and biomechanics of the human spine, to develop new implants and new surgical strategies for treating these spinal diseases.

Related works
Biomechanic models can be divided into four categories: physical, in-vitro, in-vivo and computer models. Computer models have been extensively used because they can provide information that cannot be easily obtained by other models, such as internal stresses or strains. They can also be used for multiple experiments with uniform consistency, lowering experimental cost, and simulating different situations easily and quickly. Multi-body and finite element models, or a combination of the two are the most popular simulation tools that can contribute significantly to our insight of the biomechanics of the spine.

www.intechopen.com
Finite element models (FEMs) are helpful to understand underlying mechanisms of injury and dysfunction, leading to improved prevention, diagnosis and treatment of clinical spinal problems. FEMs often provide estimates of parameters that in-vivo or in-vitro experimental studies either cannot or are difficult to obtain accurately. FEMs are divided into two categories: dynamic study and static study models. Models for static study generally are more detailed in representing spinal geometries (Greaves et al., 2008;Kumaresan et al., 1999;Natarajan et al., 2007;Teo & Ng, 2001;Yoganandan et al., 1996). They can predict internal stresses, strains and other biomechanical properties under complex loading conditions, but they generally only consist of one or two motion segments and do not provide effective insight for the whole column. Dynamic study models generally include a series of vertebrae connected by ligaments and disks modeled as springs (Maurel et al., 1997;Ng & Teo, 2005;Pankoke et al., 1998;Schmidt et al., 2008;Zander et al., 2002). They predict locally the kinematic and dynamic responses of a certain part of the spine under load. FEMs have also been widely used to study scoliosis biomechanics (Aubin, 2002;Gignac et al., 2000;Lafage et al., 2004;Schlenk et al., 2003). Understanding of the biomechanics of spine deformation helps surgeons formulate treatment strategies for surgery and design and develop new medical devices. Being complex, FEMs of scoliotic spines are usually restricted to 2D models or simplified 3D elastic beam element models. Although preliminary results achieved are promising, extensive validation is necessary before using the models clinically.
Multi-body models have advantages like less complexity, less demand on computational power, and relatively simpler validation requirements. They possess potential to simulate both kinematics and kinetics of the spine effectively. Rigid bodies are interconnected by various joints. Multi-body models can also include many anatomical details while being computationally efficient. The head and vertebrae are modeled as rigid bodies and soft tissues are usually modeled as massless spring-damper elements. Such models are capable of producing biofidelic responses and can be broken down into two categories: car collisions and whole-body vibration investigations. In the former, displacements of the head with respect to the torso, accelerations, intervertebral motions, and neck forces/moments can provide good predictions for whiplash injury (Esat & Acar, 2007;Garcia & Ravani, 2003;Jager et al., 1996;Jun, 2006;Linder, 2000;Stemper et al., 2004;Van Der Horst, 2002). In the latter, multi-body models are helpful for determining the forces acting on the intervertebral discs and endplates of lumbar vertebrae (Fritz, 1998;Luo & Goldsmith, 1991;Seidel & Griffin, 2001;Verver et al., 2003;Yoshimura et al., 2005). In both cases, multi-body models focus either on the cervical or the lumbar spine. Since spine segments are partially modeled in detail, it is impossible to investigate the kinematics of the thoracic spine region. In other words, global biodynamic response of the whole spine has not been studied thoroughly.
Although finite element and multi-body models are powerful tools used to study intrinsic properties of injury mechanisms, other techniques have been developed to obtain deeper understanding of biomechanical properties of medical diseases. One such technique is computer haptics. The word haptics was introduced in the early 20 th century to describe the research field that addresses human touch-based perception and manipulation. In the early 1990s, the synergy of psychology, biology, robotics and computer graphics made computer haptics possible. Much like computer graphics is concerned with rendering visual images, computer haptics is the art and science of synthesizing computer generated forces to the user for perception of virtual objects through the sense of touch. Simulation with the addition of haptics may offer better realism compared to only a visual interface. In recent www.intechopen.com years, haptics has been widely applied in numerous VR environments to increase the levels of realism. Especially, haptics has been investigated at length for medical education and surgical simulations, like surgical planning and laparoscopic surgical training (Basdogan et al., 1998;Forest et al., 2004;Gorman et al., 2000;Seitz et al., 2004;Williams et al., 2004).
Haptics has been widely utilized in medical fields, but little has been applied to spinal diseases. Integrating haptics into spine models means surgeons can investigate biodynamic responses of whole human spine which either have not been investigated enough in the literature or are limited to partial spine segments. Understanding biodynamic behaviour of the whole human spine is beneficia l t o w h e e l c h a i r d e s i g n applications for the disabled. When applying forces to a certain vertebra of the spine under fixed constraints on sacrum and selected vertebrae, users such as surgeons or clinicians can feel force feedback from the spine as well as examine its locomotion. These results are useful for designing suitable and comfortable wheelchairs for the disabled with specific abnormal spinal configurations. In addition, by simulating in a haptically integrated graphic environment, orthopaedic surgeons can gain insight into the planning of surgery to correct severe scoliosis. Different rod and brace designs can be experimented with using this virtual environment. Furthermore, the surgeons may be able to understand the change in force distribution following spine fusion procedures, which can also assist in post-operative physiotherapy. The objective of this chapter is to present the development of a detailed human spine model with a haptic interface, which can be useful for investigating various medical applications.

Overview of human spine structure
The spinal column ( Figure 1) extends from the skull to pelvis and is made up of 33 individual vertebrae that are stacked on top of each other. The spinal column can be divided into 5 regions: 7 cervical (C1-C7), twelve thoracic (T1-T12), and five lumbar vertebrae (L1-L5) in the lower back, five bones (joined together in adults) to form the bony sacrum, and 3-5 bones fused together to form the coccyx or tailbone.
Intervertebral discs ( Figure 2) are soft tissue structures situated between each of the 24 cervical, thoracic, and lumbar vertebrae of the spine. Their functions are to separate consecutive vertebral bodies. Once the vertebrae are separated, angular motions in the sagittal (forward, backward bending) and coronal planes (sideway bending) can occur.
Facet joints are paired joints which are found in the posterior of the spinal column ( Figure  3). The surfaces of each joint are covered by cartilage which helps to smooth the movement between two vertebrae. Certain motions are facilitated by these joints, such as: bending forward, bending backward and twisting. In addition, people can feel pain if the joints are damaged because of the connected nerves. Some experts believe that these joints are the most common reasons for spinal discomfort and pain.
The neural elements ( Figure 4) consist of the spinal cord and nerve roots. The spinal cord runs from the base of the brain down through the cervical and thoracic spine. The spinal cord is surrounded by spinal fluid and by several layers of protective structures, including the dura mater, the strongest, outermost layer. At each vertebral level of the spine, there is a pair of nerve roots. These nerves go to supply particular parts of the body.   Figure  4), which run from the skull to the base of the spine (the sacrum). In addition there are also many muscles attached to the spine, which further help to keep it stable. The majority of the muscles are attached to the posterior elements of the spine.

Finite element modeling of the spine
In context of most researchers mainly focusing on modelling partially the spinal regions, the construction of a finite element model of the human spine can be useful for clinicians to investigate clinical problems by predicting its biomechanical behaviour. A beam finite element (FE) spine model for haptic interaction is built based on a solid FE spine model, which is created in offline finite element analysis (FEA) software. The mechanical properties of the beam FE spine model are tuned so that its deformation behavior is very similar to that of the offline solid spine model. Furthermore, the online beam FE spine model is greatly simplified as compared to the offline solid FEA model and hence more appropriate for realtime simulations. Haptic feedback is provided in the real-time simulation of the beam FE spine model, in order to enhance the human-computer interaction. Based on the results of spine deformation obtained from the haptic online FE simulator, the offline FEA spine model again is used to reproduce the same deformation and hence to provide more detailed deformation and vertebrae's stress/strain information, which the haptic beam FE model is not capable to provide. Figure 6 shows the architecture of the system. Effective haptic interaction requires a high updating rate of around 1kHz and only a simple FE spine model can be used without significant latency. FE solvers for haptic interfaces are often restricted to linear or simple non-linear solutions. Most research on FE simulation of the spine uses tetrahedral or brick element types in order to obtain thorough representations of spinal deformation. However, a haptic interface's requirement for high updating rate makes it difficult to use solid elements in the haptic virtual environment, because these often result in huge numbers of elements and nodes. Therefore, a beam element was chosen for this real-time haptic simulator. The beam element type provides very limited information concerning stress and strain of vertebrae or inter-vertebral discs, which is potentially important for representing realistic behaviour. In order to provide stress/strain information of spinal deformation, commercial FEA software (in this case ABAQUS) is employed using a tetrahedral element and computed offline. The complete simulation is achieved in two steps. Firstly, simulation in the online haptic simulator obtains quick, intuitive but relatively rough www.intechopen.com results. The displacement and deformation of the spine from the haptic simulator is then input to the offline FEA application for further detailed analysis. This offline computation is carried out when there is no direct haptic interaction and sufficient computational time.
Initially, two types of FE models, beam element and tetrahedral element models, were created based on 3D models of scanned vertebrae. Material properties of FE models should be validated before simulation. In the haptic simulator, a user can push or pull any vertebra in any arbitrary direction to observe the resulting deformation of the whole spine. After the user is satisfied with the deformation of spine model in haptic simulator, the displacements of some focused vertebrae are input to the offline commercial FEA application for more detailed simulation and the results of stress and strain information can be obtained.

Geometric modeling of vertebrae
In order to create a geometric model, a healthy spine is preferable to a specific diseased one. For our system, a resin spine prototype (Budget Vertebral Column CH-59X Life Size 29" Tall), which is cast from a Chinese-Singaporean cadaver, was digitized to create the geometric model of the spine (Figure 7(a)). A common method used to build computer models of the spine is by stacking computed tomography (CT) images sequentially to develop and discretize the 3D solid model. However, with this method, the computer model does not properly represent the surface geometry of the highly irregular vertebra, especially its posterior part. For this reason we used laser scanning of the individual surfaces of the model vertebrae to generate the initial point-cloud data set as can be seen in Figure 7(b). Spline (NURBS) surface models were then constructed based on the polygonal meshes generated from the scan data (Figure 7(c) and (d)). These NURBS models can then be used as templates for CAD model generation. Customized models can be obtained by modifying these templates according to dimensional specifications taken from individual patients (Mastmeyer et al., 2006). Once the geometric models of vertebrae have been prepared, FE models can be constructed according to these geometric models.

Finite element modeling for haptic rendering
A major concern is the high updating rate that is required for haptic rendering. Low rates can result in instability or vibration when using the haptic device. The FE model for haptic www.intechopen.com interaction is simplified and optimized so that the number of elements and nodes is kept as small as possible. In the haptic interface, spines are represented by beam elements. A similar strategy of simulating spines with beam elements can be found in (Lafage et al., 2004). Figure 8 shows the construction of a rigid beam element FE model of a vertebra. Points A and B are the center points of the upper and lower endplates of the vertebral body. Points C and D are the left and rightmost points of the transverse processes. Point E is the rearmost point of the spinal process. Vertebrae are represented using these rigid beams. Point A and B are then connected to intervertebral discs. Point C, D and E are to be connected to specific ligaments. Figure 9(a) and (b) show side and front views of a section of the assemblage of the spine finite element model and explanation is expressed in Figure 9(c). Our vertebrae are considered rigid and only the intervertebral discs (IVDs) and ligaments are deformable. The beam element of IVDs is of cylinder type. The flexural moments of inertia (I y , I z ) of beam elements for IVDs are calculated from the intervertebral discs' geometries. An analogy with the beam theory allows Young's modulus to be expressed as a function of the disc stiffness in bending (K f ), obtained from in-vitro experiments (Lafage et al., 2004). Similarly, the torsion moment of inertia was expressed as a function of the disc stiffness in torsion (K t ), obtained from in-vitro experiments (Lafage et al., 2004). Ligaments were modeled using tension-only truss elements. Two ligaments for connecting transverse processes and one ligament for connecting the spinous process were created. The sacrum is also included as a quasi-rigid body composed of a set of stiff beams. The material properties of ligaments and disks can be found by reference to (Lafage et al., 2004).

Finite element models for offline FEA
The purpose of the offline FEA is to gain information on the stress/strain relationships of vertebrae and IVDs. The FE model should be sufficiently delicate and precise in order to achieve this. Normal vertebral body models have a cancellous core covered by a cortical shell of approximately 1.5mm thickness. The thickness of most IVDs is around 10mm. The nucleus pulposus of the IVD is modeled as an incompressible material comprising 40% of the total disc volume. The models of the nucleus, the annulus and the vertebral bodies are assembled with tie constraints between the interacting surfaces. Besides vertebral bodies and IVDs, ligaments play an important role in spinal biomechanics and stability. The seven intervertebral ligaments were incorporated and modeled with truss elements. The ribcage is currently not included in this FE model.
Since two FE models are employed in this system, their deformation under the same boundary conditions should be coherent. Therefore, the mechanical properties of spine FE models must be carefully tuned in order to ensure deformations from both models match. A similar approach of material property adjustment can be found in (Lafage et al., 2004). Since the offline FEA is considered as more accurate than the online FE simulator, its results can serve as a standard, meaning the mechanical properties of the beam element model should be adjusted according to the behaviour of a solid FE spine model. The iterative process of mechanical personalization is engaged in the commercial FEA package. Firstly, the offline element spine model is applied with a force of 40N on the T1 vertebra. The force direction is from the rear to the front. The displacement of each vertebra is recorded and serves as the template. Then, the same load and boundary conditions are applied to the online beam element spine model. Differences between the beam element model and the template solid element model are quantified. Mechanical characteristics of soft tissues of the beam element spine are then tuned manually and locally until the difference between the two simulations tend to less than 5 degrees for the vertebral orientation and 10mm for the vertebral bodyline. The mechanical properties to be tuned here include: the stiffness of ligaments; the Poisson's ratio, the shear stiffness and the torsional stiffness of intervertebral discs.

The reproduction of haptic simulation result in the offline FEA simulator
Although the haptic real-time simulation of the beam FE spine model helps us to understand better the kinematics of the whole thoracolumbar spine, it provides very limited information concerning stress and strain of vertebrae. If a full picture of the deformation behavior of the spine is required, it is important to study and obtain the information of stress/strain of vertebrae. Therefore, offline simulations following the haptic simulation are necessary to carry out and provide this useful information. Our system allows the user to pick interesting moments of particular spine shape during haptic simulation and reproduce www.intechopen.com the similar deformation of spine in the offline FEA simulator for greater detail. By this means, the user can obtain desired spine deformations conveniently in the haptic simulator and then observe every detail of stress/strain of vertebrae in the offline FEA simulator.
A straightforward way to reproduce the same deformation in the offline FEA simulator is to record the position and orientation of each vertebra along the spine and then apply the positional constraints to every vertebra of the offline spine model. However, this greatly increases the risk of finite element computation error because of too many constraints involved. In practice, we found this method almost always makes the FEA solver fail to converge. Another plausible way is to record the force vector and the vertebrae under the load in the haptic simulator and reproduce the same load condition in the offline FEA simulator. However, the deformation of spine is not solely determined by the load due to the high non-linearity of the spine. Plus, during surgery, the forces are often not directly applied to a vertebra; instead, surgeons often manipulate instruments, which are linked to one or multiple vertebrae. Therefore, this method is not feasible.

Introduction to LifeMOD simulation software
Recently, many software applications have been developed for impact simulation, ergonomics, comfort study, biomechanical analysis and surgical planning. Such software enables users to perform human body modeling and interaction with the environment where the human motion and muscle forces can be simulated. These tools are very useful for simulating the human-machine behavior simultaneously. LifeMOD from Biomechanics Research Group is a leading simulation tool that has been designed for this purpose.
The LifeMOD Biomechanics Modeler is a plug-in module to the ADAMS (Automatic Dynamic Analysis of Mechanical Systems) physics engine, produced by MSC Software Corporation to perform multi-body analysis. It provides a default multi-body model of the skeletal system that can be modified by changing anthropometric sizes such as gender, age, height, weight etc. The created human body may be combined with any type of physical environment or system for full dynamic interaction. The results of the simulation are the human motion, internal forces exerted by soft tissues (muscles, ligaments, joints) and contact forces at the desired location of the human body. Full information on the LifeMOD Biomechanics Modeler can be found online (LifeMOD).
In following subsections, the developing process of a detailed musculo-skeletal spine model is presented thoroughly. This process includes five main stages such as generating a default human body model, discretizing the default spine segments, implementing ligamentous soft tissues, implementing lumbar back muscles and adding intra-abdominal pressure.

Generating a default human body model
The usual procedure of generating a human model is to create a set of body segments followed by redefining the fidelity of the individual segments. The body segments of a complete standard skeletal model are first generated by LifeMOD depending on the user's anthropometric input. The model used in this study was a median model with a height of 1.78 m and a weight of 70 kg created from the internal GeBod anthropometric database. By default, LifeMOD generates 19 body segments represented by ellipsoids. Then, some kinematic joints and muscles are generated for the human model. Figure 10 shows the base model in this study.

Discretizing the default spine segments
To achieve a more detailed spine model, the improvement of the default spine model mentioned above is required and can be done in the three following steps: refining the spine segments, reassigning muscle attachments and creating the spinal joints.
From the base human model, the segments may be broken down into individual bones for greater model fidelity. Every bone in the human body is included in the generated skeletal model as a shell model. To discretize the spine region, the standard ellipsoidal segments representing the cervical (C1-C7), thoracic (T1-T12) and lumbar (L1-L5) vertebral groups are firstly removed. Based on input such as center of mass location and orientation of each vertebra, the individual vertebra segment is then created. Figure 11 shows all ellipsoidal segments of 24 vertebrae in the cervical, thoracic and lumbar regions after discretizing.
The muscles are attached to the respective bones based on geometric landmarks on the bone graphics. With the new vertebra segments created, the muscle attachments to the original segment must be reassigned to be more specific to the newly created vertebra segments. The physical attachment locations will remain the same. Figure 12(a) and (b) shows the anterior and posterior view of several muscles in neck/trunk regions. Table 1 lists attachment locations of these muscles.  It is necessary to create individual non-standard joints representing intervertebral discs between newly created vertebrae. The spinal joints are modeled as torsional spring forces and the passive 3 DOFs jointed action can be defined with user-specified stiffness, damping, angular limits and limiting stiffness values. These joints are used in an inverse dynamics analysis to record the joint angulations while the model is being simulated. The properties of the joints can be found in the literature (Moroney et al., 1988;Panjabi et al., 1976;Schultz et al., 1979;Schultz & Ashton-Miller, 1991).

Creating the ligamentous soft tissues
To stabilize the spine model, interspinous, flaval, anterior longitudinal, posterior longitudinal and capsule ligaments are created. Figure 13 displays various types of ligaments attached to vertebrae in the cervical spine region

Implementing lumbar muscles
Multifidus muscle: The multifidus muscle is divided into 19 fascicles on each side according to descriptions by the group of Bogduk (Bogduk et al., 1992a;Macintosh & Bogduk, 1986). The multifidus can be modeled as three layers with the deepest layer having the shortest fibres and spanning one vertebra. The second layer spans over two vertebrae, while the third layer goes all the way from L1 and L2 to posterior superior iliac spine (Zee et al., 2007). The rather short span of the multifidus fascicles makes it possible to model them as line elements without via-points (Figure 14(a)).
Erector spinae muscle: According to (Macintosh & Bogduk, 1987;, there are four divisions of the erector spinae: longissimus thoracis pars lumborum, iliocostalis lumborum pars lumborum, longissimus thoracis pars thoracis and iliocostalis lumborum pars thoracis. The fascicles of the longissimus thoracis and iliocostalis lumborum pars lumborum originate from the transverse processes of the lumbar vertebrae and insert on the iliac crest close to the posterior superior iliac spine (Zee et al., 2007). The fascicles of the longissimus thoracis pars thoracis originate from the costae 1-12 close to the vertebrae and insert on the spinous process of L1 down to S4 and on the sacrum. The fascicles of the iliocostalis lumborum pars thoracis originate from the costae 5-12 and insert on the iliac crest. Since muscles of the two pars thoracis are automatically generated by L i f e M O D , o n l y m u s c l e s o f t h e t w o p a r s lumborum need to be added to the model as shown in Figure 14(b).

www.intechopen.com
Psoas major muscle: The psoas major muscle is divided into 11 fascicles according to different literature sources (Andersson et al., 1995;Bogduk et al., 1992b;Penning, 2000). The fascicles originate in a systematic way from the lumbar vertebral bodies and T12 and insert into the lesser trochanter minor of the femur with a via-point on the pelvis (iliopubic eminence) (Figure 14(c)). Bogduk found that the psoas major had no substantial role as a flexor or extensor of the lumbar spine, but rather that the psoas major exerted large compression and shear loading on the lumbar joints (Bogduk et al., 1992b). This implies that the moment arm for the flexion/extension direction is small and hence the via-points for the path were chosen in such a way that the muscle path ran close to the centre of rotation in the sagittal plane. Quadratus lumborum muscle: For modelling the quadratus lumborum, the description given by Stokes and Gardner-Morse was followed (Stokes & Gardner-Morse, 1999). They proposed to represent this muscle by five fascicles. The muscle originates from costa 12 and the anterior side of the spinous processes of the lumbar vertebrae and has in the model a common insertion on the iliac crest (Figure 14(d)).
Abdominal muscles: Two abdominal muscles are included in the model: obliquus externus and obliquus internus. Modelling of these muscles requires the definition of an artificial segment with a zero mass and inertia (Zee et al., 2007). This artificial segment mimics the function of the rectus sheath on which the abdominal muscles can attach (Figure 15(a)). The obliquus externus and internus are divided into 6 fascicles each (Stokes & Gardner-Morse, 1999). Two of the modeled fascicles of the obliquus externus run from the costae to the iliac crest on the pelvis, while the other four originate on the costae and insert into the artificial rectus sheath as can be seen in Figure 15(a). Three of the modelled fascicles of the obliquus internus run from the costae to the iliac crest, while the other three originate from the iliac crest and insert into the artificial rectus sheath (Figure 15(b)).

Adding intra-abdominal pressure
Since LifeMOD and ADAMS provide tools that only generate concentrated or distributed forces, it is not possible to implement directly intra-abdominal pressure into the spine model. To overcome this difficulty, a new approach to intra-abdominal pressure modeling is proposed. Initially, an equivalent spring structure able to mimic all mechanical properties of intra-abdominal pressure such as tension/compression, anterior/posterior shear, lateral shear, flexion/extension, lateral bending and torsion is created (Figure 16). After that, the translational and torsional stiffnesses of the spring structure are determined. Finally, since adding this spring structure into the spine model is quite troublesome, a bushing element that can specify all stiffness properties of the structure is used instead (Figure 17).

Validation of the detailed spine model
To validate the detailed spine model presented above, two approaches are used and presented as follows:  With the same extension moment generated in upright position, axial and shear forces in the L5-S1 disc calculated in the model are compared to those obtained from Zee's model (2007) and experimental data (McGill & Norman, 1987).  While a subject holds a crate of beer weighing 19.8 kg, the axial force of the L4-L5 disc is computed and compared with in-vivo intradiscal pressure measurements (Wilke et al., 2001).
In the first approach, a gradually increasing horizontal force was applied onto the vertebra T7 of the spine model from posterior to anterior in the sagittal plane. From this force, axial and shear forces as well as the moment about the L5-S1 disc were calculated. Zee's model estimated an axial force of 4520 N and shear force of 639 N in the L5-S1 disc at a maximum extension moment of 238 Nm. Meanwhile, to obtain the same extension moment, the external force that needs to be applied in the present model is 1260 N. Corresponding with this force, axial and shear forces obtained in the model were 4582 N and 625 N, respectively. This is in accordance with the results presented by McGill et al. (1987) who found axial forces in the range of 3929-4688 N and shear forces up to 650 N.
In the second approach, a comparison was made with in-vivo intradiscal pressure measurements of the L4-L5 disc as reported by Wilke et al. (2001). They measured a pressure of 1.8 MPa in the L4-L5 disc while the subject (body mass: 70 kg; body height: 1.74 m) was holding a full crate of beer (19.8 kg) 60 cm away from the chest. The disc area was 18 cm 2 and based on this the axial force was calculated to be 3240 N. The same situation was simulated using the spine model in this research. The estimated axial force was 3161.6 N. This is a good match considering the fact that no attempt was made to scale the model to the subject in this study. Body mass and body height of the subject in this study are quite similar to the body mass and height used in the model.

Haptic interface for spine modelling 6.1 Haptic rendering
Haptic rendering is the process of applying forces to give the operators a sense of touch and interaction with physical objects. Typically, a haptic rendering algorithm consists of two parts: collision detection and collision response. Figure 18 illustrates in detail the procedure of haptic rendering. Note the update rate of haptic rendering has to be maintained at around 1000 Hz for stable and smooth haptic interaction. Otherwise, virtual surfces feel softer. Even worse, the haptic device vibrates.
To observe the locomotion and study dynamic properties of the spine model quickly, conveniently and more realisticly, haptic technique can be integrated into a spine simulation system. In this study, for two types of spine models developed using finite element method as well as multi-body method in LifeMOD software, the haptic rendering process has two main stages: the rigid stage and the compliant stage. Without pressing the stylus button of the PHANToM device, the users can touch and explore the whole spine model since it is considered to be rigid throughout. After the users locate a specific vertebra where he/she wishes to apply force, they can then press the PHANToM stylus button and push or drag the vertebra in any direction to make the whole spine model deform. Once the stylus button is pressed, the system switches from rigid stage to compliant stage. The haptic rendering algorithms in these two stages will be clearly presented in the following sections. Figure 19 shows the complete haptic simulation process in the system.

Haptic rendering method for the rigid stage of the spine models
In real-time haptic simulation, users can only interact with the spine model by manipulating a rigid virtual object considered as a probe on the computer screen. At present, a simple probe such as a sphere is used in this study. In the rigid stage, a common haptic rendering method is used for the aforementioned types of spine models. Since this interaction carries out at a high update rate of 1 kHz, the chosen haptic rendering method needs to be reasonable and effectively computational. As presented above, the haptic rendering method includes two phases: collision detection and collision response. For the collision detection, the simple algorithm proposed by James Arvo (1990) is utilized to check intersection between the probe and the spine model through axis-aligned bounding boxes (AABBs). Figure 20 illustrates all AABBs of a vertebra intersecting with the spherical probe.
For collision response, the algorithm developed by Gao and Gibson (2006) is used to determine force feedback. An important prerequisite in this step is that intersecting points with the spherical probe have to be specified. Figure 21 shows intersecting points during colliding process between the probe and a vertebra.   ( 2) where n is the number of surface sampling points inside probe, N i as the normal vector of surface sampling point and s i as the depth.

Haptic rendering method for the compliant stage of the spine models
After the users locate a specific vertebra where he/she wishes to apply force, they can then press the PHANToM stylus button and push or drag the vertebra in any direction to make the whole spine model deform. Once the stylus button is pressed, the system switches to another haptic rendering algorithm that uses the stretched-spring model. A virtual spring is set up connecting the vertebra and the haptic cursor. The spring has two hook points: one is on the vertebra and the other is the haptic cursor. Both of the hook points displace during spine deformation. The force magnitude is determined by the length of the virtual spring and stiffness while the force direction depends on the vector of the virtual spring. www.intechopen.com

Finite element spine model
In this stage, the haptic interface has two components: haptic rendering and an FE solver. The function of the haptic rendering component is to provide haptic feedback to the user when the user virtually deforms the spine model. The FE solver evaluates the deformation of the FE spine model according to the user's force input.
To solve the FE deformation problem, an explicit Newmark integration method is employed. The theory of elasticity consists of a set of differential equations that describe the state of stress, strain and displacement of each point within an elastic deformable body. A relationship between the deformation of the object and the exerted forces can be established by synthesizing those equations. The equation governing the resulting deformation from external forces is: 33 nn  where u is the 3n-dimensional nodal displacement vector; and are the velocity and acceleration vectors respectively; F is the external force vector; M is the mass matrix; D is the damping matrix; R(u) represents the internal force vectors due to deformation; n is the number of nodes in the FE model.
If deformation is assumed to be small, a linear elastic model can be used to approximate linear strain: where x, y and z are the independent variables of the Cartesian frame, and u, v and w are the corresponding displacement variables at the given point.
Because the internal force vector is linear with respect to the nodal displacement vector, this leads to a constant stiffness matrix, which can be pre-computed. However, this solution does not consider geometric non-linearity and is only fit therefore for linear and small deformation situations. The linear strain model introduces distortions for large global deformations. In the model the deformation of the spine is rather large and therefore the linear strain approximation is not suitable. We model the deformation using the Green strain as follows: Zhuang's mass concentration method is employed to approximate the distributed mass in real-time (Zhuang & Canny, 2000). To avoid inverting a large sparse matrix at each time step, Zhuang's method approximates the force matrix M with a diagonal matrix. Diagonalization results in approximating the mass continuum as concentrated masses at each nodal point of the finite element model. This simplifies the nonlinear system of equations used to form a set of independent algebraic equations that can be solved with relatively high updating rate.
A prototype system has been implemented. Although the FE model of the spine is a simplified beam model, the updating rate of FEM evaluation is still far below the 1000 Hz required for effective vibration-free haptic rendering. The FE model has 472 elements and 473 nodes and the time step of the explicit integration is 0.02 second. The updating rate is approximately 200 Hz using a PC graphic workstation. To avoid vibration due to the updating rate discrepancy, a force interpolation method (Zhuang & Canny, 2000) was adopted to smoothen the force feedback. Haptic feedback between two simulated states is linearly interpolated at a suitably high frequency. However, at time t n , we do not have the information of t n+1 . Zhuang's solution uses an intentional delay of an FEM deformation time step, up to 1/100s. User experimentation found that such a small lag in time is within the tolerance of human perception for virtual interaction with soft objects.

Multi-body spine model
Different from the finite element spine model whose locomotion is directly calculated in the haptic simulation, the multi-body spine model in the compliant stage utilize the preinterpolated displacement-force functions of all vertebrae to conveniently compute the movement of the spine. Based on the magnitude of the haptic external forces applied by the users, dynamic properties of all vertebrae can be easily calculated via these functions and the locomotion of the whole spine model can be rapidly observed.
The purpose of interpolating displacement-force functions of all vertebrae is to facilitate computational process of dynamic properties of the spine model during haptic simulation. To do that, some constraints are imposed on the spine model. The pelvis is fixed in 3D space. Then, constant forces are applied on each specific vertebra in the thoracic region in each axisaligned direction during simulation. The force magnitude is gradually increased with an equal increment in subsequent simulations. Correspo n d i n g t o e a c h v a l u e o f f o r c e , d y n a m i c characteristics of all vertebrae (e.g. translation, rotation) can be automatically obtained using the plots in LifeMOD as a reference. The dynamic properties are recorded after the spine model is stabilized. Based on these recorded dynamic properties, displacement-force www.intechopen.com relationships can be interpolated using the least-squares method and expressed in terms of polynomial functions. Figure 22 illustrates the graph of a dynamic property of vertebra T1 under forward forces in sagittal plane of the spine model. In these figures, each series of markers presents each dynamic property of one vertebra obtained in the simulations and the continuous line closely fit to that series is the corresponding interpolated line. Initially, without pressing the stylus button, the user can touch and explore the spine model since it is considered to be rigid throughout. Then, the user can control the probe to apply external forces in any arbitrary direction onto any vertebra by pressing the stylus button to observe the locomotion and study dynamic properties of the spine model. Figure 23 shows the haptic real-time simulation of the finite element spine model under lateral bending. During this real-time simulation, the displacement, rotation and deformation of the spine can be recorded and then input to the offline FEA application for further detailed stress/strain and deformation behaviour analysis of IVDs as well as vertebrae.

Studying stress/strain of the finite element spine model in a prone position
As presented above, all dynamic properties of the spine model recorded during the haptic simulation can be provided as input of offline FEA simulator for further detailed analysis. In this section, an example of spine manipulation is demonstrated. A human body is prone on a horizontal table and a downward force is pushed on T12 vertebra. Therefore, the spine is also in prone position. All 6 DOFs of the pelvis are constrained and 2 DOFs of the T1 on the horizontal plane are constrained in order to simulate the situation where the prone posture torso is constrained by a horizontal plane underneath the torso. In the haptic simulator, downward pressing forces are applied to the T12 vertebra by the user through the PHANToM device. Figure 24(a) shows the beam element spine model at initial status without any deformation.  Figure 24(b) shows the spine model under deformation. The deformation is then recorded and the displacement of the focused vertebra, T12 in this case, is obtained and input into the offline FEA solver. Figure 25 shows the result from the ABAQUS offline FEA solution. Only T9 to the pelvis is shown because this region is the most interesting region if a pressing load is applied on T12. The offline FEM model consisted of cortical bones, cancellous bones, bony end plates and the intervertebral discs. The disc's annulus was represented as a composite of collagenous fiber embedded in a matrix of base substance. The nucleus was modelled as an incompressible inviscid fluid. Because intervertebral discs often experience large deformation and strains under compressive loading, the materials of the FEM spine model were non-linear. Figure 25(a) shows the tetrahedral FE model without deformation. Figure  25

Investigating dynamic properties of the spine model under external haptic forces
Since the finite element spine model for the haptic simulation is represented with beam element which is in fact highly simplified, it is difficult to describe precisely the dynamic properties of the real spine. Compared to this finite element spine model, the detailed musculo-skeletal multi-body spine model developed in LifeMOD can provide better information of biodynamic behaviour of the spine. Since LifeMOD software takes account of other components (such as head, ribcage, muscles, ligaments so on), the locomotion of the spine model will become much more realistic. Furthermore, by simulating the spine model in a haptically integrated graphic environment, the users (such as surgeons, trainers) can quickly and conveniently observe the locomotion as well as gain insight into the dynamic behaviour of the spine.

www.intechopen.com
The long-term goal of this study is to develop a haptically integrated simulation platform for investigating post-operative personal spine models constructed from LifeMOD. The first step in this study is to develop a haptic interface for a detailed normal spine model whose dynamic properties are obtained from LifeMOD. Figure 26 illustrates the analysis of some dynamic properties of the spine when applying force on vertebra T1 in x-axis direction.
The next step in the study is to develop a post-operative detailed spine model. The operation conducted on the spine model can be spinal fusion or spinal arthroplasty. By using this haptically integrated graphic interface, assessing biomechanical behavior between natural spines and spinal arthroplasty or spinal fusion becomes more convenient. This will help orthepaedic surgeons understand the change in force distribution following spine fusion procedures, which can also assist in post-operative physiotherapy. In addition, the surgeons can gain useful insight from this system in the planning of surgery to correct severe scoliosis. Different designs of rods and braces can for example be experimented with using this virtual environment. It is also considered that this system may assist physiotherapists and other practitioners when designing cushioning and supports for wheelchair-bound and other clients suffering from asymmetric spinal disorders.

Conclusion
In general, we have presented the process of modelling human spine. Initially, a beam FE spine model for haptic interaction is built based on a solid FE spine model, which is created in offline finite element analysis software. Although there are some limitations on not fully taking into account soft tissues, this beam FE model is able to simulate global dynamic behaviour of the whole spine, Subsequently, an entirely detailed musculo-skeletal mutibody spine model using LifeMOD Biomechanics Modeler is completely developed to solve the limitations aforementioned and reliably validated by comparing simulation results with experimental data and in-vivo measurements. Moreover, to enhance more realistic interaction between users (such as trainers, clinicians, surgeons) and the spine model during real-time simulation, haptic technique is integrated into a graphic interface. Based on this, exploration of the spine model becomes much more realistic since the users can manipulate the haptic cursor to directly touch, grasp and even apply external forces in any arbitrary direction onto any vertebra to observe the deformation of the spine. In such a simulation interface, users can quickly and more conveniently study the locomotion and dynamic behaviour of the spine model. There has been significant progress in haptic technologies but the incorporation of haptics into virtual environments is still in its infancy. A wide range of the new society's human activities including communication, education, art, entertainment, commerce and science would forever change if we learned how to capture, manipulate and reproduce haptic sensory stimuli that are nearly indistinguishable from reality. For the field to move forward, many commercial and technological barriers need to be overcome. By rendering how objects feel through haptic technology, we communicate information that might reflect a desire to speak a physicallybased language that has never been explored before. Due to constant improvement in haptics technology and increasing levels of research into and development of haptics-related algorithms, protocols and devices, there is a belief that haptics technology has a promising future.