UFDC Home  myUFDC Home  Help 



Full Text  
A TWODIMENSIONAL HUMAN SPINE SIMULATION AND THREE DIMENSIONAL SPINE MODEL CONSTRUCTION By JIANZHI LIU A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2001 ACKNOWLEDGMENTS Many thanks are in order for Dr. Andrew J. Kurdila for his continual assistance, expertise, and leadership. I would like to thank Dr. Andrew J. Rapoff and Dr. B. J. Fregly for their assistance and expertise in helping me on this project. I would also like to thank all the faculty and students who helped at different times during my research. TABLE OF CONTENTS A C K N O W L E D G M E N T S ................................................................................................... ii LIST OF TABLES ............................... .................................. LIST OF FIGURES ..................................... ................ ............... ... ..... .. vii A B STR A C T .................................................... ........ ........... ................. .. x CHAPTERS 1 INTRODUCTION AND LITERATURE REVIEW....................................................1 1.1 In tro du ctio n ...................................... ................................. ............... 1 1.2 Literature Review .......................... ........... ..... 3 1.3 G o al an d O utlin e ............................................ ........................... 6 2 MULTIBODY DYNAMICS THEORY....... ...... ..............................................8 2.1 TwoDimensional Formulation Example................ ..... ............. 8 2.2 Identification Problem in the TwoDimensional Example .................................... 15 3 NUM ERICAL EXPERIM ENT ............................................................ ............. .22 3.1 Multibody Dynamics Simulation ToolsPromechanica..................................... 22 3.2 Forward Simulation with Custom Load.......................................................... 24 3.3 Identification of Pin Joint Position Using Inverse Simulation............................ 36 4 SCOLIOSIS MODEL DATABASE CONSTRUCTION ..............................................48 4 .1 G eom etry ..................................................... 4 8 4.2 Stiffness............................. ....... ............................... 58 4.3 A Simple Demonstration of the Stiffness Database...................... .............. 75 5 CONCLUSIONS AND RECOMMENDATIONS ............................... ................77 APPENDICES A GEOMETRY OF MODELING NODES .....................................................................78 B ATTACHING POSITIONS OF SCREW & HOOKS.................................................95 L IST O F R E FE R E N C E S ...................................................................... ..................... 112 BIOGRAPHICAL SKETCH ........................................................................114 iv LIST OF TABLES Table Page 31 8 Springs identification........................................................................... ....................34 32 16 Springs identification.......................................................................... ...................34 41 Scoliosis curve data ................... .................................... .. 50 42 Coordinates of origin nodes............................................... .........................................51 A l Lum bar 5 ................................................................78 A 2 Lum bar 4 .......... ...................................................................79 A 3 Lum bar 3 ................................................................80 A 4 Lum bar 2 ........................................................................ 81 A 5 Lum bar 1 ..................................................................................82 A 6 Thoracic 12 ................................................................83 A 7 Thoracic 11 ................................................................84 A 8 Thoracic 10 ................................................................85 A 9 Thoracic 9 ................................................................86 A 10 Thoracic 8 ................................................................87 A l Thoracic 7 ................................................................88 A 12 Thoracic 6 ................................................................89 A 13 Thoracic 5 ................................................................90 A 14 Thoracic 4 ..............................................................................91 V A 15 T h oracic 3 ................................................................92 A 16 T h oracic 2 ................................................................9 3 A 17 T h oracic 1 ................................................................94 B l Screw & H ooks on L 5 (m m ) ................................................. ............................... 95 B2 Screw & H ooks on L4 (m m ) ................................................. ............................... 96 B 3 Screw & H ooks on L3 (m m ) ................................................. ............................... 97 B4 Screw & H ooks on L2 (m m ) ................................................. ............................... 98 B 5 Screw & H ooks on L (m m ) ................................................. ............................... 99 B 6 Screw & H ooks on T 12 (m m ) ................................................ .............................. 100 B 7 Screw & H ooks on T (m m ) .............................. ... ........................................101 B 8 Screw & H ooks on T 10 (m m ) ................................................ .............................. 102 B 9 Screw & H ooks on T9 (m m ) ................................................. ............................... 103 B 10 Screw & H ooks on T8 (m m ) ................................................ .............................. 104 B 11 Screw & H ooks on T7 (m m ) ........................................ ...................................... 105 B 12 Screw & H ooks on T6 (m m ) ................................................ .............................. 106 B 13 Screw & H ooks on T5 (m m ) ................................................ .............................. 107 B 14 Screw & H ooks on T4 (m m ) ................................................ .............................. 108 B 15 Screw & H ooks on T3 (m m ) ................................................ .............................. 109 B 16 Screw & H ooks on T2 (m m ) ............................................... .................................... 110 B 17 Screw & H ooks on T (m m ) .................................... ............................................ LIST OF FIGURES Figure Page 11 H um an Spine .................................................. 1 21 Tw o vertebral body system ......... .................................... .................... ............... 9 22 Tw obody exam ple illustration.......................... ................................... ............... 11 23 Identification procedure illustration.................................................... .. ................. 17 24 N um erical sim ulation procedure illustration ........................................ .....................21 31 Human spine modeling procedure .............. ..................................... 24 32 A two vertebral body model built in Promechanica.............. .......................... 25 33 Generalized forces in twodimensional plane.................................. .......... ......... 26 34 Ytranslational trajectory numerical results from Promechanica.............. ...............28 35 Ytranslational trajectory numerical results from Matlab model............... .................28 36 Ztranslational trajectory numerical results from Promechanica...............................29 37 Ztranslational trajectory numerical results from Matlab model................. ............29 38 Xrotational trajectory numerical results from Promechanica................ .. ...............30 39 Xrotational trajectory numerical results from Matlab model............... ................ 30 310 Ytranslational trajectory from Promechanica model ...............................................31 311 Ytranslational trajectory from Matlab model ....................................... ...........31 312 Ztranslational trajectory from Promechanica model............................... ...............32 313 Ztranslational trajectory from M atlab model .................................... ............... 32 314 Xrotational trajectory from Promechanica model........................................................33 315 Xrotational trajectory from Matlab model ....................................... ............... 33 316 Translation m ovem ent sim ulation....................................................... ............... 35 317 R otation m ovem ent sim ulation............................................. ................................... 35 318 vertebral body connected by pin......................................................... ............... 36 319 Probability density distribution of spring stiffness.....................................................38 320 Cumulative distribution of spring stiffness......................................... ............... 38 321 Probability density distribution of spring stiffness.....................................................39 322 Cumulative distribution of spring stiffness......................................... ............... 39 323 Probability density distribution of spring stiffness...................................................40 324 Cumulative distribution of spring stiffness......................................... ............... 40 325 Probability density distribution of spring stiffness.....................................................41 326 Cumulative distribution of spring stiffness......................................... ............... 41 327 Probability density distribution of spring stiffness...................................................42 328 Cumulative distribution of spring stiffness......................................... ............... 42 329 Probability density distribution of spring stiffness...................................................43 330 Cumulative distribution of spring stiffness......................................... ............... 43 331 Probability density distribution of spring stiffness...................................................44 332 Cumulative distribution of spring stiffness......................................... ............... 44 333 Probability density distribution of spring stiffness...................................................45 334 Cumulative distribution of spring stiffness......................................... ............... 45 335 Probabilty density distribution of different spring numbers........................................46 336 The cumulative distribution of different spring numbers............................................46 337 Error vs. different num ber of springs.................................................. ..... .......... 47 41 Spine M odel in Prom echanica M otion........................................................ ... ........... 52 42 A Vertebral Model with Complex Surface...................................................................53 43 A Sim ple Geom etry V ertebral M odel........................................ .......................... 53 44 Right Lateral View of the Primary Nodes ........................................... ............... 54 45 Top View of the Prim ary N odes ...................................................................... 55 46 Bottom V iew of the Prim ary N odes ................................... ........................................... 56 47 Sim plified G eom etry Spine M odel ...................................................................... ...... 57 48 Illustration Stiffness Property in Threedimensional Coordinates System....................58 49 Flextionextension between L1 and L2...................................... ......................... 60 410 Flextionextension between L2 and L3................................ ................................. 61 411 Flextionextension betw een L3 and L4...................................... ....................... 62 412 Flextionextension between L4 and L5................................................... ............... 63 413 Flextionextension betw een L5 and S1...................................... ........ ............... 64 414 Right and Left Axial Torque between L1 and L2 ....................................................65 415 Right and Left Axial Torque between L2 and L3 ............... ............................... 66 416 Right and Left Axial Torque between L3 and L4 ....................................................67 417 Right and Left Axial Torque between L4 and L5 ............... ............................... 68 418 Right and Left Axial Torque between L5 and S ................................................69 419 Left and Right Bending Moment between L1 and L2...............................................70 420 Left and Right Bending Moment between L2 and L3 ...........................................71 421 Left and Right Bending Moment between L3 and L4...............................................72 422 Left and Right Bending Moment between L4 and L5 ...........................................73 423 Left and Right Bending Moment between L5 and S .............................................74 424 LoadDisplacement Compare between Database and Reconstruction .......................76 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science A TWODIMENSIONAL HUMAN SPINE SIMULATION AND THREE DIMENSIONAL SPINE MODEL CONSTRUCTION By Jianzhi Liu December 2001 Chairman: Dr. Andrew Kurdila Major Department: Aerospace Engineering, Mechanics, and Engineering Science This thesis describes two different, but complementary modeling tasks for the human spine. Methodologies for modeling and simulation of the motion of the human spine vary dramatically in complexity. While providing detailed stress and strain field representations, full finite element modeling of the human spine can be computationally expensive. Alternatively, nonlinear multibody dynamics representations are often used because of their simplicity. These formulations employ rigid models of vertebrae interconnected via conventional mechanical joints. However, it is well documented that intervertebral motion can depart significantly from conventional mechanical joint constraints. We present an identification methodology that employs a relaxation technique in which joint mechanical properties are represented via a probability measure. In addition we describe an effort to obtain accurate threedimensional models of the human spine. These efforts are motivated by considering scoliosis. Scoliosis is defined as abnormal lateral curvature of the spine. It is usually considered as a three dimensional deformity, because axial rotation will always accompany the lateral curvature. The correction of the deformity is required when the patient risks severe deformity. A threedimensional spine model is constructed and a computer simulation tool is provided. CHAPTER INTRODUCTION AND LITERATURE REVIEW Chapter 1 first talks about some fundamental concepts and issues related to the human spine modeling, then discusses some basic information about scoliosis. An outline of the following chapter is given last. 1.1 Introduction Because each spine is unique and reacts in a slightly different way, the dynamics of the human spine are difficult to quantify. Therefore, in cervical  order to study the dynamics of the spine, it is necessary to vertebrae . cervical lordosis develop a mathematical model of the human spine that is Y general enough to encompass the possible variations in , thoracic 6 form and function, but specific enough to react realistically vertebrae thoracic kviphosis under loading and stress. kvhosis Human vertebral column (Figure 1.1) is a complex structure and mechanism. It is made up of bony vertebrae lumbar lordosis interconnected by discs and ligaments. The normal spine consists of seven cervical, twelve thoracic, five lumbar, sacrum five sacral and five coccygeal vertebrae. The sacral and coccyx coccyx coccygeal vertebrae are usually fused. In the lateral plane, the normal spine exhibits physiologic cervical lordosis, Figure 11 Human Spine thoracic kyphosis, and lumbar lordosis. They increase the flexibility and shock absorbing capacity of the vertebral column, while maintaining the adequate stiffness and stability. Although a vertebra is not actually a rigid body, usually it can be considered to behave as such, because deformations within the vertebra are insignificant compared to the deformations of the connective tissue. Two adjacent vertebrae, the ligaments, inter vertebral disc and the facets form a functional spinal unit, or a motion segment. It exhibits biomechanical characteristics similar to those of the entire spine. Each of the motion segments has six degreeoffreedom (DOF), and the whole spine can produce three movements: flexion extension, lateral bending, and axial rotation. We simply summarize them as the combinations of rotational and translational movements. The motions in human spine are coupled, where translation or rotation along or about one single axis is consistently associated with translation or rotation along or about a second axis. Some literature has shown these coupled motions clearly [Oxa92], [Pan94]. The coupled motion can lead to computation difficulties while modeling a complex structure like the human spine. Another fundamental issue related to our research is scoliosis. Scoliosis is defined as abnormal lateral curvature of the spine. It is usually considered as a three dimensional deformity, because axial rotation will always accompany the lateral curvature. Most scoliotic curves will stabilize if not treated, but a few will progress. Mild scoliotic curves result in cosmetic deformities of the trunk. These can create psychological problems for the patients who have them. Severe scoliotic curves may alter balance and coordination, interfere with the function of internal organs, allow premature degeneration of the spinal column, and cause deterioration of neurological function. Thus, the correction of the deformity is required when the patient risks severe deformity. The purpose of this correction is to stabilize the spine until a solid fusion occurs in the treated area. Once the spine has been fused, it can once again perform its function without the aid of the instruments. Hence, it is necessary to generate a general model that can represent these problems related to the human spine. Many efforts haven been made in the pass. They can be categorized into different sub areas according to the researchers' interests. 1.2 Literature Review Multiple biomechanical models have been employed for the spinal modeling. Most can be grouped into two classes. These classes are computational models and experimental models. The computational models can be used to provide insight into the mechanisms by which a particular motion occurs or can be used to investigate the complex dynamics of multilinked segments during motion. Also, these models are used to undertake analyses of activities, which would otherwise be impossible to carry out experimentally. Computational modeling can be classified into two categories: one includes the numerical models, which usually use finite element methods. The other includes the analytical models, where the geometry is simplified and requires few boundary conditions. Both categories of models have their strengths and weaknesses, but the inherent simplicity of analytical models often leads to an intuitive understanding of human spine mechanisms. Not surprisingly, the computational cost and accuracy of these methods vary greatly. Detailed finite element discretisations of intervertebral discs can provide fine scale stress and strain field approximations. Such highly resolved fields are obtained, of course, at significant computational cost [Bel73] [Bel74] [Kul75] [Kul76] [Sch73]. A recent study reports on the development and validation of a nonlinear viscoelastic FE model that qualifies the mechanical responses of the L2/L3 motion segment subset to time varying external loads. The vertebral bodies are represented as rigid bodies, and intervertebral discs are modeled via nonlinear finite elements. An alternative choice makes use of nonlinear multibody dynamics formulations [Buc96] [Hur94]. In these formulations, the vertebral bodies are assumed to be rigid and interconnected by a collection of conventional mechanical joints. The advantages of this formulation include the fact that these models are relatively low dimensional. It can often provide accurate predictions of large amplitude articulated motion. Optimization techniques are usually used in the nonlinear dynamics formulation, where objective functions are chosen, and a set of parameters is found that minimizes that function. One earlier example is the model developed by Schultz [Sch73]. In his study, Schultz described his two objective functions, which minimized the compression on the motion segment. In addition to computational modeling that determines the physics of spinal motion, a significant number of clinical experimental studies of spinal motion have appeared. Some of them are relatively old and represent early efforts to validate specific spine models [Kul75]. Other earlier experimental studies have estimated effective compliance parameters for motion segments. These are subsequently employed to provide more accurate lumped simulation models [Wil83]. There are some other experimental studies of spinal motion that have been carried out independent of numerical simulation, solely to examine and characterize the physics of spinal motion. Panjabi et al [Oxa91] [Pan94] investigated a threedimensional human spine structure that exhibits nonlinear viscoelastic mechanical behavior. As advanced measurement technologies have become more and more common, the experimental techniques and protocols have likewise matured. A hybrid approach called EMG assisted optimization (EMGAO) was developed by using both optimization and an EMG assisted approach. In other studies [Har98] [Ros98], lateral radiographs are employed to characterize spinal kinematics and frequency domain characterization of the spinal motion is studied. Despite the progress made in computational modeling and experimental methodologies to refine models, significant technical issue remains unresolved. As noted by M. Klenberger in [Kle98] "... Biological material data remains one of the limiting factors for biomechanical modeling. Existing data are comprised mostly of quasistatic forcedisplacement measurements with calculations of linear stiffness. Normal variance of biological tissues typically yields a standard deviation for stiffness values of at least 20 percent. No matter how refined an element mesh might be, the accuracy of the model can never be better than 20 percent." Thus, the uncertainty in data derived from biological specimens has lead to increased interest in identification and inverse methods, because identification and inverse methods by their nature are developed to accommodate some level of noise uncertainty and "poorly conditioned" data. Research [Wol94] [Ris97] [Kle98] [Xu99] provides recent examples in which identification and inverse methods have been developed in biomechanics. 1.3 Goal and Outline The objective of our research is to derive an identification method for human spinal motion that is based on a relaxation of the original governing equations. This approach is based on a simple observation. While some human joints are well represented by simple conventional mechanical joints, the relative motion of human vertebrae is not amenable to such a straightforward representation. In the following chapters, we present an identification methodology that employs a relaxation technique in which joint mechanical properties are represented via the probability measure. Numerical studies are used to study the performance of this approximation framework. Chapter 2 discusses the fundamentals of the multibody dynamics formulation. The study of the dynamics of the multibody system includes a development of the governing equations of the motion. In the model, each spine motion segment has six degree of freedom (DOF), and the whole human spine can undergo the translational and rotational movements. Based on this observation, the governing equations of motion (EOM) of the human spine are developed by using Lagrange's Equation. Based on the identification theory and algorithm introduced in chapter 2, the numerical experiments are implemented in chapter 3. In this chapter, several cases of identification are studied. A forward numerical simulation and an inverse identification numerical experiment are discussed. Most of the simulation work has been achieved using Promechanica, which is a multibody dynamics simulation software. In chapter 4, some basic tools for scoliosis research issues have been developed. These tools include a spine (L5T1) model constructed using the simulation software Promechanica in chapter 4. To facilitate the future development of spinal modes incorporating rod, hook and screw assembles, a spinal simulation database was developed that includes both geometry and attachment parts. The attachment points of pedicle screw and hooks are also included in this geometry database. In addition, a stiffness database has been constructed from the published literature. At the end of this chapter, a numerical simulation experiment is studied to verify the stiffness properties. The detailed geometry database and stiffness database can be found in the Appendix. Chapter 5 provides concluding remarks and suggestions for future work. CHAPTER 2 MULTIBODY DYNAMICS THEORY In the first part of chapter 2, a quick review of fundamental multibody dynamics theory is given. The governing equations for the model are derived based on Lagrange's Equations. In the second part, we discuss the basic procedures of numerical simulation. 2.1 TwoDimensional Formulation Example Before talking about several aspects of the theoretical background for the relaxed formulation used in the multibody identification problem, we derive the equations of motion. In this section we first show that classical Lagrangian methods of dynamics are reasonable to model the dynamics of the human spine. Generally, we can describe the configuration of a rigid body in space by using six coordinates. Three of them describe the body translation; three other coordinates define the orientation of the body. In twodimensional analysis, the configuration of a rigid body system can be identified by using three coordinates; two of them define the translation of the points on the body, and one coordinate defines the orientation of the body with respect to a selected inertial frame of reference. Consider one rigid body in figure 21. From the viewpoint of rigid body kinematics, the global position of arbitrary point Pi on the rigid body I can be defined as r' = 1R+y (21) where r' is the global position of point Pi, Rc is the global position vector of the origin Oi Figure 21 Two vertebral body system of the body I in the fixed reference frame Bi, and / is the position vector of point Piwith respect to Oi. So the whole body configuration can be completely defined, provide the vector components of R,, y' in the right side of the equation are known. The vector r' and R, define in the global coordinates, and / is defined in the body fixed local coordinates. We need to express the vector y in terms of the global systems. This is done by defining the orientation transformation of the body fixed coordinates with respect to the global frame of reference. For example, the model depicted in figure 21 is comprised of two vertebral bodies, denoted as I and J. Each body has an associated body fixed reference frame Bi and Bj, whose orientation is defined via a change of basis expressed in terms of an orientation matrix [L(,0)] as follows. b, I = [L(0,) b*, I S b ,,2 [L(O j2 (22) b = [L 0i k Le3 In this equation, bk is the kth basis vector defining body i, e~ is a basis of the inertial frame E, and O, is the angle measures the rotation of the body frame relative to the inertial frame. The location of the center of mass of body Bi and Bj is given by R 1e + (23) R = Xe+ Y,2 (24) R c=X +Y = (24) respectively. Similar expression for velocity and acceleration of the arbitrary point on the rigid body can be obtained. Thus the position, velocity, and acceleration of an arbitrary point on the rigid body can be written in terms of these coordinates. Equation (21) implies that the general motion of a rigid body is equivalent to a translation of one point, usually the original point of the body fixed frame, plus a rotation. From the viewpoint of kinetics, we can describe the kinetic energy of a rigid body as the sum of translational kinetic energy and the rotational energy. The total kinetic energy of a system is the sum of its components' kinetic energy. In the absence of any explicit constraint between the two bodies, the kinetic energy of the system is easily shown to be 1 2 .2 .2 .2 .2 .2 T=m X +YY +Jr6 +m X +Y +J 6 2 2 2 j J J+ 2 J J X m * m, (25) 1 2 2 2 2 2 J 0 = X Y,0 X0 YJ Yi where m, is the mass of body I and J, is the inertia of body I about the 3axis passing through its center of mass. X, Y and 0 are the translational and rotational velocity of the body I. They are all defined in the global coordinates system. The same description can be applied to the other body J. e2 f \P ' ."Ri Figure 22 Twobody example illustration If we sought to enforce a conventional revolute joint at the point denoted as "s" in figure (22), and assume it is frictionless, so no work is done by this constraint, the kinematics constraint could be expressed as D(s) = P(s) RP (s) X, + [1(0, )TY 1(S) X [(O y,,, (s) (26) In this 2D and twobody system, y,,1 and yY, are defined as the translational position of the arbitrary points Pi and Pj in X direction, and Y,2 and y,,2 as the translational position of the arbitrary points Pi and Pj in the Y direction. Both of them are the functions depending on the position of the revolute joint "s". One method to solve this problem is to use Lagrange multiplier equation. A classical Lagrange multiplier formulation of the equations would be d IT IT T V m  +y '2 = Qk k = ...N (27) dt 8 qk aqk I=1 ak where Qk are the generalized forces, q, are the generalized coordinates. By defining potential energy V, = I (s)c(s)s(ds) (28) 2 instead of using Lagrange multiplier, we could use the standard Lagrange's equation to solve an unconstrained problem. d 8T 8T 8V 8V, T  + + =Q Qk k=l1.N (29) dt qk) c/ c O qk where yu is the probability measure. In this case, the kinetic energy is the same as that in the classical Lagrange's multiplier formulation. We just need to calculate the potential energy part SY 1 T d (s)(s )j(ds)2 q 2 (210) = (f (T (s)) (s)u(ds) + (s) (S(ds)) 2 Oq 8q Alternatively, this equation can be written in terms of the summation convention as IVY I a D, (s)2l (s)t(ds) 8q, 2 aq, ( 1 Lj0 (s) BO(s) f (s)p(ds) + (s) (ds) 2 F8q, O8q V ( (0, (s) = (s)u(ds) (211) Oq, 8q, In this form, it would be very difficult, if not impossible, to solve the above relaxed formulation. The probability measure / in equation (210), (211) could be an element of an infinite dimensional space. To make the problem tractable, we employ a simple approximation of the probability measure yu. N /(.)= ZK1 () (212) l=1 In this equation, s is the dirac measure concentrated at s,. That is S1 s, EA s, (A) = (213) From a mathematical standpoint, we approximate p by /pk which is the weighted sum of dirac functions concentrated at sk,k = 1.. m. From a practical standpoint, this approximation can be thought as seeking a distribution of spring elements along the curve between the two rigid bodies. The expression in equation (211) can be simplified to aV 5K (SI) 8q, Z=1 aq, rT (214) Recall the constraint function is given as x(s)R_ (s) R_ (s) +X, y[(O ,{ (s){ X [ ]T (s) (215) Y, Y,,2 (S)j [ ) Yj, 2 (S) Recall orientation transformation of the basis, the transform matrix [L(O)] can be expressed in the two dimensional plane as follows, &b,i cosO, sin6, 0 ^1 Lbo,2 sin 0 cosO, O 02J (216) b,,3 0 0 1 le3 rel\ cos0, sin, 0 b,,, e = sin 0, cos 0H bb (217) e3 0 0 1 b],3 Plug the basis transform equations (216) (217) into the constraint equation, the constraint vector D becomes X, +cosO, sinO, l ,(s) X ~rX _cosO sin0 Jy,(s) D(S)= Lsine, cosO, [Z,2(S)J l Ls in, cos, jy,2(S) (218) So, we can describe the constraint function in the inertial frame E, and then calculate the Jacobian in equation (214) by definition, Oqj a0, (s) ax, OY, a0, aQ x ayj a0 q 02 002 002 02 02 2(219) ax, aY, ao, aOx aY a0o where the generalized coordinates q are defined as q = {X, Yo, X, YJ 8, (220) They are the translational and rotational basis in a twodimensional coordinates system as we mentioned above. If we plug the constraint function into equation (219), differentiate the vector components separately, we get the Jacobian matrix S 0 (sin(60)y1cos(6)y2) 1 0 (sin(8 )y, +cos(8 )y2 ) qJ LO 1 (cos(O,)y, sin(O,)y,2) 0 1 (cos(O,)y,, +sin(8 )Y,2)y (221) At this point, we have derived both the kinetic energy and the potential energy for this twobody system using Lagrange's Equations d cT a T +V d O + OQk (222) dt 0q ) qk dqk The final set of equations for the relaxed formulation for the two vertebral bodies depicted in figure (21) is X 1 0 m, 0 1 J, 0 (sin()y cos()y2 ) (cos(O)y sin(O)y,) mI 1 1 0 m Y 0 1 J, ; (sin(O)y, +cos(O)Y,2) (cos(O)Y, +sin(O),2)) Q1 Q2 jp1cos(o) sin(OfQ f1"l 3 cos(O) sin(O.) Yi), Q LJsin(O) cos(O) y, (SI) YJ sin(QO) cos(O]) Y%,2(s)Jj K (223) 2.2 Identification Problem in the TwoDimensional Example The system identification issue in our research is mainly to derive a mathematical model for the human spine based on clinically measured data. By essentially adjusting model parameters within this given model, we try to make the output of the simulation system coincide with the experimental data as well as possible. In the above section, we have derived the governing equations of the twobody system. But we wish to know how well this equation can describe the twobody example. In other words, we want to find the a set of the suitable probability measures that can fit our model with the experimental data. As we know, the mathematical model describes the relationship between the input signals and output signals. Considering the example we have, the inputs would be the different forces and moments applied on the vetebral bodies, while the output would be the vetebral body's orientation, position and movement. The outputs are affected not only by the system characterization, but also by the input signals. We denote the inputs by u(t) and outputs by y(t), respectively. The relationship is depicted in figure 23. All these signals are functions of time, and the value of the input at time t will be denoted as u(t). Often, in the identification context, only discretetime points are considered, since in the clinical experiment, the instruments used in the experiments typically record the signal just at discrete time instants. The identification of spine model is based on these experimental input and output data. Our human spine system identification is comprised of the following steps: Design the clinical experiment and collect inputoutput data from the experimental procedures. This is the first stepping stone in going from an experiment to the theratical model. However, these experimental data often need to be polished by removing trends and outliers, and selecting useful information from the original data. Many research has been done in this field, as noted in the literature we reviewed in the introduction. Select and define a spine model mathematically to represent the experiment. There are a set of candidate system descriptions, within which a model is to be found. Compute the best spine model characteristic parameters according to the inputoutput clinical experimental data and a given criterion of fit. Then, examine the obtained model's properties. In figure 23, the output y(t) is the system response to the input signal u(t), while d(t) is the desired output. In spine modeling, d(t) can be considered as the experiment data we obtained from the clinical experiment. By comparing the output of the given model with the desired output, we can adjust the characteristic parameter within the mathematical model, then minimize the error e(t) between the desired output d(t) and the output of identification model. Sd(t) u(t) ,,;'7 y(t) S olEL > COMPARE e(t) Figure 23 Identification procedure illustration In our simple example, the identification problem is a minimization problem over the probability measure space. In this twobody vertebral segment system, we need to find a stiffness distribution that could reconstruct the clinical experimental data. Suppose we have measurements (y,y) (which stand for the displacement and velocity measurements in our case). Then, we define the performance function as J(i)= y H,(q(u),q()) + j1 H2(q (), q())2 (224) The functions H1 and H2 estimate the experimental data from the model for a given Pu. So the the performance function defines the error E as described in figure 23. In this work we try to identify the stiffness coefficients of the system that minimize the performance equation. Recall the governing function we have derived in the first session. It's difficult to get a analytical solution in the form of H1 and H2 that we can substitute it into equation (224) directly. Numerical analysis is widely used today to solve this problem. It can provide the numerical answer while the analytical answer is difficult to obtain. There are two major computing procedures employed in solving this problem. One is to find a numerical solution to the ordinary differential equations arising from solving the Lagrange's equations. Another is to use a leastsquares method to calculate the spring stiffness coefficients. We obtained the governing equation of motion in the form of equation (223). This is a set of secondorder ordinary differential equations coupled to each other. However, these secondorder equations can be reduced to a system of simutaneous first order equations. Let's express one of the secondorder differential equation as a example, d2X dX = f(t, ) (225) dt dt The initial values can be described as X(to) = X0 and X'(to)= Xo. We can convert this to a pair of firstorder equations by the simple expedient of defining the derivative as a second function. 19 dX t= Y, X(t) = X,, (226) dt dY = f(tX, YO), Y,(to)= Xo (227) dt This pair of firstoder equation is equivalent to the original equations. Similar equivalent firstorder equations can be obtained for the other second order differential equations. So, equation (226) can be covered into 1 m, 1 Y, Y2 Q2 Y6 12 Q6 Y m, YO 0 Y, 1 Y, 0 1Y 1 0 Y12 1 0 1 1 1 (228) where the vector Y, contains six generalized coordinates q, namely (x, ,,,X,Y,,,0,,X,, Y,) And, ZK,( Q,t), stands for the Jacobian matrix multipied by the constraint vector 1=1 Sqj D as follows. 1 0 0 1 (sin(6),1 cos(OQ)y,) (cos(,)y, sin(,)y,2) x 1 0 0 1 (sin(0,)y,, + cos(8,)y,2) ( cos(,)y,, + sin(O)y,,2) {i ] cos(qs ) sin(6O ) (s,) XJ cos(8,) sin(O ) g (si) IYJ sin(8O) cos(,) 'Y(s)J sin(O,) cos(O ) r,,2(sl)j (229) At this point, we can apply numerical methods to solve these pairs of firstorder equations. Another numerical issue is related to the identification procedure, where we try to fit the clinical experiment data into the given mathematical model. Obviously, the problem involved is one of fitting nonlinear data. The overall algorithm describing the numerical solution procedure is described in figure 24. First, an initial guess of the probability measure vector/ pis made. That is, the initial value of the spring stiffness vector ko is given. By using the ODE solver in Matlab, we can obtain the numerical approximation at the output data. We then compare these data with the clinical experiment data. By checking the performance function, or error measure to a certain threshold, we decide to stop or go on to the next iteration step. The last results in this iterative process are estimates of the identified characterization parameters. That is, we obtain estimates of the spring stiffness distribution between the two vertebral bodies. Detailed numerical simulation procedures and simulation results will be explained in the following chapters. Figure 24 Numerical simulation procedure illustration 4:::::::::: CHAPTER 3 NUMERICAL EXPERIMENT In this chapter, we first give an introduction to common terms for the simulation tool Promechanica. Then, two numerical experiments are discussed in detail. The first one is a forward numerical simulation, and the second is an inverse dynamics identification experiment. 3.1 Multibody Dynamics Simulation ToolsPromechanica In this chapter, we will describe the development of a computer model of the human spine using the commercial dynamics program Promechanica. It has a particular high level language to allow modeling of rigid bodies, joints and constraints. It is possible using this language to simulate the movement of mechanisms. Additionally it is possible to measure forces and moments at required points. Since these quantities will be used in the remainder of the text, we introduce the basic terms within the programming language. Terminology: Body: Representation of a rigid body within the computer model. Connection: Joint used to connect two or more bodies together (pins, cylinder, ball, bearing, universal joint, sliders, free, weld, gimbals, planar and six degree of freedom). Load: Applies loads to points, joint axis or between two points on different parts (gravity, fixed loads, moments, static and dynamic friction). Driver: Used to specify the acceleration, velocity or position of a part or part joint axis (ramps or constant functions, cosine, table or combinations of these. Measure: Used to monitor changes in model parameters during the analysis. Mass primitive: A geometric shape (brick, cone, cylinder, plate, and sphere) which can be added to a part, changing its mass, center of mass and moments/products of inertia. The geometry of bodies are created by connecting points specified in the body fixed local coordinate system (LCS). Mass and inertial properties can then be assigned to the bodies either explicitly by the user or by using mass primitives and material defined within the program. These bodies can be connected to each other and to a ground point using a variety of joints to form a mechanical system. Drivers and measures can be assigned to each joint axis if so required. Surface geometry allows a more realistic visualization of the model. Surface is attached automatically to mass primitives and can be attached manually to bodies. The mass and inertial properties of bodies are entered explicitly by the analyst. Also there are some other capabilities to tailor the kinematic and dynamic aspects of the mechanism design. One important capability is the custom load, which is created outside of Promechanica Motion by the analyst. This is convenient for a user to model a more complex load. To analyze the stiffness characteristics of the human spine and simulate its movements, we would like to describe the joint characteristics between two vertebral bodies as a distribution of springs. There are no such joints can be used directly inside the Promechanica from it lists of buildin entities. The custom load is good way that we can define our specific model by using analystdeveloped loads. 3.2 Forward Simulation with Custom Load Recall the main procedure in the human spine modeling as described in Fig 31. First, we obtained the data from the clinical experiments. By using an approximation method, we can identify the mechanical characteristic as parameters in the model we construct. Then, we input these identified parameters into Promechanica, to simulate the motion of the experiment. Figure 31 Human spine modeling procedure If we consider the approximation procedure as the inverse dynamics problem, the Promechanica simulation provides the direct dynamics problem. "The direct dynamics problem, being the inverse of inverse dynamics solutions, should provide a way to check the inverse dynamics problem."[DLC1997] That is, if the approximate solution, calculated using our mathematical model in chapter 2, is used as input to Promechanica simulation with the same model, the outputs of the simulation should match the data as the input to the identification procedure. If instead, the output of simulation and the input to the inverse dynamics problem are different, then an error exists either in the fundamental formulation or the implementation of either the inverse dynamics or direct dynamics solution techniques. A forward numerical experiment was designed for this consistency check. Figure 32 A two vertebral body model built in Promechanica Consider the twobody system example we used in driving the governing equations. We construct the same twobody system example in Promechanica Motion by using published data together with information on the spine geometry imported from a commercially available source. Figure 32 shows the two vertebral bodies with geometry attached. In this forward simulation, we simply compare two mathematical models. One is the Matlab model based on the governing equations derived in chapter. Another is the one constructed in Promechanica. The same external forces are applied to identified model and Promechanica model. Also, we input the same mechanical parameters to these two models, like the same number of springs, the same attachment position, and the same stiffness probability measurements. The outputs of both systems are time responses to the same external forces. From the clinic viewpoint, the Promechanica simulation output can be considered as the clinical experimental data. By comparing the outputs of these two models, we can find if the Matlab model matches the Promechanica model. z Q1 Q2 X Figure 33 Generalized forces in twodimensional plane The external forces were added along the three generalized coordinates in this twodimensional example. They are generalized forces Q1 and Q2 along the Z and Y translational direction, and one generalized moment Q3 about the Xaxis, as shown in the above figure 33. The generalized forces could have different expressions. For example, in a simulation using 8 discrete springs, we use Q1 = A1 cos(27rft + )+ A2 cos(2rfh + ) Q2 = Bi COS(27rft + 0)+ B2 CO (2fh +0) Q3 = C1cos(27rf, +O)+C2 COs(27rfh +0) (31) where fl, fh stand for the low frequency part and high frequency part in the external forces. In a 16spring rotational movement simulation, we use 27 Q1 = A1x2 + Bx + C1 + Dx1 + E1x2 Q2 = A2x2 + B2x+ C2 + D2x 1 + E2x2 Q3 = A3x2 + B3 + C3 + D3x 1 + E3 2 (32) With different coefficients Ai, Bi, Ci, Di, Ei, we can give different force inputs. After the simulations have been run, we compare the numerical results from both models as shown figures 34 through 315. 8sDrina numerical simulation comDarison 50 100 150 200 250 Figure 34 Ytranslational trajectory numerical results from Promechanica 50 100 150 200 250 Figure 35 Ytranslational trajectory numerical results from Matlab model 40 39 5 39 38 5 38 37 5 37 0 50 100 150 200 250 Figure 36 Ztranslational trajectory numerical results from Promechanica 40 39 5 39 38 5  38 37 5 37 0 50 100 150 200 250 Figure 37 Ztranslational trajectory numerical results from Matlab model 15 05 0 05 1 5 0 50 100 150 200 250 Figure 38 Xrotational trajectory numerical results from Promechanica 15 05 0 0 5  1 1 5 0 50 100 150 200 250 Figure 39 Xrotational trajectory numerical results from Matlab model 16sDrina numerical simulation comDarison x 10 50 100 150 200 250 Figure 310 Ytranslational trajectory from Promechanica model x 10 50 100 150 200 250 Figure 311Ytranslational trajectory from Matlab model 38 308 38 307  38 306 38 305 38 304 38 303 38 302 38 301 38 3 38 299 0 50 100 150 200 250 Figure 312 Ztranslational trajectory from Promechanica model 38 309 38 308  38 307 38 306 38 305 38 304 38 303 38 302 38 301 38 3 38 299 0 50 100 150 200 250 Figure 313 Ztranslational trajectory from Matlab model x 10 1 6 14 12 08 06 04 02 0 50 100 150 200 250 Figure 314 Xrotational trajectory from Promechanica model 0 016 0 014 0 012 0 01 0 008 0 006  0 004  0 002 0 0 50 100 150 200 250 Figure 315 Xrotational trajectory from Matlab model In many different simulation cases, the numerical trajectory from the Matlab model matches the results from the Promechanica Motion simulation, as shown in the above figures. At last, we used the forward simulation capability to solve an identification problem. That is, we use Promechanica Motion to get an output trajectory. This data can be considered as the clinical, experimental data. Then, the experimental trajectory data is supplied to the identification procedure in Matlab. We can identify the spring stiffness distribution that best represents the trajectory obtained from Promechanica. This is the same as the identification procedure we discussed in the chapter 3. We compare the original input spring stiffness coefficients use to generate the experimental trajectory and the identified stiffness coefficients in the following tables. Table 31 8 Springs identification Spring 1 2 3 4 5 6 7 8 Input 0.2 0.2 0.1 0.3 0.15 0.12 0.08 0.22 Identify 0.1962 0.1989 0.1045 0.29 0.152 0.1227 0.0818 0.2239 Table 32 16 Springs identification Spring 1 2 3 4 5 6 7 8 Input 0.3 0.15 0.12 0.08 0.22 0.25 0.16 0.3 Identify 0.2913 0.1413 0.11 0.0743 0.2151 0.2506 0.1608 0.3048 Spring 9 10 11 12 13 14 15 16 Input 0.1 0.07 0.09 0.18 0.2 0.2 0.2 0.1 Identify 0.1052 0.0776 0.0986 0.19 0.2099 0.1915 0.1919 0.1071 When we input the identified spring stiffness into the simualtion software Promechanica, we can visualize of the verte bral body movent. Figure 316 Translation movement simulation Figure 317 Rotation movement simulation 3.3 Identification of Pin Joint Position Using Inverse Simulation The last examples considered trajectory matching problems. That is, we sought the distributions to minimize the difference between observed data and calculated outputs. One of the critical properties of a relaxed formulation is that it should be able to represent the classical results as a special case. A numerical experiment is used to show this is the case. Suppose we use a pin joint to connect two vertebral bodies. We want to know if we can use the relaxed formulation to approximate this kind of the classical mechanical joint. Obviously, this is an identification problem. b1,2 e, e2 Figure 318 vertebral body connected by pin Similar to the previous numerical experiment, we use the simulation tool Promechanica to construct a twodimensional, twobody system. By importing a two vertebral segment system, we use a pin joint to connect the vertebrae together as shown in figure 318. The center of the pin is at a certain point between the two bodies in the sagittal plane. Since our spine model is a threedimensional geometry model, we use drivers to apply constraints on three coordinates (one translational and two rotational), and thus simplify the three dimensional system to a two dimensional system. After the dynamical simulation, we get the trajectories of the movement for the system connected via the classical pin joint. Then we want to find a distribution of springs that can replace the pin joint and can obtain the same system outputs. The attachment points for each spring are fixed on the surface of each vertebral body. The interval distance between these attachment points is equal. We can make the distance finer by adding more locations to connect the vertebral bodies. By using the identification method as we discuss in the chapter three, we can identify the stiffness coefficient value of each spring. 38 The 2 spring forces case Discrete Probability Density Function 2 Springs Case data ' data2 data3 [ data4 x1014 A 10 15 20 25 30 35 Position of Spring Attaching Point marked by 40 B Figure 319 Probability density distribution of spring stiffness Cumulative Distribution of Spring Stiffness x1014 2 Springs Case A 10 15 20 25 30 35 40 B Position of Springs Attaching Points Figure 320 Cumulative distribution of spring stiffness I I I I x 10 x10 39 The 4 spring forces case Discrete Probability Density Function 4 Springs Case 11 20 25 '30 35 Position of Spring Attaching Point Marked By '* Figure 321 Probability density distribution of spring stiffness Cumulative Distribution of Spring Stiffness S103 4 Springs Case 20 25 30 3 Position of Springs Attaching Points Figure 322 Cumulative distribution of spring stiffness x 10 x10 40 The 8 spring forces case Discrete Probability Density Function 8 Springs Case 05 0 A 10 1 2 B Position of Spring Attaching Point Marked By '* Figure 323 Probability density distribution of spring stiffness Cumulative Distribution of Spring Stiffness x 103 8 Springs 4 5         4 3 35 a 3 Q  Co U) S 25 2 S15 O 05 0 1 ________________ A 10 15 20 25 30 35 40 B Position of Springs Attaching Points Figure 324 Cumulative distribution of spring stiffness 41 The 16 scoring forces case Discrete Pobability Density Function 16 Springs Case x 1 0 x10 t bu A 10 15 20 25 30 35 40 B Position of Spring Attaching Point Marked By '* Figure 325 Probability density distribution of spring stiffness Cumulative Distribution of Spring Stiffness x 103 16 Springs Case 6  Co 0 C E S1. A 10 15 20 25 30 35 40 B Position of Springs Attaching Points Figure 326 Cumulative distribution of spring stiffness x 3 Figure 326 Cumulative distribution of spring stiffness 42 The 32 spring forces case Discrete Probability Density Function 32 Springs Case x 10 x10 35 o S 3 25 S 2 015 o 053 O 10 15 20 25 30 35 40 Probability Density Distribution of Springs Figure 327 Probability density distribution of spring stiffness n Pr b b lt e st iti u i n o p i g Cumulative Distribution of Springs Stiffness 32 Springs Case 0 012 0 01 fc U) S0 008 0 S0 006 O 0 004 S0 002 0 0 A 10 15 20 25 30 35 40 Position of Springs Attaching Points Figure 328 Cumulative distribution of spring stiffness . . . . . . ..I I I I I 43 The 64 sDrinu forces case Discrete Probability Density Function 64 Springs Cases x 10 x10 A 10 15 20 25 30 35 40 B Position of Spring Attaching Point Figure 329 Probability density distribution of spring stiffness Cumulative Distribution of Spring Stiffness 64 Springs Case 0 02 0 0 18    0 016 o 0 016 . 0 014 001 S0 012 0 008 S0 006 E0 004 0 002 A 10 15 20 25 30 35 40 B Position of Springs Attaching Points Figure 330 Cumulative distribution of spring stiffness x 10 44 The 128 spring forces case Discrete Probability Density Function 128 Spring Cases A 10 15 20 25 30 35 40 B Position of Spring Attaching Point Figure 331 Probability density distribution of spring stiffness Cumulative Distribution of Springs Stiffness 128 Springs Case 0 04 0 035 c 0 03 . U) 0 025 0 02 VQ 0 0 015 . o o FO 001 . 0 005 A 10 15 20 25 30 35 40 B Position of Springs Attaching Points Figure 332 Cumulative distribution of spring stiffness The 256 spring forces case Discrete Probability Density Function 256 Springs Case 0 005 0  04 C 0 03 0 02 S0 01 0 A 10 15 20 25 30 Position of Spring Attaching Point 35 40 Figure 333 Probability density distribution of spring stiffness Cumulative Distribution Function of Spring Stiffness 256 Springs Case 0 35 uE S0 25 . 0 02  015 C > 01 E 0 0 05 15 20 25 30 35 40 B Position of Springs Attaching Points Figure 334 Cumulative distribution of spring stiffness 46 If we consider the spring stiffness coefficients as the probability function, the value of stiffness can be normalized between 0 and 1. By plotting all the spring coefficients probability distribution and cumulative distribution are obtained. We can see the identified numerical results clearly. With a number of springs connecting the two bodies, the lager stiffness coefficients are more likely move to the position where the pin joint exactly is. Spring Stiffness Distribution 2,4,8,16,32,64,128,256 springs case . 0 07 006 a 005 S004 0 03 S002 001 8 Figure 335 Probabilty 0 35 03  S0 25 02 S0 15 01 E C 0 05 20 32 44 Position of Springs Attaching Points density distribution of different spring numbers Cumulative Distribution of Spring Stiffness Comparation 2,4,8,16,32,64,128,256 springs cases 20 25 30 position of springs attaching points Figure 336 The cumulative distribution of different spring numbers Instead of the pin joint between the two vertebral boides, we input these different spring stiffness distribution into a Promechanica Motion model with the same external forces, and simulated the system to get the resulting trajectories of the system. By comparing the error between the original system and the identified system, we found that better correlation is achieved using a finer distribution vertebral segment. The following figures show the relationship between the approximation error and different number of the springs. 28 5 28 27 5 27 26 5 26 25 5 25 24 5 24 23 5 0 10 20 30 40 50 60 70 Figure 337 Error vs. different number of springs The identification above shows that using a finer spring distribution yields a smaller approximation error between identified model and the ideal model. Also, the position of the pin joint can be identified more clearly. This simulation illustrates an important fact: the joint relaxation technique includes, as a special case, the classical pin joint constraint. CHAPTER 4 SCOLIOSIS MODEL DATABASE CONSTRUCTION In chapter 4, we summarize the development of a structural and geometric database that will support the future study of scoliosis. We begin with a discussion of the spine stiffness database construction. Finally, we show an example of a twodimensional vertebral segment whose stiffness characteristics are modeled using the database. The tables in that constitute the database can be found in the Appendix. 4.1 Geometry In order to construct a mathematical model of the human spine, we need to construct a threedimensional scoliosis model that can be used in Promechanica geometry. From the literature, we have found many types of human spine models that have been published. Based on these published models, we have built an electronic database. A threedimensional spine model and coordinate system is constructed for use with the dynamics simulation tool Promechanica. This simulation package is selected because scoliosis is a threedimensional deformity. It makes the spine deform in a way that requires a six DOF model. However, by adding the corresponding constraints and drivers, we can still simplify it to a twodimensional model in some specific cases. The spine can be idealized as a collection of rigid bodies interconnected by some joints or elements. The vertebrae are modeled as rigid bodies, since their deformation is negligible compared to the intervertebral tissue. The geometry of an individual vertebral body can be constructed from the published data by using the Promechanica geometry attachment function. In the global coordinates system, the sacrum will be considered as the fixed ground of the whole model. Each vertebral body above it, from L5 to T1, can be imported individually by using the subassembly function. In each subassembly system a local coordinate system is defined. Every vertebral body is imported into its corresponding subassembly system. We can describe the geometry characteristics in this local subassembly system. Then, we import subassemblies into the global system one by one, according to the published data defining the shape of scoliosis. The global coordinates of each primary node in the subassembly system can be computed by using the transform matrix we derived in chapter 3. The orientation values in the sagittal plane are derived from the published data. By inputting the mechanical properties, like the mass, and other connecting joint constraints, a spine model is constructed in Promechanica. In table 41, some geometry data are given that characterize a spine that is scoliotic. Table 41 Scoliosis curve data Vertebra Sagittal Vertebra Inferior Disc Orientation (deg) Height (mm) Height (mm) L5 8 22 15.7 L4 2 28 14 L3 11 26 12.2 L2 21 26 11.4 L1 21 24 10 T12 18 24 8.4 T11 17 20 6.8 T10 17 20 5.1 T9 15 20 4.7 T8 8 20 4.5 T7 0 20 4.0 T6 3 20 3.2 T5 4 20 2.6 T4 9 20 2.2 T3 15 18 2.7 T2 20 18 3.1 T1 24 16 4.5 Table 42 Coordinates of origin nodes Vertebra globalX globalY globalZ L5 0 0 0 L4 0 31.19 37.46 L3 0 32.05 77.73 L2 0 41.45 112.73 L1 0 51.50 145.50 T12 0 64.00 175.40 T11 0 68.02 202.95 T10 0 76.32 226.65 T9 0 83.62 250.35 T8 0 88.72 274.15 T7 0 90.72 298.05 T6 0 90.72 321.25 T5 0 89.52 343.75 T4 0 87.92 365.95 T3 0 81.28 387.22 T2 0 75.78 407.63 T1 0 66.04 427.90 Table 42 shows the threedimensional coordinates of the original nodes of each local subassembly systems. The subassembly systems are imported into the global system by pointing their original nodes and orientations. Using this information the global coordinates of all the individual points can be calculated. The relative position of any arbitrary points on the superior and inferior vertebral bodies is then found in this global system via the transformation matrix between the local coordinate systems and global coordinate system. A spine model form T1 to L5 is generated in Promechanica, as shown in Fig 41 Fig 41 Spine Model in Promechanica Motion To make a simpler geometry model, we pick up 17 primary nodes, which could describe the basic geometric outline of vertebra. And the body of vertebra is simplified to a cylinder. The other primary nodes can be picked up on the pedicle, superior articular process, transverse process, and inferior articular process. Using spline curves to link these nodes, we can get a basic shape of a vertebral body. Fig 42 A Vertebral Model with Complex Surface Fig 43 A Simple Geometry Vertebral Model 8 ~ 15  16 < ^__^' 0 3 14  1 Fig 44 Right Lateral View of the Primary Nodes I W1 %,i,fl, i~h.O nnv CiviI Maim rrF FqrTEl 16 7I,~ Fig 45 Top View of the Primary Nodes I W11%,i...., n........O ,n=,, l.....n. m r)r FqrEi 3 16 C1 .4 I   I Fig 46 Bottom View of the Primary Nodes I W1 %,i,fl, i~h.O nnv flivim.nil Mam rFI qrTE Fig 47 Simplified Geometry Spine Model A simplified geometry spine model is constructed as shown in Fig 47. One objective of this biomechanical simulation is to study the implants that attach to the bone. We have derived the position of the pedicle screws, transverse process hooks, and lumina hooks based on these primary nodes. Because the internal fixation devices used in spinal surgery are not strictly sensitive to their position, we can model the geometry and the position of these implant devices by some simple computations using our biomedical experiment database. The geometric database describing the orientation of the screws and hooks are described in the Appendix. 4.2 Stiffness Many researchers have studied on the stiffness properties of the connective tissues between the vertebral bodies. The literature review showed that different stiffness values have been reported in different clinical experiments. The diversity in results is partly because of the differing experimental methods and instruments used in these experiments. Also, the experimental data will be different when testing different specimens. Among these experiments, Panjabi et al made extensive studies of the nonlinear viscoelastic mechanical behavior of the human spine. The threedimensional behavior of the spine is presented in the form of average loaddisplacement curves for different vertebral bodies. After examining these discrete experiment data, we have identified and tabulated a nonlinear loaddisplacement function to describe the stiffness characteristics of each Y t TMY RY SRZ Txl C Fig 48 Illustration Stiffness Property in Threedimensional Coordinates System vertebral level. By inputting these simplified stiffness curves into the model we constructed in Promechanica, we can simulate the clinical experiment and reconstruct the relationship between the load and vertebrae movements using numerical simulation tools. This simple test illustrates the utility of the database. Fig 48 illustrates the threedimensional coordinates system that is necessary to describe the threedimensional physical properties of the spine. The origin of the coordinate system is located at the infer posterior corner in the midsagittal plane of the vertebral body. The Xaxis is directed to the left and is perpendicular to the sagittal plane. The Yaxis is directed superiorly, and the Zaxis is directed anteriorly. The broad arrows show the pure moments and thin, straight arrows show the translations on the three rotational axes and three translational axes respectively. In the following figures, we show the loaddisplacement curves between each level of the lumbar spine in threedimensional system. They are from the location between first and second lumbar vertebrae through location between the fifth lumbar and first sacral vertebrae. For each location, the loaddisplacement relationships are described according to the different type of loads applied on the vertebrae, pure moments of flexionextension, lateral axial torque, and bilateral bending, respectively. 60 Loaddisplacement curves are shown at each of the five locations on the vertebral lumbar due to the application of flexion and extension moments. The main motion, defined as the rotation in the direction of the applied moment, is indicted by the heavier line. The figures also show the coupled motions, defined as all rotations and translations other than the main rotation. L1L2 12 10 8 6 4 V         I I I I I I 0 2   2 ........... 0 I i i .. .... ..... ............. .............. .............. .............. .............. I .............. .............. .............. .............. .............. i .............. .............            1 1 6 8 10 19 12 10 8 6 Extension 4 2 0 2 4 6 8 10 12 Flexion MOMENT (Nm) ] !A Fig 49 Flextionextension between L1 and L2 J* 61 L2L3 12 1 0   S ....... .................................................................................. 0 4 ............. ..... ..... .............. . ... ........... ........... .............. .............. .............. .............. ............. S2 .... ....... ............. ............ ............ ... O. 0ee   12 I I I I 6 . . I ........ .... .... ... ..... ... ... ...         I     o A o 4 6 8 10 12 10 8 6 4 2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) Fig 410 Flextionextension between L2 and L3 62 L3L4 ............. ........... ......... .............. ....... T .............. r ... ........... ... ......... ...... .............. .. ................ ..... i i i I  .....i ..............................I.................................. . . . ...  .. . . .... . . ..    I                                                  10 8 6 4 2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) Fig 411 Flextionextension between L3 and L4 6 8 10 12 1  63 L4 L5 14 T. T .r A * 10 8 6 4 2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) Fig 412 Flextionextension between L4 and L5 12L 12 64 L5S1 ............. .............. .............. .............. ....... .............. .............. ............ ..............i .............. .............. ............ ........... ............. ... ........  ....... ..... .............. .............. .............. .............. .............    &        T r 1          i         i    10 8 6 4 2 0 2 4 6 8 10 12 Extension Flexion MOMENT (Nm) Fig 413 Flextionextension between L5 and S1 12 L 12 Loaddisplacement curves are next shown at each of the five locations on the lumbar spine due to the application of right and left axial torques. The main motion, defined as the rotation in the direction of the applied moment, is indicted by the heavier line. On the following figures, also shown are the coupled motions, defined as all rotations and translations other than the main rotation. E E 1 0) 0) 0L 7a I O 1 O 2 .A L1L2 _ .......... .............. ............ .............. .............. .. ............... _ ........  ^^+ J  I II 12 10 8 6 Right Torque 4 2 0 2 4 6 8 10 12 Left Torque MOMENT (Nm) Fig 414 Right and Left Axial Torque between L1 and L2 L2L3 ........ ............ .............. .............. .............. ............................ .............. .............. .............. .............. i.......... _ 3D      "S      ^ E o 0 z 12 10 8 6 4 2 0 2 4 6 8 10 12 Right Torque Left Torque MOMENT (Nm) Fig 415 Right and Left Axial Torque between L2 and L3 I 2 . Fig 415 Right and Left Axial Torque between L2 and L3 67 L3L4 i ............. ............. .............. .............. .... .............. ............................................... .............. ............. ... 1 T F F1 1 1 V4 10 8 6 Right Torque 4 2 0 2 4 6 8 10 12 Left Torque MOMENT (Nm) Fig 416 Right and Left Axial Torque between L3 and L4 68 L4L5 E 2    4    i 2 .........................................   I.............. .............. ..............I ..............  .............. .............. . i i I II I I o 4 12 10 8 6 4 2 0 2 4 6 8 10 12 Right Torque Left Torque MOMENT (Nm) Fig 417 Right and Left Axial Torque between L4 and L5 69 L5S1         ............ ............ ......... ..... ....... ....... .............. .............. .. ............ .............. .. .... ....... I I I I I I_ I I I I 2 10 8 6 4 2 0 2 4 6 8 10 12 Right Torque Left Torque MOMENT (Nm) Fig 418 Right and Left Axial Torque between L5 and S1 70 Loaddisplacement curves below correspond to each of the five locations on the lumbar spine due to the application of left and right lateral bending moments. The main motion, defined as the rotation in the direction of the applied moment, is indicted by the heavier line. The figures also show the coupled motions, defined as all rotations and translations other than the main rotation. 8 6 4 E 0 1 0) 01 2 I 4 6 L1L2 ............. .............. ...............  .............. ..............  ..............  .............. ............    i i i i 0  _ ............... I .  .....................       I .............. I  I ......... ....                   ... .. .. .. .. .. .. .. .. . 10 8 6 Left Bending 4 2 0 2 4 6 Right Bending 8 10 12 MOMENT (Nm) Fig 419 Left and Right Bending Moment between L1 and L2 71 L2L3 ............. I.............. .............. .............. .....i. ...... .. .............. TT.. ............. ........... 10 8 6 Left Bending 4 2 0 2 4 6 8 10 12 Right Bending MOMENT (Nm) Fig 420 Left and Right Bending Moment between L2 and L3 72 L3L4 6      4 E  .. .. .. .. .. .... ....... .. .... . E 0 4 .. .. .. i ..   ............   2 I I 4 6 8 12 10 8 6 4 2 0 2 4 6 8 10 12 Left Bending Right Bending MOMENT (Nm) Fig 421 Left and Right Bending Moment between L3 and L4 73 L4L5 6   E I i i i i i i i i . i S 2 .......... ....... .7.. ..... 0 I 6 4      8 12 10 8 6 4 2 0 2 4 6 8 10 12 Left Bending Right Bending MOMENT (Nm) Fig 422 Left and Right Bending Moment between L4 and L5 74 L5S1 S. ........... ............. ............. ............ .............. ............. ............. ...... ...... 4 ......... ............. I ..... ............ 6          o I i   12 10 8 6 4 2 0 2 4 6 8 10 12............. 28 Left Bending Right Bending MOMENT (Nm) Fig 423 Left and Right Bending Moment between L5 and S1 4.3 A Simple Demonstration of the Stiffness Database After the spring stiffness database for the spine (Lumbar and thoracic) has been built, we simply approximate the stiffness coefficients as the nonlinear function for each vertebral level. When we input these nonlinear functions into a model in Promechanica, we can simulate the clinical experiment with the same external moments. Take an example consisting of an L5S1 twobody system. Define X as the relative angular rotation between the bodies. The unit of angle is degree; the unit of displacement is Nm. We input the conditional spring stiffness coefficients as followings: Phasel: 5.8 < X < 5 (Degree) Initial condition = 5 (Degree) Spring stiffness = 5/0.8 (Nm/mm) Phase2: 5 < X < 1.7 (Degree) Initial condition = 1.7 (Degree) Spring stiffness = 5/3.3 (Nm/mm) Phase3: 1.7 < X < 1.7 (Degree) Initial condition (Degree) Spring stiffness = 0 (Nm/mm) Phase4: 1.7 < X < 6.8 (Degree) Initial condition = 1.7 (Degree) Spring stiffness = 2.5/5.1 (Nm/mm) PhaseS: 6.8 < X < 8.5 (Degree) Initial condition =6.8 (Degree) Spring stiffness = 7.5/1.7 (Nm/mm) Similar to the twobody system we built in the chapter 4, we construct the two body S1L5 system in Promechanica, fixing S1 on the ground, using a 6 DOF joint to connect S1 and L5. We add a spring on the rotational joint axis between the two bodies. By adding the appropriate constraints to the system, we can keep the movement of L5 with admissible limits defined in the database. It should be noted that the nonlinear stiffness profiles may bring very large forces into the system. Also, at the point where we change the stiffness coefficient suddenly, we may also bring large forces or high frequency forces to bear on system. So, we add a damper into the system, (c = 10) to make the system more stable. The reconstruction result is depicted in the Figure 424. Other stiffness database input code can be found at appendix. 1 0 T   i  i     1  6. Ssimulition Angle  10 8 6 4 2 0 2 4 6 8 10 Extension/Flexion MOMENT Fig 424 LoadDisplacement Compare between Database and Reconstruction CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS This thesis presented a method for simulation of the human spine and described the creation of a database for human spine simulation. By using fundamental multibody dynamics theory and computer simulation software, a forward numerical simulation and an inverse dynamics identification experiment are described. The identification method for human spinal motion is based on a relaxation of the original governing equations derived using Lagrange's method. The results of the numerical experiment show that the position of the pin joint can be identified as a special case of the identified probability measure. The reconstructed simulation results match the original simulation successfully. Using computer simulation program Promechanica, we also constructed a structural and geometric human spine database based on some published human spine data. The stiffness characteristics of intervertebral bodies were incorporated into the database. Experimental results from the literature were reconstructed from the database by using Matlab. The database should provide the foundation for future scoliosis model simulation and analysis. APPENDIX A GEOMETRY OF MODELING NODES Table Al Lumbar 5 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 0 0 2 0 16 11 0 14.3 13.1 3 0 32 0 0 31.7 4.5 4 0 16 11 0 17.4 8.7 5 46 43 3 46 42.2 9 6 46 43 3 46 42.2 9 7 3 51 0 3 50.5 7.1 8 3 51 0 3 50.5 7.1 9 3 52 11 3 53 3.6 10 3 52 11 3 53 3.6 11 21 47 9 21 45.3 15.5 12 21 47 9 21 45.3 15.5 13 21 47 12 21 48.2 5.3 14 21 47 12 21 48.2 5.3 15 0 61 5 0 61.1 3.6 16 0 72 15 0 73.4 4.8 17 0 61 20 0 63.2 11.3 Table A2 Lumbar 4 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 3 38.5 2 0 20 14 0 17.5 52.1 3 0 40 0 0 37 37.4 4 0 20 14 0 16.5 24.1 5 45 43 8 45 52.2 44.9 6 45 43 8 45 52.2 44.9 7 3 55 0 3 55 36.8 8 3 55 0 3 55 36.8 9 3 58 14 3 58.5 22.6 10 3 58 14 3 58.5 22.6 11 21 61 14 21 58.5 50.7 12 21 61 14 21 58.5 50.7 13 21 61 16 18 45.3 15.5 14 21 61 16 18 45.3 15.5 15 0 73 1 0 69.9 35.3 16 0 92 17 0 88.3 18.6 17 0 73 21 0 69.2 15.3 Table A3 Lumbar 3 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 0.3 77.7 2 0 19 13 0 21.4 84.8 3 0 37 0 0 36.6 70.6 4 0 19 13 0 16.5 61.3 5 52 62 5 52 62.1 70.8 6 52 62 5 52 62.1 70.8 7 3 52 0 3 51.3 67.8 8 3 52 0 3 51.3 67.8 9 3 55 10 3 52.4 57.4 10 3 55 10 3 52.4 57.4 11 16 58 13 16 59.7 79.4 12 16 58 13 16 59.7 79.4 13 16 58 23 18 58.4 50.7 14 16 58 23 18 58.4 50.7 15 0 76 0 0 74.3 63.2 16 0 97 9 0 93.3 50.4 17 0 76 18 0 71.5 45.5 Table A4 Lumbar 2 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 9.7 113.7 2 0 18 13 0 31.2 119.4 3 0 36 0 0 43.3 100.8 4 0 18 13 0 21.9 95.1 5 48 60 7 48 68.2 98.8 6 48 60 7 48 68.2 98.8 7 3 51 0 3 57.3 95.5 8 3 51 0 3 57.3 95.5 9 3 53 9 3 56 86.3 10 3 53 9 3 56 86.3 11 16 58 16 16 69.6 107.9 12 16 58 16 16 69.6 107.9 13 15 58 21 13 59.7 79.4 14 15 58 21 13 59.7 79.4 15 0 71 4 0 77.4 92 16 0 91 9 0 91.4 72.7 17 0 71 14 0 71 75.2 Table A5 Lumbar 1 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 22.2 146.4 2 0 18 12 0 43.3 151.2 3 0 36 0 0 55.8 133.5 4 0 18 12 0 34.7 110.5 5 42 54 6 42 74.8 132.7 6 42 54 6 42 74.8 132.7 7 3 50 0 3 68.9 128.5 8 3 50 0 3 68.9 128.5 9 3 51 7 3 67.3 121.6 10 3 51 7 3 67.3 121.6 11 16 56 16 16 80.2 141.3 12 16 56 16 16 80.2 141.3 13 14 56 20 13 69.6 107.9 14 14 56 20 13 69.6 107.9 15 0 69 4 0 88.1 125.4 16 0 89 3 0 104.2 111.8 17 0 69 12 0 82.3 110.5 Table A6 Thoracic 12 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 34.7 176.3 2 0 16 12 0 53.6 182.8 3 0 32 0 0 65.1 166.4 4 0 16 12 0 46.2 160 5 25 55 5 25 88.5 164.1 6 25 55 5 25 88.5 164.1 7 3 50 4 3 83.5 164.7 8 3 50 4 3 83.5 164.7 9 3 48 4 3 79.1 157.7 10 3 48 4 3 79.1 157.7 11 15 46 15 15 83.1 176.4 12 15 46 15 15 83.1 176.4 13 12 50 20 13 80.2 141.3 14 12 50 20 13 80.2 141.3 15 0 62 0 0 93.7 157.2 16 0 76 7 0 104.8 146.2 17 0 62 15 0 89 142.9 Table A7 Thoracic 11 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 43.6 203.7 2 0 17 10 0 62.8 208.3 3 0 34 0 0 76.1 193.8 4 0 17 10 0 56.9 189.2 5 30 60 10 30 103.9 195.7 6 30 60 10 30 103.9 195.7 7 3 47 8 3 90.9 197.6 8 3 47 8 3 90.9 197.6 9 3 48 1 3 89.2 188.7 10 3 48 1 3 89.2 188.7 11 15 44 14 15 89.8 204.2 12 15 44 14 15 89.8 204.2 13 15 48 13 15 86.1 176.4 14 15 48 13 15 86.1 176.4 15 0 62 1 0 102.6 184.6 16 0 77 10 0 114.3 171.6 17 0 62 14 0 98.8 172.2 Table A8 Thoracic 10 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 51.9 227.4 2 0 16 10 0 70.1 232.3 3 0 32 0 0 82.5 218 4 0 16 10 0 64.3 213.2 5 30 62 4 30 112.4 213.1 6 30 62 4 30 112.4 213.1 7 3 45 5 3 96.4 219 8 3 45 5 3 96.4 219 9 3 46 2 3 95.3 212 10 3 46 2 3 95.3 212 11 13 44 15 13 98.4 228.9 12 13 44 15 13 98.4 228.9 13 15 46 10 15 92.8 204.2 14 15 46 10 15 92.8 204.2 15 0 63 8 0 109.8 201.3 16 0 77 21 0 119.4 184.8 17 0 63 21 0 106 188.9 Table A9 Thoracic 9 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 59.2 251 2 0 16 10 0 77.2 256.5 3 0 32 0 0 90.1 242.7 4 0 16 10 0 72.1 237.2 5 33 66 0 33 122.9 233.9 6 33 66 0 33 122.9 233.9 7 3 50 5 3 108.8 242.9 8 3 50 5 3 108.8 242.9 9 3 52 2 3 108.9 235.6 10 3 52 2 3 108.9 235.6 11 12 45 14 12 106.3 252.9 12 12 45 14 12 106.3 252.9 13 13 48 10 13 101.4 228.9 14 13 48 10 13 101.4 228.9 15 0 68 12 0 121.8 221.8 16 0 83 24 0 133.1 206.3 17 0 68 24 0 118.7 210.2 Table A10 Thoracic 8 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 64.3 274.9 2 0 16 10 0 81.5 282.6 3 0 32 0 0 96 270.4 4 0 16 10 0 78.7 262.8 5 33 64 3 33 128.1 268.9 6 33 64 3 33 128.1 268.9 7 3 50 5 3 114.5 272.9 8 3 50 5 3 114.5 272.9 9 3 52 2 3 115.5 265.7 10 3 52 2 3 115.5 265.7 11 12 46 13 12 111.7 281.3 12 12 46 13 12 111.7 281.3 13 12 48 12 12 109.3 252.9 14 12 48 12 12 109.3 252.9 15 0 68 16 0 129.4 249.6 16 0 82 35 0 140.6 228.8 17 0 68 26 0 128 239.7 Table Al Thoracic 7 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 66.3 298.8 2 0 15 10 0 81.3 308.8 3 0 30 0 0 96.3 298.8 4 0 15 10 0 81.3 288.8 5 34 58 5 34 124.3 303.8 6 34 58 5 34 124.3 303.8 7 3 45 5 3 111.3 303.8 8 3 45 5 3 111.3 303.8 9 3 47 2 3 113.1 296.8 10 3 47 2 3 113.1 296.8 11 11 43 13 11 109.3 311.8 12 11 43 13 11 109.3 311.8 13 12 45 9 12 114.7 281.3 14 12 45 9 12 114.7 281.3 15 0 63 17 0 129.3 281.8 16 0 76 32 0 142.3 266.8 17 0 63 28 0 129.3 270.8 Table A12 Thoracic 6 node localX localY localZ globleX globleY globleZ 1 0 0 0 0 66.3 322 2 0 15 10 0 80.8 322.8 3 0 30 0 0 96.3 323.6 4 0 15 10 0 81.8 312.8 5 34 58 6 34 123.9 331 6 34 58 6 34 123.9 331 7 3 46 6 3 111.9 330.4 8 3 46 6 3 111.9 330.4 9 3 47 2 3 113.3 322.4 10 3 47 2 3 113.3 322.4 11 10 43 13 10 108.6 337.4 12 10 43 13 10 108.6 337.4 13 11 45 8 11 112.3 311.8 14 11 45 8 11 112.3 311.8 15 0 63 13 0 129.9 312.3 16 0 77 29 0 144.7 297 17 0 63 24 0 130.5 301.3 