Abstract
A 3D dynamic computer model for the movement of the head-neck complex is presented. It incorporates anatomically correct information about the diverse elements forming the system. The skeleton is considered as a set of interconnected rigid 3D bodies following the Newton-Euler laws of movement. The muscles are modeled using Enderle's linear model, which shows equivalent dynamic characteristics to Loeb's virtual muscle model. The soft tissues, namely, the ligaments, intervertebral disks, and facet joints, are modeled considering their physiological roles and dynamics. In contrast with other head and neck models developed for safety research, the model is aimed to study the neural control of the complex during fast eye and head movements, such as saccades and gaze shifts. In particular, the time-optimal hypothesis and the feedback control ones are discussed.
1. Introduction
The computer modeling of the head-neck complex has been widely addressed for more than thirty years. Special interest arises from the number of deaths related to injury of the head-neck system during vehicle accidents [1–4]. Thus, most of the models of the head-neck system are intended to represent the behavior of the system from an external point of view; that is, many studies focus on the response of the complex to external forces, such as those encountered during car crashes or in extreme loading conditions [1, 2, 4–6]. Other models are intended to study the stability and conditions of the cervical spine to assess the stage of patients with spinal cord injury or related pathologies [7].
Tien and Huston report that the first models of the head-neck complex developed in the sixties were based on simple pendulum models and evolved to multibody models in the eighties [4]. In these multibody models, the skull and cervical vertebrae are considered rigid bodies connected to muscles which exert forces over them [2–4].
The muscles of the neck, which generate head movements and help maintain the stability of the cervical spine, have received most of the attention. Diverse models have been proposed to represent the dynamics of the neck muscles [8]. However, many of the models presented up to the middle nineties did not take into account the important role of soft tissues and their complex biomechanical behavior.
The soft tissues of the cervical spine consist of a variety of structures connecting and surrounding the osseous elements of the cervical vertebral column [9]. The main categories of soft tissues include the spinal ligaments, intervertebral discs, facet joints, and uncovertebral clefts. Their general function is to enable and limit movement of the cervical spine. Due to the complexity and nonrigidity of these structures, many authors have chosen to model them by using finite elements approach [7, 9, 10]. Finite element (FE) methods are used in the numerical solution of field problems, where the spatial distribution of one or more dependent variables is to be established. FE has a number of advantages over other numerical analysis methods (e.g., see Cook [11]). FE models have been used for the study of the cervical spine mechanics during the last two decades [12]. Anatomical details are often based on data from computed tomography (CT). The FE method is useful, not only to study kinematics but also stresses and strains in cervical tissues. Many FE models are parameterized, which allows for studies with different cervical geometries.
As shown above, there is a good deal of research on the modeling of the head-neck complex for safety research. In consequence, most models have been developed to evaluate extreme conditions, such as the loading and impacts during automotive crashes. These models are useful to understand the conditions of failure in the systems and consider large deformations in the full complex, which make the use of FE appropriate.
However, no previous reports have been found in the modeling of the complex for studying the response during tracking tasks such as gaze shifts. This paper reports on the development of an anatomically correct 3D dynamic model of the head-neck complex to study and analyze the neural control of the muscles during fast eye movements. The model takes into account current knowledge about the dynamic behavior of each one of the structures forming the complex, that is, bony components, muscles, and soft tissues, with the aim of representing the dynamics of the whole complex during saccades and fast gaze shifts. In this sense, due to the small range of displacements, previous models for extreme conditions incorporating FE may not be extrapolated directly for these normal circumstances. In addition, the computational cost of using FE for a 3D dynamic simulation escalates, which creates added burden in forthcoming optimization studies on neural control of the complex [5]. A previous preliminary report considered including FE for the soft tissues that are now not in place [13].
Since vision is the most important of the senses, the oculomotor system has been the most widely addressed by researchers. Dynamic models for the oculomotor system have been proposed since 1954, with the work by Westheimer. Enderle et al. propose a linear model of the oculomotor muscle that has the static and dynamic characteristics of the rectus eye muscle within its operating range [14–16]. Recently, this linear model was parametrized for the muscles in the neck used during head movements [17], based on the nonlinear model by Zajac [18]. Simulated annealing was used to estimate the linear model’s parameter values from the virtual muscle nonlinear model, developed by Chen et al. [19] and Song et al. [20]
The ultimate goal of this work is to address the generation of the neural signals controlling the contraction of the head-neck muscles during fast eye movements. It has been suggested that the neural control of the eye muscles during saccade movements follow a time-optimal control approach [21] and a physiologically limited time-optimal control approach [14, 22].
In addition, the availability of an anatomically correct model for normal conditions allows the study and understanding of the neural paths associated with the control of this system. In particular, this model may provide information about the nature of the strategy used for the control of the head and eye responses during saccades. There are two opposing viewpoints: time-optimal control [14, 21, 23] and integrated feedback control [24, 25]. The development of this model provides an appropriate testing framework to conduct optimization analyses over these different strategies for the control of the fast eye and head movements.
The paper is organized considering first a description of the methods for the development of the morphologic representation of the head-neck complex, followed by the general framework of the model. Then, each one of the constituents of the complex is described with their corresponding modeling. Lastly, simulation results are presented, followed by discussion, with special emphasis on the neural control of the head and neck for maintaining stationary equilibrium.
2. Morphologic Representation
In order to have a physiologically correct model of the head-neck complex, a geometrical representation of the physical constituents of the system was created, using images from the visible human Project [26].
For each one of the cervical vertebrae and the skull images, contours are defined using a graphic user interface (GUI) in Matlab, shown in Figure 1. The user is asked to provide a number of points in the contour, and then, they are filtered and resampled using the forward elliptical Fourier transform [27]. The number of points to use in a representation was chosen depending on the smoothness of the structure and size. These contours are saved in a structure containing a reference to the associated bony structure.
Finally, a second GUI is used to assemble each one of the bodies from the saved contours (Figure 2). These body representations are used to determine the relative center of gravity and orientation of each one of the constituent bodies in the model. The location and size of the eyes were also verified using the Visible Human Project images. A representation of the full bony structure of the head-neck complex plus the eyes is shown in Figure 3. Finally, the attachment location of muscles, ligaments, and geometric configuration of the zygapophysial joints and intervertebral disks are determined from these images and contrasting with an anatomy atlas [28].
3. General Structure of the Model
Once the geometric description of the elements is performed, the effect of main muscles and soft tissues in the system are integrated to form the model. Figure 4 presents the general structure of the model in Simulink. To take into account the different constituents of the head-neck system, different modeling strategies are used. The skeleton (head and cervical vertebrae) are modeled as rigid bodies and muscles by means of Enderle’s dynamic model [16], with parameter values obtained through simulated annealing as reported elsewhere [17]. Finally, soft tissues (ligaments, intervertebral disk, and facet joints) are modeled using nonlinear viscoelastic representations. The use of FE for soft tissues was disregarded due to the added computational cost for running dynamic simulations and expected low deformations during physiological conditions [5]. The next sections specify the modeling approaches for each type of structure in the complex.
Each simulation is performed using the Runge Kutta method with simulation step 0.2?ms and lasting 200?ms. Dynamic states such as position, orientation, and velocities of the rigid bodies are saved. Forces generated by muscles and ligaments are saved as well to keep track of their response during the simulation. In the following subsections, specific approaches for the modeling of the different constituents of the complex are presented.
4. Bony Structures
During normal loading, bones do not significantly deform. Consequently, they are modeled as rigid bodies. Other studies for car crashes have modeled the bones as rigid bodies as well [1, 6]. Each rigid body has six degrees of freedom, that is, three translational and three rotational; its dynamics are simulated by taking into account the Newton-Euler laws by computing the net force and moment [30]. Since the creators of Matlab recently developed SimMechanics as a package for mechanical systems simulations under the Simulink approach, the implementation of the bony structures dynamics and intervertebral disks uses this approach. Each bony structure is represented as a rigid body, with the mass and inertia properties presented in Table 1. The center of gravity of the head is assumed to be 1?cm anterior and 2?cm superior to the ear [3], and the center of gravity for each vertebra is located at the middle of posterior wall of the spinal canal [3].
The initial orientation of each vertebra is measured and presented in Table 1. A body coordinate system is attached to each bony structure and located at the center of mass of each structure. A global, fixed coordinate system is assigned at the center of gravity of T2. The location of the axes for each bony element is consistent with the recommendations of the International Society of Biomechanics for the reporting of kinematic data [31]. The axis lies upward and parallel to the gravity field, the axis in the forward direction, and the axis pointing lateral to the right, completing a right-hand orthogonal triad.
Each rigid body in SimMechanics receives the net wrench (moment and force vector, [29]) from muscles, ligaments, and facet joints attaching to it. This net wrench acts at the center of gravity of the rigid body and is passed to the rigid body model by means of body actuator blocks (see Figure 5). Interactions from the intervertebral disks and the atlanto-occipital joints are modeled directly in SimMechanics.
One critical factor during the simulation of the full systems is the orientation of the bodies, for which conventional approaches use Euler/Cardanic angles. However, a number of serious issues have been raised by some authors [29, 30]. In addition to being sequence dependent and defined in a different way for each joint, parameters used in the model have singularities that affect computation of angular velocity and angular acceleration vectors. Other authors use quaternions to express the orientation and computation of the kinematics of the bodies [29, 30]. In consequence, the orientation of the skull and cervical vertebrae provided by SimMechanics are converted from rotation matrices to quaternions. The quaternion algebra is implemented using the quaternion toolbox for Matlab developed by Pierre [32].
5. Intervertebral Disks
Intervertebral discs differ from ligaments in that they respond to or experience multiple load vectors. Under any external loading, with the exception of direct uniaxial tension, disks carry compressive forces in association with other components. The existence of these other components is mainly due to the eccentric location of the head with respect to the cervical column, which imposes moments on the disks. During physiologic and traumatic load applications, cervical intervertebral discs respond to a variety of load vectors including compression, bending, and tension [9]. The disks are held in some degree of compression during normal physiological conditions due to the weight of the head [6].
The dynamics of intervertebral disks is modeled using SimMechanics as bushing joints (see Figure 6). Bushing joints allow movement in any of the prismatic or rotational axis (six degrees of freedom, DOF). The response in each one of the DOF is computed by means of an S-function. Intervertebral disk responses are linear for the three rotational axis, lateral bending and rotation, with stiffness characteristics of intervertebral disks as presented in Table 2. A nonlinear viscoelastic response in the flexion-extension DOF is implemented using Camacho et al.’s response functions [33]. Following Van Lopik and Acar implementation, the flexion extension response functions presented by Camacho et al., shown in (1), are divided by two [6], where is the response moment (in ) and is the deformation angle (in radians). The corresponding parameter values and as presented in Table 2 for both flexion and extension responses. Finally, translational damping coefficients of the disks are set to 1000?Ns/m and rotational coefficients to 1.5?Nms/rad [1, 6].
6. Ligaments
Ligaments are uniaxial structures that resist only tensile or distractive forces. However, depending on their mechanical properties (i.e., attachment points and structure), ligaments can resist tensile forces in a range of directions due to their orientation [9].
Images from the Visible Human Project [26], and data from literature [10, 34–37] were used to determine the dimensions and attachment points of the cervical spine ligaments. The following ligaments are included in the model: apical ligament, transverse ligament, alar ligaments, tectorial membrane, anterior and posterior membranes, and left and right capsular ligaments for the upper portion of the complex. To represent the dynamics of membranous ligaments, three spring elements are used.
For each motion segment (i.e., unit formed by a pair of vertebrae and their connecting structures), the following ligaments are included: anterior and posterior longitudinal ligament, flava ligament, interspinous ligament and capsular ligaments. In addition to these ligaments, included in the previous models by de Jager [1] and Van Lopik and Acar [6], the ligament nuchae is considered in the present model. Figure 7 presents a representation of the complex with the different ligaments (76 segments in total).
The ligaments are modeled as linear viscoelastic spring elements in Simulink responding only in a tensile model. The stiffness associated with each ligament is taken from Yoganandan et al. [10]. Each ligament is assumed to be at its resting length for the initial state of the model. A viscous component in the dynamics of the ligaments is added using the value of 300?Kg/s as used in previous models [1, 6].
7. Zygapophysial or Facet Joints
Like the intervertebral discs, zigapophyseal joints respond to multiple load vectors, helping them in the resistance of compressive forces. These joints provide a complementary functionality to the intervertebral discs. Because of the oblique orientation of the facet processes, the external load is resisted by normal and shear forces in the joint. These joints are fundamental for the stabilization of other tissues in the re,gion of the neck and act to limit the torsion of the disc [9]. Finite element simulations of the zigapophyseal joint include facet bone, capsular ligament, and air gap between the two cartilages. While facet bone is modeled as a solid element, the capsular ligament is modeled as explained above. The space between the two cartilages is defined using sliding or contact gap elements. The synovial fluid is modeled using fluid elements and synovial membrane using membrane elements [9].
The facet joint surfaces are approximated by slices of spheres as reported by Van Lopik and Acar [6]. The position of the facets is determined by images from visible human [26], facilitating and improving the accuracy of the process. The contact between the facets is defined and modeled as frictionless rigid body contact allowing the facets to slide relative to each other without friction. Rigid body contact reaction was modeled following de Jager [1], as (2) describes
8. Muscles Modeling
We incorporate Enderle et al.’s model for the rectus eye muscle; individual parameters are estimated and parameterized for each one of the neck muscles [16, 17]. Anatomical data about the attachment sites and force-generating parameters for the muscles are presented in [8].
In Enderle’s model, each muscle is represented as a viscoelastic parallel combination connected to a parallel combination of active state tension generator, viscosity element, and length tension elastic element. Each of the elements is linear and their existence is supported with physiological evidence [16]. The force exerted over each bone (skull or vertebrae) is the tension from the muscles attaching to it and the forces exerted by soft tissues. In case of muscles, the total length of the muscle and the shortening velocity are determined by the dynamics of the skull and the vertebrae.
Enderle’s linear model was parameterized for the muscles in the neck used during head movements based on the nonlinear model by Sierra and Enderle [17]. Simulated annealing was used to estimate the linear model’s parameter values from the Virtual Muscle nonlinear model, developed by Chen et al. [19] and Song et al. [20]. This previous report presents on the determination of the parameters of linear muscle models for the head and neck region [17].
9. Early Stopping of The Simulation
In order to avoid simulation issues due to nonphysiological states in the complex, the range of motion (ROM) of each cervical unit is considered during the computing of the response in each intervertebral disk. If the linear or angular displacement in any segment is detected to be larger than 150% of their physiological value, as previously reported [6, 9, 38], a stop signal is generated to halt the simulation and an early stopping report is saved. This feature in the model allows aborting simulations with inadequate muscle activation patterns that provide motions out of the natural ROM for the complex.
10. Simulation Results
The model is simulated first to maintain a stationary position looking in the straight-ahead direction (i.e., primary position). In order to estimate the level of muscle activation needed to maintain primary position, an optimization problem was defined. The cost function to be minimized is stated in (3) and considers the mean square error in gaze location and the mean square linear and angular displacement of each rigid body Gaze stands for the location of the point of gaze as defined below. stands for location of the center of mass of each rigid body in the model. stands for the angle of the orientation quaternion of each rigid body in the model superscript stands actual value of the parameter during the simulation. Superscript stands desired value of the parameter during the simulation. refers to the number of non fixed rigid bodies forming the complex (i.e., 9). and are adhoc selected gains to ponderate the relative components of the error. is the final simulation time.
The point of gaze (Gaze) is defined assuming that the eyes are looking straight ahead to a point located on a fixed plane board 1?m ahead of the center of the eyes. This definition of gaze point is done to be consistent with the experimental protocol to be followed in future research.
In order to penalize simulations which have been halted before the intended stopping time (), the cost function is adjusted as follows where is the desired simulation time, is the time at which the simulation was halted and is scaling factor (selected as to penalize by a factor of 3.5 a halt at 90% of simulation time)
Note that (4) does not change the Cost value if the simulation is performed without interruption for the full intended simulation time. However, in case of an early halt in the simulation due to the violation of ROM, the earlier the halt occurs, the stiffer the penalty that is introduced. An example of the penalty function is presented in Figure 8.
Using the cost gains of in (5) and the adjusted cost in (4) as objective function, the optimization problem is stated as where is the vector of active forces in the head-neck muscles.
The optimization is performed by the Genetic Algorithms (GA) using the Genetic Algorithm and Direct Search Toolbox 2.4.2 for Matlab. This routine can take advantage of the multiple cores in the current computing systems or in a cluster of computers, which speeds up the full optimization process.
GA works by taking a set of different guesses of the parameter set ( for the current problem) and computing the associated cost or fitness value. The set of guesses is called a population because once all the costs are computed (simulating each one of the scenarios and computing the cost using (3) and (4)), a selection of the best children is performed and a number of crossovers and mutations generate a new population. This new population is tested and checked again for a number of generations or until a stopping criterion is reached. For a more complete description and analysis of the metaheuristic, the reader is referred to [39].
The optimization was performed using 20 individuals per population and a maximum of 100 generations. In consequence, at the most, the algorithm would need to evaluate 2000 sets of active force values. The search was run in a Dual Core AMD Opteron Processor 285 2.61?GHz (2 Processors) with 4?GB of RAM memory and running Windows Vista 64?bits. Processing settings allowed the search to run using the four cores of the machine in parallel.
To define a good starting point with enough randomness to do a broad search, ten of the initial individuals were assigned using the stationary tension in 8 major muscles estimated by Bernhardt et al. [38] and equal activation levels for the other muscles in increases of 5% of the total active force for the muscle (from 0% to 45%). The other ten initial individual active forces were generated randomly for each muscle in the 0–100% range.
Figure 9 shows one individual used in the search. This is a set of muscle tensions to act as initial conditions to the model. Note that this individual provides very large muscle activation estimates and ends up with nonphysiological configuration of the complex. This particular simulation was aborted and a large associated cost function discouraged the usage of these estimates in the creation of offspring for further generations of the search.
(a)
(b)
After 35 generations, the best set of active forces for the muscles in the head neck complex was taken. Figure 10 shows the images of the head-neck complex for different times during the running of the simulation. Figure 11 shows the values of the tension associated with each one of the muscles in the complex. Note that for maintaining primary position, the total tension needed in each muscle is low. The bottom of this figure presents the dynamic response of ligaments during the simulation. Note that there are oscillations in the response, yet these are not strong enough as those of Figure 9.
(a)
(b)
Table 3 presents a comparison of the search time using parallel processing and serial processing. The data takes into account that in the actual search, using 4 cores in parallel, a generation of 20 individuals took on average 80 minutes to be computed. The serial processing time is estimated by taking into account the average single test run for an individual without use of parallel capabilities. This table makes clear the advantage of using parallel processing with genetic algorithms and multiple cores. The expansion of the search to be run on a cluster is transparent and will further reduce the search time.
11. Discussion
The implementation of a computer model should consider at first its intended use. This is crucial to determine the requirements associated with it and guide the decisions. Yet, there are a good number of models that use the FE modeling of the head-neck complex; they are aimed at research on passive response of the complex to car crashes [7, 9, 10, 12, 40]. Since the intended use of the model presented here is the research on the neural pathways controlling saccades and gaze shifts, the computational cost associated with FE discourages their implementation for the soft tissues [5]. In addition, the nature of normal fast orientation movements leads to the conclusion that no significant deformation takes place in soft tissues. In consequence, the implemented model considers the dynamics of the different soft tissues in the complex using a lumped parameter approach already used in recent developments for safety research [1, 6].
In addition, the bony structures are modeled using SimMechanics, which proved to be a very efficient platform for the testing of dynamic responses in mechanical systems, and these results are consistent with a previous implementation based in classical dynamics of rigid bodies [30].
A major novelty of this work is that the linear muscle models used in the head-neck region, recently determined by means of simulated annealing [17], were implemented successfully. These linear models replace previous developments using the nonlinear virtual muscle model, created by Chen et al. [19] and Song et al. [20].
The determination of the muscle activation levels for the static condition, shown in the previous section, makes use of genetic algorithms as an efficient technique when dealing with very costly optimization tasks such as the one at hand. Due to the complexity of the model, the simulation time for each set of conditions might take several minutes to evaluate. The parallel paradigm associated with the evaluation of all the member of one population (i.e., all the candidate sets) makes straightforward the implementation of parallel tasks to reduce the optimization time. One challenge remains in the determination of a feasible fitting function that addresses the problem under consideration. A second challenge is the determination of adequate initial populations considering a level of variability among them in order to take advantage of the crossovers and mutations over generations.
The proposed cost function in (3) and (4) is suitable to be used in the determination of neural pathways for the generation of saccades. For this case, gains can be set to consider only the gaze component of the cost. A second alternative is to leave a reduced component from displacement errors in the bony structures of the complex to optimize as well for minimum variability in the complex.
The major difference in this case is the set to optimize. Saccades and gaze shifts can be generated by pulse-step profiles of activation [15, 41, 42]. Consequently, the set of parameters to optimize change to be the times of the steps and the activation levels for muscles associated with the movement.
12. Conclusions
In this paper, the development of the 3D dynamic model of the head neck complex for research on fast orientation movements is presented. Specific details about the implementation, ranging from the determination of morphologically correct locations for the structures to their dynamic are discussed. An assessment of the condition of the model components to avoid states outside of the physiological range of movement is considered as well.
In order to make the model useful for research of neural control of the complex during fast orientation movements, a global optimization problem is stated. The implementation of the solution using genetic algorithms is discussed; also, the advantages of using parallel processing to speed up the solution of the problem are highlighted. These efforts leave a working model suitable for performing the evaluation of different strategies of neural control.